Mechanistic Modelling of Biomass Growth, Glucose Consumption and Ethanol Production by Kluyveromyces marxianus in Batch Fermentation

This paper presents results concerning mechanistic modeling to describe the dynamics and interactions between biomass growth, glucose consumption and ethanol production in batch culture fermentation by Kluyveromyces marxianus (K. marxianus). The mathematical model was formulated based on the biological assumptions underlying each variable and is given by a set of three coupled nonlinear first-order Ordinary Differential Equations. The model has ten parameters, and their values were fitted from the experimental data of 17 K. marxianus strains by means of a computational algorithm design in Matlab. The latter allowed us to determine that seven of these parameters share the same value among all the strains, while three parameters concerning biomass maximum growth rate, and ethanol production due to biomass and glucose had specific values for each strain. These values are presented with their corresponding standard error and 95% confidence interval. The goodness of fit of our system was evaluated both qualitatively by in silico experimentation and quantitative by means of the coefficient of determination and the Akaike Information Criterion. Results regarding the fitting capabilities were compared with the classic model given by the logistic, Pirt, and Luedeking–Piret Equations. Further, nonlinear theories were applied to investigate local and global dynamics of the system, the Localization of Compact Invariant Sets Method was applied to determine the so-called localizing domain, i.e., lower and upper bounds for each variable; whilst Lyapunov’s stability theories allowed to establish sufficient conditions to ensure asymptotic stability in the nonnegative octant, i.e., R+,03. Finally, the predictive ability of our mechanistic model was explored through several numerical simulations with expected results according to microbiology literature on batch fermentation.


Introduction
Alcoholic fermentation is an anaerobic process that transforms sugars like glucose or fructose into ethanol and carbon dioxide. Several yeast species are used commonly in this process, e.g., Kloeckera, Hanseniaspora, Candida, Pichia, Kluyveromyces, and Saccharomyces among others. The growth rate of these microorganisms has an ultimate effect on the sensorial characteristics of the final product, which can be positive or negative depending on the yeast used [1].
Overall, yeasts are indispensable for biotechnological processes such as wine and beer production [2]. In this research, we focus on investigating glucose consumption and ethanol production from several strains of Kluyveromyces marxianus (K. marxianus). This yeast has a great potential for alcoholic fermentation due to its intraspecific characteristics such as higher specific growth rates, the ability to grow on various substrates, and tolerance to high temperatures [3][4][5]. Further, Kluyveromyces sp. produces aromatic compounds such as fruity esters, carboxylic acids, ketones, furans, and other alcohols in liquid fermentation such as 2-phenyl ethanol whose sensorial characteristics can influence the quality of wine, distilled drinks, and fermented foods [6]; refer to Fonseca et al. for an extensive review on the biotechnological potentials of K. marxianus [4,6].
Concerning industrial applications, fermentation is commonly performed in batch culture, which brings certain advantages such as the reduction of contamination risk, in addition to the fact that a large capital investment is not necessary since high-priced production equipment is not required compared to a continuous culture process [7]. Batch process implies that yeasts are incubated in a closed container under controlled conditions with a culture medium composed of the necessary nutrients [8]. Hence, biomass cannot grow indefinitely and four phases have been identified in its dynamics, i.e., lag phase, exponential growth, stationary state, and death phase. While this process is carried out, the substrate is consumed and converted into the product, e.g., ethanol produced by sugars such as glucose [9]. Therefore, properly identifying the time interval of these phases as well as predicting the maximum product concentration that could be produced from the initial concentrations of both substrate and biomass may help to optimize production costs on the resulting product of several applications. The latter may be achieved by mechanistic modeling through predictive microbiology, which can be considered a powerful tool to investigate and summarize the overall effects of varying conditions and environmental factors within food formulation and processing [10]. Further, mathematical models could aid in gaining insights concerning microbial food safety and quality assurance of increasingly complex food products [11,12], as well as estimating shelf life and forecasting food spoilage [13,14].
Mathematical models in predictive microbiology can be classified according to different criteria, uses, and functionalities that are not mutually exclusive. Based on the type and number of variables, models are classified into primary, secondary, and tertiary; they can also be differentiated on the basis of their mathematical background as mechanistic or empirical [15], and they can be categorized into structured and unstructured conforming to the complexity of the chemical compounds of the biomass [16]. Primary models are those that represent biomass growth dynamics as a function of time, the main equations in the literature are the exponential functions of Gompertz [17] and Vazquez-Murado [18], the logarithmic function of Baranyi et al. [19] and the cubic model of Garcia et al. [20]. All models are described by parameters such as maximum growth rate [µ max ], lag time [L], and both initial [X 0 ] and maximum biomass [X max ] Concentrations, while secondary models relate to the latter with environmental conditions such as temperature and pH, and other variables such as substrate and product concentrations over time, e.g., equations of Monod [21], Teisser [22], Haldane [23], and Moser [24], which aim to describe biomass growth dynamics as a function of substrate concentration and have been widely used to investigate bacterial growth [25]. Tertiary models are the result of combining primary and secondary models through the use of computer tools that allow predictions regarding the growth or death of microorganisms in food when different environmental conditions are combined [26]. Concerning the second classification mentioned, mechanistic models are formulated by means of theoretical bases and provide an interpretation of microbial growth in terms of known processes and empirical models are usually composed of polynomials of the first or second degree and pragmatically describe the data with convenient mathematical relationships, this does not usually give information on precise responses of microorganisms, because they do not take into account known processes [27]. Finally, according to the third category described, unstructured models consider biomass only as a chemical compound in a culture and its dynamics is described by simple models, while structured models also take into account changes in the internal cellular structure of biomass in terms such as the content of RNA, enzymes, reagents, metabolism and products [28]. The Gompertz, Vazquez-Murado, Baranyi and Garcia models, mentioned above, are also classified as unstructured models since biomass is considered a variable described only by its concentration. Mathematical models used by Sansonetti [29], Lei [30], and Steinmeyer [31] are classified as structured because they describe the growth of biomass considering the intracellular reactions produced by its metabolism.
Thus, it is important to highlight that in a batch fermentation process, multiple reactions occur, so the adaptability and evolution of microorganisms in short periods and changes in environmental conditions usually characterize this type of process, consequently, the modeling of these systems is complex due to time-varying characteristics of biological systems, resulting in nonlinear systems dynamics [28]. Hence, a mathematical model formulated from a system of nonlinear differential equations will allow the application of nonlinear systems control methods to optimize the process so that the characteristics of the final product can be predicted when the environmental conditions of the culture are controlled and the initial conditions of biomass, substrate, and product values are known. It is worth mentioning that most of the models found in the literature focus on the yeast Saccharomyces cerevisiae since it is one of the most used in the industry; however, biotechnological opportunities have been found in non-Saccharomyces yeasts because they have metabolic characteristics that lead to the production of compounds of interest. Therefore, it is important to model the growth of K. marxianus because of the great potential in the production of esters compounds of industrial importance [32]. Thus, in this paper, we applied mechanistic and computational modeling to formulate a system of three coupled nonlinear first-order Ordinary Differential Equations (ODEs) that describe dynamics between biomass, glucose (substrate), and ethanol (product) concentrations over time.
Mechanistic modeling allowed us to provide both qualitative and quantitative descriptions concerning the relationships of biomass growth, glucose consumption, and ethanol production from 17 strains of K. marxianus, while computational modeling was used to fit experimental data from these three variables and establish numerical values for each parameter of the mathematical model. Further, nonlinear theories such as the Localization of Compact Invariant Sets (LCIS) method and Lyapunov's Stability Theory were applied to provide a complete analysis of the local and global dynamics of our proposed biological system [33].

Materials and Methods
This section provides all the information concerning the experimental data of biomass growth, substrate consumption and ethanol production, i.e., karyotypes of the K. marxianus strains with identifiable chromosomal differences among them, environmental conditions, chemical characteristics of the medium, lab equipment used for measurements, and periods for each measurement, then the mathematical model is formulated and each equation as well as values and units of parameters are described. This section concludes by describing the procedure to approximate the experimental data and to fit the numerical values of each parameter by designing an algorithm in Matlab.

Experimental Data: Culture Medium and Analytical Techniques
Experimental data was obtained from alcoholic fermentation in batch culture by K. marxianus, 17 strains with different genetic profiles were incubated in 20 g/L of yeast extract peptone dextrose agar at 30 • C in order to study their kinetic growth, glucose consumption and ethanol production. Codification and origin of studied karyotypes of K. marxianus are identified by Páez et al. in [34], where 15 strains were obtained from different places of México, and they were isolated from agave fermentation for mezcal production, in addition to 2 reference strains that were isolated from pozol (CBS6556) in México, and from yoghurt (CBS397) in Netherlands.
Biomass concentration was measured with a spectrophotometer UV-VIS DR 6000 (HACH, Loveland, CO, USA) by optical density at 600 nm, values in g/L were obtained relating optical density with a calibration curve of the dry weight of K. marxianus. For glucose consumption and ethanol production by High-Performance Liquid Chromatography (HPLC series 1200, Agilent Technology, Palo Alto, CA, USA), a BIORAD HP-87H + (8%) ion exchange column was used, in an AGILENT® 1200 series equipment, with H 2 SO 4 0.005 N as mobile phase, at a flow of 0.5 mL/min, the column temperature was 60 • C, and the Refractive Index detector temperature was 60 • C. The injection volume of 5 µL, calibration curves were made with glucose and ethanol Sigma Aldrich at 99% purity or higher, and a determination coefficient higher than 0.99 for each compound [36,37].
Fermentation was made in duplicate for every strain and samples were taken each hour for 13 consecutive hours. Two samples were taken every hour for each variable in the time interval of the process where t goes from 0 to 13, then the average value of the two measurements was computed. Therefore, each variable, i.e., biomass [x(t)], glucose [y(t)], and ethanol [z(t)], has 14 observations with a total of 42 experimental data points (n) for each K. marxianus strain.

The KM Mechanistic Model
The KM mechanistic model is proposed to describe the dynamics of alcoholic fermentation. This is a biochemical process carried out by yeasts (also known as biomass), to transform sugars such as glucose into ethyl alcohol, otherwise known as ethanol (main product) and other byproducts. In this case, the alcoholic fermentation is taken in a batch fermentation process with established laboratory conditions of temperature and known initial glucose concentrations (substrate). Our mathematical model describes the relationships between biomass concentration [x(t)], glucose consumption [y(t)], and ethanol production [z(t)] over time. The set of three first-order ODEs is presented beloẇ where each state variable x(t), y(t) and z(t) are measured in g/L, and the time unit is given in hours. Now, by considering results from Leenheer and Aeyels (see Section II.A in [38]), all solutions with nonnegative initial conditions [x(0), y(0), z(0) ≥ 0, ] will be located in the nonnegative octant as indicated below i.e., each positive half trajectory of the system will be positively forward invariant in R 3 +,0 . The latter also considers the biological sense of each variable as there is no meaning for negative values of biomass, glucose or ethanol in the scope of the KM system (1)-(3). It is important to mention that variables cannot grow exponentially indefinitely, and they must have biologically feasible limits which will be discussed in the next section. Values and units of each parameter of the KM system (1)-(3) are given in Table 1. Now, let us describe our mechanistic model based on the experimental data described in the previous section and the following biological assumptions. Biomass growth dynamics is described by Equation (1), where the first term uses the classical Monod form for the growth of microorganisms [21], where ρ 1 is the biomass maximum growth rate (also found in the literature as µ max ), and ρ 2 is the affinity or half-velocity constant for glucose consumption. The second term describes biomass death due to ethanol accumulation toxicity by the law of mass action (see Section 2.3 in [39]) with a rate ρ 3 . This term is negative because ethanol accumulation increases the membrane fluidity and negatively affects the membrane protein's function, which can lead to cell growth inhibition or even death [40,41]. The third term represents the natural yeast death rate [ρ 4 ] mainly due to environmental resources depletion [42]. Glucose dynamics is formulated in Equation (2) as a decrescent function where the law of mass action gives the first two terms. The first one represents glucose consumption to support biomass growth. In contrast, the second term accounts for the glucose consumption used for ethanol production, with rates ρ 5 and ρ 6 , respectively. The third term represents the spontaneous decomposition rate of glucose [ρ 7 ] [43]. The latter stems from the fact that the culture medium is placed in a sealed container in batch fermentation, and no other nutrients (primarily glucose) are supplied into the system. Ethanol dynamics is described in Equation (3). The first term represents ethanol production associated with biomass growth. It is due to ethanol being recognized as a primary metabolite, a product obtained from reactions required or cellular growth [44,45]. The second term represents the glucose conversion to ethanol not directly linked to cellular growth, attributed to the need for Gibbs's free energy for cellular maintenance [44,46]. In both cases, terms are formulated by the law of mass action with respective rates ρ 8 and ρ 9 . Finally, the third term represents ethanol degradation with a rate ρ 10 . The flow diagram shown in Figure 1 was constructed to illustrate the dynamics of the system. It should be noted that fixed parameters values were estimated for the 17 K. marxianus strains, particularly for ρ 2 , ρ 3 , ρ 4 , ρ 5 , ρ 6 , ρ 7 and ρ 10 . Further, concerning the death rate of biomass [ρ 4 ], spontaneous decomposition rate of glucose [ρ 7 ], and degradation rate of ethanol [ρ 10 ], one can see that they are in a different order of magnitude and the following constraint is formulated for these three parameters: Now, concerning the equilibrium points of the KM system, Equations (1)-(3) have a unique biologically meaningful equilibrium point in the domain R 3 +,0 given by Another set of five equilibria with at least one negative value is shown in Appendix A; therefore, these equilibrium points are discarded from the biologically meaningful dynamics of the system. Further, from the biological characteristics of each variable the following can be stated with respect to each solution as time increases lim t→∞ x(t) = lim t→∞ y(t) = lim t→∞ z(t) = 0, due to the eventual death of microorganisms, glucose consumption and ethanol degradation [47], asymptotic stability of trajectories is discussed in the next section.

Parameter Value Estimation
First, let us compute the glucose decomposition rate [ρ 7 ] by assuming a first-order kinetics [48] for glucose dynamics, and considering a half-life [t 1/2 ] of 96 years [43]. Then, ρ 7 can be computed from the next equation where y 0 is the glucose initial concentration, i.e., y(0); hence Now, in order to determine the numerical values of parameters ρ i , i = 1, 2, 3, 4, 5, 6, 8, 9, 10; the computational model of Equations (1)-(3) was formulated as follows by applying Euler's method (see Section 1.7 in [49]) where ∆ t was set to 1 × 10 −5 . Then, an algorithm was formulated in Matlab 2022b with the lsqcurvefit function from the optimization toolbox as its core [50] (initial points were set as 1 × 10 −1 for each parameter [ρ 1 , ρ 8 , , and optimotions of the function were set as follows: Max Function Evaluations = 1 × 10 3 , Max Iterations = 1 × 10 3 , and Function Tolerance = 1 × 10 −9 ). This allowed us to establish a fixed value for parameters ρ j , j = 2, 3, 4, 5, 6, 10; by averaging the corresponding values for each of the 17 strains, these results are shown in Table 1. However, this procedure was not applicable for parameters ρ 1 , ρ 8 and ρ 9 , as it was expected that each strain of K. marxianus will have its own biomass growth rate [ρ 1 ], and its corresponding ethanol production rates (ρ 8 and ρ 9 ), this is directly linked to the chromosomal differences among the strains affecting their growth kinetics. Hence, the main algorithm was redesigned to fit these three parameters and consider the others as fixed constants. Overall results are shown in Table 2 with their corresponding standard error (SE), and 95% confidence interval (CI). These two statistics allow us to establish that estimates for parameters ρ 1 , ρ 8 , and ρ 9 in the 17 strains are statistically significant. The latter follows from the fact that each SE(ρ k ) < ρ k /2, k = 1, 8, 9; i.e., the value of the SE is less than half of the value fitted for each parameter, thus, the null hypothesis [ρ k = 0] can be rejected (see Section 5.2.8 from Koutsoyiannis [51]). Furthermore, both the lower and upper limit of the 95% CI of all fitted values are positive, hence, as there is no change in the sign of the bounds, this implies that the value of the null hypothesis is excluded, and one can conclude that all P-values are less than 0.05 (see Chapter 17 from Motulsky [52]). Strain  Finally, it should be noted that the in silico experimentation performed in this research was done on a high-end desktop computer with a Ryzen 9 5950X CPU, 128 GB of RAM DDR4 CL18, a 12 GB GPU NVIDIA GeForce RTX 3080, and 1 TB Samsung 980 Pro Gen 4 NVMe M.2. The complete algorithm that was designed to fit the numerical values of parameters [ρ 1 , ρ 8 , ρ 9 ], and to determine results concerning the statistics and goodness of fit can be found in the Supplementary Materials.

Results
In this section, the in silico experimentation is performed by means of several numerical simulations, and results relating to the nonlinear analysis of the system are derived, i.e., bounds for the localizing domain, asymptotic stability, and existence and uniqueness for all solutions of our model in the nonnegative octant R 3 +,0 .

In Silico Experimentation and Goodness of Fit
First, qualitative results are illustrated by means of numerical simulations. For the sake of simplicity, the strains were clustered in groups of four from strain 1 to the 16 (see respectively), and results concerning only for the strain 17 are shown in Figure 6. In all panels, the × green marker represents the average value for the two experimental data measurements for each variable, i.e., biomass [x(t)], glucose [y(t)], and ethanol [z(t)], whilst the blue continuous line represents the approximated value given by the KM system (1)-(3) when is solved by means of Equations (6)-(8) with ∆ t = 1 × 10 −5 . The time units are given in hours and the concentration for each variable is measured in g/L as indicated in each axis. Values for all ten parameters corresponding to each strain are shown in Tables 1 and 2. Now, let us provide a quantitative measure of the fitting capabilities of the KM mechanistic model (1)-(3), thus, the coefficient of determination R 2 is calculated for each variable with results shown in Table 3.       Further, the Akaike Information Criterion (AIC) [53][54][55] was computed by considering a small sample relative to the number of parameters (n/K < 40) with a bias correction as indicated below where n is the total number of experimental data points; f e the experimental data and f a the approximated value for the residual sum of squares (RSS); and K the number of parameters of the system; therefore, n/K = (3 × 14)/10 = 4.2. Results, including RSS, AIC and R 2 , for the complete trajectory of the system, i.e., φ(x, y, z) for the total of 42 experimental points (14 for each variable) are summarized in Table 4. Table 4. In order to provide overall measures for the fitting capabilities of our mathematical model, i.e., the KM system (1)- (3), values were calculated for the RSS, the AIC, and the R 2 to estimate and describe the dynamics between the three variables φ(x, y, z), where x(t), y(t) and z(t) represent, respectively, the evolution of biomass, glucose, and ethanol.  The AIC yields a value that relates the amount of information that our model loses when approximating the experimental data. Hence, one can compare the capabilities of the model to estimate the concentrations over time of biomass [x(t)], glucose [y(t)], and ethanol [z(t)] among the 17 K. marxianus strains while providing a statistical measured for the quality of the KM system (1)-(3).

Nonlinear Analysis: Localizing Domain, Asymptotic Stability, Existence and Uniqueness
The localizing domain can be determined by computing the upper bounds for all variables of the KM mechanistic model (1)-(3), the lower bounds are given by the boundary of the domain R 3 +,0 , i.e., {x inf = 0, y inf = 0, z inf = 0}. The latter is achieved by means of integration and the LCIS method [56]. Within the localizing domain, one may find all biologically meaningful dynamics of the system, i.e., compact invariant sets such as equilibrium points, periodic orbits, limit cycles and chaotic attractors (see Section 3 in [57]), among others.
First, in order to find the upper bound for the glucose concentration [y(t)], Equation (2) is integrated as follows dy by considering ρ 7 > 0 and x(t), z(t) ≥ 0 from the domain R 3 +,0 . Then, with y(0) ∈ R +,0 . Therefore, all solutions with nonnegative initial conditions will be bounded as indicated below hence, any upper bound for x(t) and z(t) depending on K y will be directly related to the glucose initial concentration [y(0)], which is expected as biomass and ethanol production over time is directly related to glucose dynamics. Now, let us provide the mathematical background that allows us to compute a localizing domain where all compact invariant sets of a nonlinear dynamical system are located. The General Theorem concerning the LCIS method was formalized by Krishchenko and Starkov (see Section 2 in [58]) and it states the following: Each compact invariant set Γ oḟ x = f (x) is contained in the localizing domain: From the latter f (x) is a C ∞ −differentiable vector function where x ∈ R n is the state vector. h(x) : R n → R is a C ∞ −differentiable function called localizing function, h| S denotes the restriction of h(x) on a set S ⊂ R n with S(h) = x ∈ R n | L f h(x) = 0 , and L f h(x) = (∂h/∂x) f (x) is the Lie derivative of f (x).
Hence, one can define h inf = inf{h(x) | x ∈ S(h)} and h sup = sup{h(x) | x ∈ S(h)}. Furthermore, if all compact invariant sets are contained in the set K(h i ) and in the set K h j then they are contained in K(h i ) ∩K h j as well. The nonexistence of compact invariant sets can be considered for a given set Λ ⊂ R n if Λ ∩ K(h) = ∅, then the systemẋ = f (x) has no compact invariant sets located in Λ.
Following the LCIS method, one can explore the next localizing function then, the Lie derivative may be written as follows and the set S(h 1 ) = L f h 1 = 0 is given by and the next two conditions are formulated where (9) is directly fulfilled by (4). Now, let us apply the Iterative Theorem in order to find an upper bound for the localizing function from the latter, the upper bound for the biomass concentration [x(t)] may be written in terms of the parameters and the initial glucose concentration [y(0)] as follows Now, an upper bound for the ethanol concentration [z(t)] can be determined by the following localizing function h 2 = β 1 x + β 2 y + z; β 1 , β 2 > 0, whose Lie derivative is computed as indicated below and at this step, the following conditions are formulated then, set S(h 2 ) = L f h 2 = 0 , can be written as follows where the next condition is formulated and it holds by (4). Then, the Iterative Theorem is applied to get the following result then, the upper bound for the localizing function h 2 is derived as follows now, from the latter one can get the upper bound for ethanol concentration [z(t)] over time in terms of the parameters, the initial glucose concentration [y(0)] and the upper bound of biomass x sup as given below Results shown above allow us to conclude the following regarding the boundedness of the KM system (1) Now, let us briefly provide the mathematical background concerning the stability theory in the sense of Lyapunov, particularly the direct method where it is necessary to formulate a Lyapunov candidate function, which is usually denoted as V(x) : R n → R, a continuously differentiable function whose temporal derivative is given bẏ V(x) = (∂V/∂x) f (x). This function must be positive definite, i.e., V(0) = 0 and V(x) > 0 for x = 0, whilst a negative definite function is also V(0) = 0 but V(x) < 0 for x = 0. Further, function V(x) is said to be radially unbounded if V(x) → ∞ as x → ∞. The latter allows the formulation of the Global Asymptotic Stability Theorem (see Chapter 4 in [59] and Chapter 2 in [60]) which states the following: The equilibrium point x * is globally asymptotically stable if there exists a function V(x) positive definite, radially unbounded and decrescent such that its temporal derivativeV(x) is negative definite. A function V(x) satisfying properties of this theorem is called Lyapunov function.
Following the latter, let us explore the next Lyapunov candidate function then, the time derivative is computed as shown beloẇ which can be rewritten as followṡ where it is evident thatV(0, 0, 0) = 0, therefore the following constraints on coefficients γ 1 and γ 2 are formulated to ensureV(x, y, z) < 0 thus, as parameters ρ i , i = 1, 2, 3, 5, 6, 8, 9; in both conditions are different for each term, then it is possible to assume that there exists a set of solutions that satisfies (14) and (15). Hence, the following result can be concluded: (14) and (15) are fulfilled, then the KM mechanistic model (1)- (3) is asymptotically stable and all trajectories will go to the equilibrium point (x * 0 , y * 0 , z * 0 ) = (0, 0, 0).
Concerning the existence and uniqueness of solutions for the KM system (1)-(3), let us introduce the following notations for the sake of simplicity x, y, z) = ρ 8 xz + ρ 9 yz − ρ 10 z, and compute the Jacobian matrix [∂ f /∂u](t, u) (see [49] at Section 7.4) with results shown below for f i (t, u), i = 1, 2, 3; and and it is evident that f i (t, u) and [∂ f /∂u](t, u) are continuous and exist on the domain [33]. Hence, the latter implies that f i (t, u) is locally Lipschitz in u on Ω (see Lemma 3.2 by Khalil in [59] at Section 3.1). Further, each element of (16) is bounded by Theorem 1. Thus, the following can be concluded:

Discussion
The KM mechanistic model (1)-(3) was formulated by considering the biological relationships between each variable in a controlled batch fermentation where concentrations in g/L were measured for biomass [x(t)], glucose [y(t)], and ethanol [z(t)] over 13 consecutive hours. Then, by means of the lsqcurvefit function, an algorithm was developed in Matlab to approximate the experimental data from the 17 K. marxianus strains discussed at Section 2; both qualitative (see  and quantitative (see Tables 3 and 4) results were shown in Section 3. The in silico experimentation illustrates the capabilities of the system to approximate the experimental data of each strain, whilst both the R 2 and the AIC provide a value for the goodness of fit of the model to each set of data. In Table 4, one can see that R 2 values range from 0.955 to 0.994, and AIC from −43.478 to 33.184, these values are for strains 7 and 9, respectively. Now, it should be noted that the dynamics between biomass growth, substrate consumption and product generation have been modeled before by means of the logistic growth law [62], the Pirt Equation [63], and Luedeking-Piret Equation [64] as indicated below in Equations (17)- (19), respectively: where µ max is the biomass maximum growth rate, this parameter is equivalent to ρ 1 in our mathematical model; X max the maximum concentration value of biomass in the experimental data set for the time-interval of the process being observed; Y X/S the biomass/substrate yield; m is the maintenance coefficient; α is the growth-associated coefficient for the product; and β is the non-growth-associated coefficient for the product. Our algorithm was applied to approximate the experimental data of the 17 K. marxianus strains with overall results shown in Table 5. The main comparison between the KM system (1)-(3) and Equations (17)- (19) is performed with respect to the biomass maximum growth rate, given by ρ 1 and µ max , respectively. Tables 2 and 5 show that estimated values of ρ 1 are on average ∼0.717 smaller than those estimated for µ max . The latter is a direct consequence of the biological assumptions on which each mechanistic model was formulated. The KM system (1)-(3) was constructed by considering interactions between the three variables as illustrated in the flow diagram of Figure 1, whilst the logistic, Pirt, and Luedeking-Piret Equations (17)- (19) are constructed by only assuming a logistic growth for biomass without taking into account the overall effect of ethanol production over the entire system as well as the death rate of biomass [x(t)], decomposition rate of glucose [y(t)], and degradation of ethanol [z(t)]. Further, the in silico experimentation concerning Equations (17)- (19) illustrated in Figures A1-A5 at Appendix B shows that approximated values for substrate [S(t)], i.e., glucose, becomes negative as time increases, which is not biologically possible for this variable. Further, one can see from the experimental data that ethanol production does not follow a smooth sigmoidal growth, the data even illustrates degradation among some strains, which is better approximated by our model as it is shown in the lower panels of Figures 2-6.
When comparing the goodness of fit by computing the AIC and R 2 , it is evident that the KM system (1)-(3) had overall better results than the logistic, Pirt, and Luedeking-Piret Equations (17)- (19). Although the latter has fewer parameters than ours (six and ten, respectively) and the AIC penalizes a model with more parameters to be fitted, results for the RSS were lower for the KM system which ultimately worked in our favor. Further, the capabilities of the KM mechanistic model may extend beyond its ability to approximate experimental data and estimate the biomass maximum growth rate, in Appendix C the in silico experimentation illustrates the dynamics for t ∈ [0, 39], i.e., three times the period for the experimental data. Figures A6-A10 show that as time increases and the substrate is no longer added into the system, then the death of biomass and degradation of ethanol begins to take over the system. The latter was expected from the asymptotic stability results of Section 3, particularly Theorem 2 and Corollary 1, as these state that the concentration of all variables will eventually be zero, i.e., both biomass [x(t)] and ethanol [z(t)] concentrations are going to be depleted. Additionally, it is important to note that all solutions of the KM system are bounded from above, which is consistent with the localizing domain results of Theorem 1.
Regarding the values of parameters m and β, our algorithm yielded results in the magnitude of 10 −14 for m in all strains; in fact, setting m to zero does not affect the ultimate results for the other parameters [µ max , Y X/S , α, and β] which may allow us to completely disregard this term [−mX] from Equations (17)- (19). Concerning β, values for 12 strains were in the same order of magnitude 10 −14 , however, the following results were determined for strains 3, 4, 11, 12, and 17: 87.633 × 10 −3 , 38.660 × 10 −3 , 51.050 × 10 −3 , 53.813 × 10 −3 , 37.702 × 10 −3 , respectively. Hence, the non-growth-associated coefficient for the product may influence the dynamics in some karyotypes of K. marxianus.

Conclusions
Mechanistic modeling has proven to be a powerful tool capable of describing the relationships between different variables in the dynamics of biological systems when considering assumptions based on scientific principles underlying the phenomenon being modeled. In this work, a set of three coupled first-order ODEs was formulated which can approximate experimental changes over time of alcoholic fermentation in batch culture by 17 different strains of K. marxianus. The KM mechanistic model (1)-(3) describes biomass growth [x(t)], glucose consumption [y(t)], and ethanol production [z(t)] in concentrations of g/L per hour. The parameter values of the system were estimated through a nonlinear curve-fitting algorithm in Matlab with the experimental data of each batch culture fermentation described in Section 2. The latter allowed us to conclude that seven parameters have the same numerical value for the dynamics observed in the 17 strains, particularly the affinity with substrate constant [ρ 2 ], inhibition rate of biomass growth due to product accumulation [ρ 3 ], biomass death rate [ρ 4 ], consumption rates for biomass growth and ethanol production [ρ 5 and ρ 6 ], glucose spontaneous decomposition rate [ρ 7 ], and ethanol degradation rate [ρ 10 ]; these values are shown in Table 1. However, the biomass maximum growth rate [ρ 1 ], ethanol production associated with biomass growth [ρ 8 ], and glucose converted in ethanol [ρ 9 ] parameters have specific values for each strain, results are shown in Table 2 with a 95% confidence interval that gives us the margin of error for each parameter value estimation.
As predictive microbiology establishes, mathematical models must be simplified until measurable parameters can be obtained, the KM mechanistic model successfully achieves this with ρ 1 , ρ 8 , and ρ 9 as the main parameters that describe the overall dynamics of the batch fermentation process under study in this research. The biomass growth rate is a very specific value for each strain that must be as high as possible. Ethanol production with respect to biomass growth represents the fermentative capacity of each strain, and the concentration of glucose converted to ethanol is directly related to these rates. It should be noted that in batch culture the latter requires high sugar concentrations to achieve alcoholic fermentation.
Further, the in silico experimentation illustrates that our model may be able to accurately predict the concentration of each variable as it is shown in Appendix C; nonetheless, further experimental data are needed to properly validate this assessment. One can see in Figures A6-A10 that when no more substrate is added to the culture, then biomass growth goes into the death phase, and ethanol degradation begins to happen in the system. This behavior is to be expected as the nonlinear analysis of the system allowed us to conclude that all concentrations will eventually go to zero in the absence of glucose, i.e., the asymptotic stability of the equilibrium (5) [(x * 0 , y * 0 , z * 0 ) = (0, 0, 0)] by Theorem 2 and Corollary 1; further, concentrations over time of all variables are bounded by the Localizing Domain Theorem 1. The latter is illustrated in all panels for the predictions of biomass growth [x(t)], glucose consumption [y(t)], and ethanol production [z(t)].

Appendix C. Predictive Ability of the KM Mechanistic Model
This appendix presents results concerning the in silico experimentation when solving the KM mechanistic model (1)-(3) for a time interval of t ∈ [0, 39] in order to illustrate its ability to predict the dynamics of the three variables, i.e., the concentration in g/L over time between biomass [x(t)], glucose [y(t)], and ethanol [z(t)], after the last experimental data point without further substrate addition into the batch. It should be noted that at this stage of the research there is not available data to validate if the model is be able to accurately predict the evolution of both biomass and ethanol in the system, further experimental data points could be helpful to better fit the values of biomass death rate [ρ 4 ], and ethanol degradation rate [ρ 10 ]. However, Figures A6-A10 allow us to illustrate results concerning Theorem 2 and Corollary 1 from the asymptotic stability analysis. x(t) [g=L] x(t) [g=L] x(t) [g=L] x(t) [g=L]  x(t) [g=L] x(t) [g=L] x(t) [g=L] x(t) [g=L]