Coupled Logistic Models to Describe Culture pH, Growth, Product Formation, and Substrate Consumption in Batch Cultures of Lactic Acid Bacteria

1. Introduction

Lactic acid bacteria (LAB) are widely used in food fermentations and biotechnological applications, where the efficiency of biomass production and substrate utilization directly affects product yield and process optimization [1,2,3].

These bacteria produce organic acids (mainly lactic acid) through the catabolism of sugars present in the substrate, leading to a decrease in culture pH [3,4]. Biomass accumulation critically depends on culture pH, as it influences cellular metabolism, intracellular homeostasis, enzyme activity, and nutrient uptake [4,5,6,7]. Previous studies have shown that biomass production reaches a maximum within an optimal pH range (e.g., near neutral or mildly acidic conditions) for many LAB strains. In contrast, excessively low pH values significantly reduce biomass accumulation and may even lead to cell death, as acid stress inhibits metabolism and cell division [8,9,10] and can impair LAB growth more severely than nutrient depletion [4,11].

Acid accumulation increases the proportion of undissociated acid species, which can diffuse back into the cell, dissipate the proton motive force, and drain cellular energy as bacteria expel excess protons—ultimately inhibiting growth and reducing biomass yield [12,13,14,15]. Thus, in LAB fermentations, biomass formation is directly linked to the availability of growth-limiting nutrients, particularly nitrogen sources, and is strongly modulated by culture acidification resulting from lactic acid production [4,16,17].

LAB are also capable of producing bacteriocins—antibacterial peptides often considered pH-dependent primary metabolites—which are highly sensitive to the time course of culture pH during fermentation [18]. Therefore, changes in culture pH can modulate bacteriocin synthesis [19,20], likely through regulation of biosynthetic genes [21,22], modulation of the enzymatic activity responsible for post-translational processing of prebacteriocins into their active forms [15,23], and adsorption of bacteriocins to the cell wall of the producer strain [24].

Furthermore, total sugar and nitrogen consumption are governed by biomass production, as these nutrients serve as energy sources to support growth and cellular maintenance, are incorporated into new biomass, and are utilized for metabolite production [17,25,26].

Mathematical modeling of LAB fermentations is an essential tool for understanding, predicting, and optimizing microbial growth, metabolite production, and substrate utilization in complex culture systems [25]. Among the available empirical approaches, logistic-type equations have been widely and successfully applied to describe batch biomass growth, lactic acid formation, antimicrobial metabolite production, substrate consumption, and pH evolution during LAB fermentations due to their simplicity and good fitting performance across different media and strains [27,28,29,30]. These models assume self-limiting behavior and capture the sigmoidal dynamics typically observed in batch cultures [27,28].

However, a major limitation of these equations is that they describe culture variables solely as functions of time, implicitly assuming that microbial growth, product formation, and nutrient consumption in batch cultures are purely autocatalytic processes. Moreover, such models do not explicitly account for the underlying biological and biochemical mechanisms driving growth, metabolite production, and nutrient utilization in LAB fermentations.

Consequently, treating kinetic parameters—such as growth, production, and consumption rates—as independent constants may limit the predictive power and physiological relevance of these models.

By linking the kinetic parameters of the logistic equations through biologically meaningful relationships—accounting for pH effects, nutrient limitation, and growth–production coupling—the proposed approach provides a more precise, comprehensive, and integrated description of LAB fermentation dynamics while retaining the mathematical simplicity of logistic models. Such models not only enhance our understanding of microbial physiology but also provide practical benefits for the design and control of fermentation processes. This approach may improve the predictive accuracy of growth models and support the optimization of culture conditions, ultimately contributing to more efficient fermentation processes and a deeper understanding of microbial physiology.

To address the limitations of logistic-type equations in describing the kinetics of LAB fermentation, the present study focuses on developing an improved modeling framework based on coupled logistic models that incorporate the relationships among biomass growth, lactic acid production, total antibacterial activity, pH evolution, and nutrient consumption (total sugars and nitrogen).

2. Materials and Methods

2.1. Bacterial Strains

The bacterial strains used in this study were Lactococcus lactis subsp. lactis CECT 539, Lactobacillus casei subsp. casei CECT 4043 (both obtained from the Spanish Type Culture Collection, CECT), and Pedicoccus acidilactici NRRL B-5627 (obtained from the Northern Regional Research Laboratory, NRRL, Peoria, IL, USA). These strains were used as producer of lactic acid (LA) and total antibacterial activity (TAA) in the different culture media.

The three bacterial strains were initially grown at 30 ºC in MRS (de Man, Rogosa, and Sharpe) broth for 24 h, and an aliquot (1 mL) was used to inoculate 50 mL of the corresponding fermentation medium (Table 1). These pre-cultures were incubated at 30 ºC and 200 rpm in an orbital shaker for 12 h. Subsequently, the main fermentations were carried out by inoculating 1 mL of the respective pre-culture into the corresponding fermentation medium (Table 1).

2.2. Culture Media, Fermentation Conditions and Data Collection

The culture media used in the batch cultures of L. lactis CECT 539, Lact. casei CECT 4043, and Ped. acidilactici NRRL B-5627 are summarized in Table 1. During these fermentations, the concentrations of biomass (X(t), in g/L), lactic acid (LA(t), in g/L), and total antibacterial activity (TAA(t), in activity units (AU)/mL), along with culture pH, total sugars (TS(t), in g/L) and nitrogen (TN(t), in g/L), were measured at regular intervals [31,32,33,34,35,36].

The values of these culture variables were used to develop coupled logistic equations from the corresponding uncoupled logistic equations, providing a more accurate description of the time courses of all six variables in the different batch cultures.

The applicability of the coupled logistic equations was further tested using data on culture pH, biomass, and bacteriocin (a component of TAA) produced by Pediococcus acidilactici LB42–923, Lactococcus lactis subsp. lactis ATCC 11454, Leuconostoc carnosum Lm1, and Lactobacillus sakei LB 706 in 18–h batch cultures in TGE medium [23]. Each culture was incubated at its respective optimum temperature. The culture pH, biomass, and bacteriocin data were extracted from the publication by Yang and Ray [23] using GetData Graph Digitizer 2.24 software.

The preparation of the different fermentation media has been described in previous studies by our research team [31,32,33,34,35,36] or in the publication by Yang and Ray [23].

2.3. Analytical Determinations

The concentrations of biomass, X(t); lactic acid, LA(t); total antibacterial activity, AA(t); total sugars, TS(t); and total nitrogen, TN(t), in batch cultures of Lactococcus lactis subsp. lactis CECT 539, Lactobacillus casei CECT 4043, and Pediococcus acidilactici NRRL B-5627 were determined using the methods described in previous studies [31,32,33,34,35,36]. The initial concentrations of TS and TN in TGE medium were measured according to the procedures reported in the same references [31,32,33,34,35,36].

Analytical determinations of biomass and bacteriocins produced by Pediococcus acidilactici LB42–923, Lactococcus lactis subsp. lactis ATCC 11454, Leuconostoc carnosum Lm1, and Lactobacillus sakei LB 706 have been previously described by Yang and Ray [23].

2.4. Models Improvement

2.4.1. Uncoupled Models

Batch biomass growth, X(t), as well as the productions of lactic acid, LA(t), and antibacterial activity, TAA(t), by lactic acid bacteria (LAB) in various culture media have commonly been described using the simplest form of the logistic equation [27,28,29,37], as shown below:

$$X(t)={\displaystyle \frac{K_X}{1+a_X·e^{-b_X·t}}}$$

(1)

$$LA(t)={\displaystyle \frac{K_{LA}}{1+a_{LA}·e^{-b_{LA}·t}}}$$

(2)

$$TAA(t)={\displaystyle \frac{K_{TAA}}{1+a_{TAA}·e^{-b_{TAA}·t}}}$$

(3)

Where K_X, K_LA, and K_TAA represent the maximum biomass concentration (g/L), lactic acid concentration (g/L), and total antibacterial activity (AU/mL), respectively. The constants a_X, a_LA, and a_TAA are dimensionless, whereas b_X, b_LA, and b_TAA denote the rates (h⁻¹) of biomass growth, lactic acid production, and total antibacterial activity synthesis, respectively. The variable t denotes the fermentation time (h).

The decrease in culture pH, pH(t), has usually been described using the von Bertalanffy equation (4) [27].

$$pH\left(t\right)={pH}_f+{\displaystyle \frac{{(pH}_i-{pH}_f)}{e^{b_{pH}·t}}}$$

(4)

Where pH_i and pH_f represent the initial and final culture pH values, respectively; and b_pH denotes the rate constant (h⁻¹) of pH decrease.

In this case, at t = 0 h, pH₀ = pH_i, and as t → ∞, pH_∞ = pH_f.

Because the pH evolution in LAB cultures generally follows a decreasing pattern [23,27,37], similar to that observed for nutrients (total sugars, nitrogen, proteins, and phosphorus) [27,28,37], equation (5), originally proposed for nutrient consumption [28], can also be used to describe pH evolution in LAB cultures:

$$pH\left(t\right)={\displaystyle \frac{{pH}_f}{{1-a_{pH}·e}^{{-b}_{pH}·t}}}$$

(5)

In this case, at t = 0 h:

${pH}_0={\displaystyle \frac{{pH}_f}{1-a_{pH}}}$

and as t → ∞, pH_∞ = pH_f.

In the case of DW medium, equations (4) and (5) adequately fitted the experimental pH data (Figure 1), which exhibited a C-shaped profile. However, these equations did not always provide an optimal description of pH decrease, particularly when the pH followed a typical inverted S-shaped profile, as observed in MRS and TGE broths (Figure 1).

The differences between the pH profiles in DW medium and those in MRS and TGE broths may be attributed to variations in the buffering capacity of these culture media [38], which reflects their resistance to changes in the initial culture pH due to the presence of proteins, amino acids, and salts with buffering properties [39]. Commercial broths contain higher protein concentrations (13.32 ± 0.021 g/L in MRS broth [32] and 10.22 ± 0.61 g/L in TGE broth) than deproteinized diluted whey (2.04 ± 0.08 g/L) [33]. Consequently, the buffering capacity of the commercial media (MRS and TGE) is higher.

Therefore, in this study, model (4) was modified by introducing a constant (BC) in the numerator to account for the buffering capacity of the culture medium on the pH time course. In addition, the denominator of the model was modified to smooth the rate of pH decrease, as follows:

$$pH\left(t\right)={pH}_f+{\displaystyle \frac{{(pH}_i-{pH}_f)·BC}{{1+a_{pH}·e}^{b_{pH}·t}}}$$

(6)

Where a_pH is a dimensionless constant.

With this approach, model (6) describes pH changes in LAB fermentations as a function of organic acid production (mainly lactic acid) and the buffering capacity of the culture medium [40].

At t = 0 h:

${pH}_0={pH}_f+{\displaystyle \frac{{(pH}_i-{pH}_f)·BC}{1+a_{pH}}}$

Considering that at t = 0 h, pH₀ = pH_i, then:

$BC=1+a_{pH}$

Substituting BC into model (6) gives:

$$pH\left(t\right)={pH}_f+{\displaystyle \frac{{(pH}_i-{pH}_f)·(1+a_{pH})}{{1+a_{pH}·e}^{b_{pH}·t}}}$$

Similarly, as t → ∞, pH_∞ = pH_f.

As shown in Figure 1, model (6) provided the best fit to the experimental pH data, accurately describing both C-shaped and inverted S-shaped profiles in the three culture media, with higher R² values. Therefore, this model was used to describe the pH profiles in the different culture media (CW; unsupplemented DW; DW supplemented with glutamic acid or glucose; unsupplemented and glutamic acid-supplemented MPW; and MRS and TGE broths).

To describe the time courses of total sugars, TS(t), and total nitrogen, TN(t), the following models were used:

$$TS(t)={TS}_f+{\displaystyle \frac{{TS}_i-{TS}_f}{e^{b_{TS}·t}}}$$

(7)

$$TN(t)={TN}_f+{\displaystyle \frac{{TN}_i-{TN}_f}{e^{b_{TN}·t}}}$$

(8)

Where TS_i and TS_f, and TN_i, TN_f, represent the initial and final concentrations (g/L) of total sugars (TS) and total nitrogen (TN), respectively. The constants a_TS and a_TN are dimensionless, whereas the parameters b_TS and b_TN represent the rate constants (h⁻¹) for TS and TN consumption, respectively.

2.4.2. Development of Coupled Models

Since biomass growth, the synthesis of both lactic acid and antibacterial activity, culture pH evolution, and the consumption of total sugars and nitrogen are measured simultaneously during fermentation, time-independent relationships among all six variables can be derived by eliminating the fermentation time variable. To achieve this, t is isolated from models (1), (2), (3), (6), (7), and (8).

Starting from model (1):

$$1+a_X·e^{-b_X·t}={\displaystyle \frac{K_X}{X(t)}}$$

$$a_X·e^{-b_X·t}={\displaystyle \frac{K_X}{X(t)}}-1={\displaystyle \frac{K_X-X(t)}{X(t)}}$$

$$e^{-b_X·t}={\displaystyle \frac{K_X-X(t)}{a_X·X(t)}}$$

Taking the natural logarithm yields:

$$-b_X·t=ln\left({\displaystyle \frac{K_X-X\left(t\right)}{a_X·X\left(t\right)}}\right)$$

Therefore:

$$t=-{\displaystyle \frac{1}{b_X}}·ln\left({\displaystyle \frac{K_X-X\left(t\right)}{a_X·X\left(t\right)}}\right)$$

(9)

Applying the same procedure to models (2) and (3) gives:

$$t=-{\displaystyle \frac{1}{b_{LA}}}·ln\left({\displaystyle \frac{K_{LA}-LA\left(t\right)}{a_{LA}·LA\left(t\right)}}\right)$$

(10)

$$t=-{\displaystyle \frac{1}{b_{TAA}}}·ln\left({\displaystyle \frac{K_{TAA}-TAA\left(t\right)}{a_{TAA}·TAA\left(t\right)}}\right)$$

(11)

For models (6) to (8):

$$t={\displaystyle \frac{1}{b_{pH}}}·ln\left({\displaystyle \frac{{pH}_i-pH\left(t\right)+a_{pH}\left({pH}_i-{pH}_f\right)}{a_{pH}·(pH\left(t\right)-{pH}_f)}}\right)$$

(12)

$$t={\displaystyle \frac{1}{b_{TS}}}·ln\left({\displaystyle \frac{{TS}_i-{TS}_f}{TS\left(t\right)-{TS}_f}}\right)$$

(13)

$$t={\displaystyle \frac{1}{b_{TN}}}·ln\left({\displaystyle \frac{{TN}_i-{TN}_f}{TN\left(t\right)-{TN}_f}}\right)$$

(14)

Since LAB growth depends on the evolution of culture pH [21,27,31,32,40], the parameter b_X in model (9) can be corrected to account for this effect. For this purpose, the relationships between b_X and b_pH was derived by equating model (9) with model (12).

Solving for b_X gives:

$$b_X=b_{pH-X}·{\displaystyle \frac{ln\left(\frac{a_X·X\left(t\right)}{K_X-X\left(t\right)}\right)}{ln\left(\frac{{pH}_i-pH\left(t\right)+a_{pH}\left({pH}_i-{pH}_f\right)}{a_{pH}·(pH\left(t\right)-{pH}_f)}\right)}}$$

(15)

Where b_pH–X represents the rate of biomass formation as affected by pH.

Substituting b_X into model (1) gives:

$$X(t)={\displaystyle \frac{K_X}{1+a_X·e^{\left(-b_{pH-X}·\frac{ln\left(\frac{a_X·X\left(t\right)}{K_X-X\left(t\right)}\right)}{ln\left(\frac{{pH}_i-pH\left(t\right)+a_{pH}\left({pH}_i-{pH}_f\right)}{a_{pH}·\left(pH\left(t\right)-{pH}_f\right)}\right)}·t\right)}}}$$

(16)

Proceeding in the same way for lactic acid and antibacterial activity, and considering that lactic acid is produced by LAB cells and that antibacterial activity is a primary metabolite whose production depends on the pH time course [23,27,31,37,41], the synthesis of both compounds can be expressed as:

$$LA(t)={\displaystyle \frac{K_{LA}}{1+a_{LA}·e^{\left({-b}_{X-LA}·\frac{ln\left(\frac{K_{LA}-LA\left(t\right)}{a_{LA}·LA\left(t\right)}\right)}{ln\left(\frac{K_X-X\left(t\right)}{a_X·X\left(t\right)}\right)}·t\right)}}}$$

(17)

$$TAA(t)={\displaystyle \frac{K_{TAA}}{1+a_{AA}·e^{\left(-b_{X-TAA}·\frac{ln\left(\frac{K_{TAA}-TAA\left(t\right)}{a_{TAA}·TAA\left(t\right)}\right)}{ln\left(\frac{K_X-X\left(t\right)}{a_X·X\left(t\right)}\right)}·t+b_{pH-TAA}·\frac{ln\left(\frac{K_{TAA}-TAA\left(t\right)}{a_{TAA}·TAA\left(t\right)}\right)}{ln\left(\frac{{pH}_i-pH\left(t\right)+a_{pH}\left({pH}_i-{pH}_f\right)}{a_{pH}·\left(pH\left(t\right)-{pH}_f\right)}\right)}·t\right)}}}$$

(18)

Where b_X–LA represents the rate of lactic acid production as affected by biomass growth, and b_X–TAA and b_pH–TAA represent the rates of total antibacterial activity synthesis as affected by biomass growth and pH, respectively.

Since the reduction in culture pH is caused by lactic acid production by growing LAB cells, the relationship between b_pH and b_X can be obtained by equating models (9) and (12). Therefore, the culture pH evolution could be expressed as:

$$pH(t)={pH}_f+{\displaystyle \frac{{(pH}_i-{pH}_f)·(1+a_{pH})}{1+a_{pH}·e^{\left(-b_{LA-pH}·\frac{ln\left(\frac{{pH}_i-pH\left(t\right)+a_{pH}\left({pH}_i-{pH}_f\right)}{a_{pH}·\left(pH\left(t\right)-{pH}_f\right)}\right)}{ln\left(\frac{K_{LA}-LA\left(t\right)}{a_{LA}·LA\left(t\right)}\right)}·t\right)}}}$$

(19)

Where b_LA–pH represents the rate of pH decline as influenced by lactic acid production.

Considering that total sugars and nitrogen are consumed by the growing LAB cells, the coupled logistic models for nutrient consumption can be expressed as:

$$TS(t)={TS}_f+{\displaystyle \frac{{TS}_i-{TS}_f}{e^{\left(-b_{X-TS}·\frac{ln\left(\frac{{TS}_i-{TS}_f}{TS\left(t\right)-{TS}_f}\right)}{ln\left(\frac{K_X-X\left(t\right)}{a_X·X\left(t\right)}\right)}·t\right)}}}$$

(20)

$$TN(t)={TN}_f+{\displaystyle \frac{{TN}_i-{TN}_f}{e^{\left(-b_{X-TN}·\frac{ln\left(\frac{{TN}_i-{TN}_f}{TN\left(t\right)-{TN}_f}\right)}{ln\left(\frac{K_X-X\left(t\right)}{a_X·X\left(t\right)}\right)}·t\right)}}}$$

(21)

Where b_X–TS and b_X–TN represent the rates of total sugars and total nitrogen consumption, respectively, as influenced by biomass growth.

Although more complex relationships among the six culture variables could be derived, the simple coupled models (16)–(21) provide a tighter relationship among the variables, indicating that the change in one variable is driven by the effect of another.

2.5. Model Parameters Determination and Model Evaluation

Experimental data on the kinetics of biomass growth, lactic acid production, total antibacterial activity synthesis, culture pH, and the consumption of total sugars and nitrogen from triplicate fermentations in different culture media [31,32,33,34,35,36] were used to develop the coupled models and estimate the corresponding model parameters. The values of the constants were initially obtained by minimizing the sum of squared differences between the experimental and model-predicted values using the nonlinear least-squares (quasi-Newton) method implemented via the Solver tool in Microsoft Excel 2016 (Microsoft Corp., Redmond, WA, USA). To avoid overestimation of the initial values of each culture variable (X, LA, TAA, pH, TS, and TN) during the fitting of the models to the experimental data, a constraint was applied in all cases requiring that the model-predicted values at time zero (Yp₀) match the corresponding experimental values (Yexp₀).

The corresponding R² values for the models were determined using the same Microsoft Excel spreadsheet.

The uncoupled models were independently fitted to the corresponding experimental data by minimizing the sum of squared differences (SSD) between the experimental values (Yexpᵢ) of biomass (X), lactic acid (LA), total antibacterial activity (TAA), culture pH, total sugars (TS), and nitrogen (TN), and the corresponding predicted values (Ypᵢ) generated by the models:

$$SSD=\sum _{i=1}^n{\left({Yexp}_i-{Yp}_i\right)}^2$$

In the case of coupled models, which share some parameters (models (16)–(21)), the model constants were simultaneously estimated by minimizing the total sum of squared differences (SSD_T) between the experimental values (expᵢ) of the six culture variables (X, LA, TAA, pH, TS, and TN) and their corresponding predicted values (pᵢ) generated by the models:

$${SSD}_T={\sum }_{i=1}^n{\left({Xexp}_i-{Xp}_i\right)}^2+{\sum }_{i=1}^n{\left({LAexp}_i-{LAp}_i\right)}^2+{\sum }_{i=1}^n{\left({TAAexp}_i-{TAAp}_i\right)}^2+{\sum }_{i=1}^n{\left({pHexp}_i-{pHp}_i\right)}^2+{\sum }_{i=1}^n{\left({TSexp}_i-{TSp}_i\right)}^2+{\sum }_{i=1}^n{\left({TNexp}_i-{TNp}_i\right)}^2$$

Lower SSD values and higher R² values indicate a better fit and greater accuracy of the models in describing the trends observed in the experimental data.

3. Results and Discussion

The efficacy of the coupled logistic models (16)–(21) in describing the time courses of six culture variables—namely biomass growth (X), product formation (LA and TAA), culture pH, and nutrient consumption (TS and TN)—was demonstrated using experimental data from three lactic acid bacteria strains (CECT 539, CECT 4043, and NRRL B–5627) cultivated in different culture media [31,32,33,34,35,36].

To further validate the accuracy of the proposed models in describing the six culture variables, the results obtained from the coupled models (16)–(21) were compared with those generated by the corresponding unmodified logistic models (1), (2), (3), (6), (7), and (8).

In addition, the models were validated using data from batch cultures of Pediococcus acidilactici LB42–923 (a pediocin AcH producer), Lactococcus lactis subsp. lactis ATCC 11454 (a nisin producer), Leuconostoc carnosum Lm1 (a leuconocin Lcm1 producer), and Lactobacillus sakei LB 706 (a sakacin A producer), grown in TGE broth at uncontrolled pH [23].

As only experimental data on growth, bacteriocin production, and pH time courses were available in this case [23], only models (16), (18), and (19) were validated, and the corresponding results were compared with those obtained using the unmodified models (1), (3), and (6), respectively.

Considering that lactic acid is primarily a growth-associated product, as it is predominantly synthesized by LAB strains during the active growth phase [30,42], model (19) was modified by substituting the lactic acid time-course data, LA(t), and its corresponding constants (K_LA and a_LA) with those of biomass growth (X(t), K_X, and a_X), as indicated below:

$$pH(t)={pH}_f+{\displaystyle \frac{{(pH}_i-{pH}_f)·(1+a_{pH})}{1+a_{pH}·e^{\left(-b_{X-pH}·t·\frac{ln\left(\frac{{pH}_i-pH\left(t\right)+a_{pH}\left({pH}_i-{pH}_f\right)}{a_{pH}·(pH\left(t\right)-{pH}_f)}\right)}{ln\left(\frac{K_X-X\left(t\right)}{a_X·X\left(t\right)}\right)}\right)}}}$$

(22)

Where b_X–pH represents the rate of pH decline as influenced by biomass growth.

3.1. Modeling Batch Growth, Product Formation, and Nutrient Consumption by Lactococcus lactis CECT 539 in Concentrated Whey and MRS Broth

3.1.1. Uncoupled Models

Figure 2A shows the time courses of experimental biomass growth, lactic acid production, total antibacterial activity synthesis, culture pH, and consumption of total sugars and nitrogen by strain CECT 539 in concentrated whey (CW) and MRS broth, together with the predictions of models (1), (2), (3), (6), (7), and (8).

For the CW culture (Figure 2A and Table 2), models (1), (3), and (6) provided an adequate description of the experimental data for biomass growth (R²_X = 0.9846, SSD–X = 2.29 × 10^–3), total antibacterial activity production (R²_TAA = 0.9838, SSD–TAA = 1.18), and culture pH (R²_pH = 0.9952, SSD–pH = 6.12 × 10^–3). However, models (2), (7), and (8) did not satisfactorily describe the time courses of lactic acid (R²_LA = 0.9500, SSD–LA = 2.58 × 10^–2), total sugars (R²_TS = 0.9638, SSD–TS = 3.34), and total nitrogen (R²_TN = 0.9512, SSD–TN = 8.27 × 10^–3), respectively (Figure 2A and Table 2).

The failure of model (2) to adequately describe the kinetics of lactic acid synthesis in CW medium can be attributed to the fact that lactic acid production did not follow the typical logistic-type profile observed in MRS broth (see the upper right panel of Figure 2A). Even when the constraint requiring that the model-predicted values at time zero (Yp₀) match the corresponding experimental values (Yexp₀) was removed, the R²_LA value increased to 0.9793 and the SSD–LA decreased to 1.18 × 10⁻². However, the predicted initial lactic acid concentration increased to 0.082 g/L, which is considerably higher than the experimental value (0.016 g/L).

The relatively low R² obtained for models (7) and (8) can be explained by the fact that the time courses of total sugars (TS) and total nitrogen (TN) did not strictly follow a C-shaped profile during CW fermentation (see the lower left and lower right panels of Figure 2A).

With regard to the culture in MRS broth, the modeling results were slightly better than those obtained in CW medium (Figure 2A), with higher R² values (0.9979, 0.9981 and 0.9809) obtained for models (1), (2), and (7), respectively (Table 2). However, the predictions of models (3), (6), and (8) for the kinetics of total antibacterial activity (TAA) synthesis, culture pH, and total nitrogen (TN) were similar in both culture media, yielding comparable R² values.

The differences observed in the predictions of the models (1), (2), (3), (6), (7), and (8) between fermentations conducted in the two culture media can be attributed solely to the degree of agreement between the experimental time-course data for the corresponding culture variables and the mathematical structure of each model.

3.1.2. Coupled Models

Uncoupled models assume that each variable evolves independently with its own logistic time constant. In LAB fermentations, the growth, the synthesis of lactic acid and total antibacterial activity, pH drop and consumption of nutrients (TS and TN) are not independent processes evolving only with time, because, fermentation is a closed biochemical network, not six independent time–driven phenomena [27,29,30,31,32]. Then, developing of coupled models could reflect this biological interdependence.

In this study, the following simple biologically linked through cause–effect relationships between the six variables were assumed: i) pH affects cell growth kinetics, ii) lactic acid production decreases pH, iii) primary metabolite production (LA) depends on growth, iv) antibacterial activity depends on both biomass and pH evolution and, v) biomass formation drives sugar and nitrogen consumption.

The coupled models (16)–(21) were fitted simultaneously to the experimental data for the six culture variables (X, LA, TAA, pH, TS, and TN), incorporating biologically meaningful cause–effect relationships among them. The statistical performance of these models in concentrated whey (CW) and MRS broth is summarized in Table 3, and the corresponding fits are presented in Figure 2B.

In CW medium, the coupled formulation significantly improved the overall description of the fermentation process compared with the uncoupled approach (Table 2 and Table 3). The coupled biomass model provided an excellent fit (R²_X = 0.9962; SSD–X = 6.25 × 10⁻⁴), outperforming the uncoupled model (R²_X = 0.9846; SSD–X = 2.29 × 10⁻³). The approximately one-order-of-magnitude reduction in SSD–X indicates that incorporating the pH-dependent correction of the growth parameter (b_pH–X = 0.29 h⁻¹) enhanced the representation of growth dynamics. This confirms that biomass evolution in CW is strongly influenced by culture acidification.

In contrast, the coupled model did not substantially improve the description of lactic acid production (R²_LA = 0.9492; SSD–LA = 3.99 × 10⁻²) compared with the uncoupled model (R²_LA = 0.9500; SSD–LA = 2.58 × 10⁻²). The lower goodness-of-fit reflects the deviation of lactic acid production in CW from a strictly logistic profile (Figure 2A and 2B). Although coupling with biomass (b_X–LA = 0.51 h⁻¹) ensures biological consistency, it does not fully capture the irregular kinetic pattern observed experimentally. This suggests that additional metabolic constraints or substrate limitations may influence acid synthesis in CW.

A marked improvement was observed for TAA production (R²_TAA = 0.9943; SSD–TAA = 0.47) compared with the uncoupled model (R²_TAA = 0.9838; SSD–TAA = 1.18). The simultaneous dependence on biomass (b_X–TAA = 0.23 h⁻¹) and pH (b_pH–TAA = 0.12 h⁻¹) enhanced the predictive capacity of the model. The lower SSD indicates that antibacterial activity is better described when both growth and acidification effects are considered, supporting its classification as a growth-associated metabolite influenced by environmental pH.

Although the fit remained acceptable (R²pH = 0.9750), the coupled pH model showed a lower R² and higher SSD–pH (4.13 × 10⁻²) than the uncoupled formulation (R²pH = 0.9952; SSD–pH = 6.12 × 10⁻³). This slight deterioration is likely due to the imposed dependence of pH on lactic acid production (b_LA–pH = 0.83 h⁻¹), which did not follow a strictly logistic time profile (Figure 2A and 2B). Nevertheless, the coupled model ensures mechanistic coherence between acid production and pH decline.

A substantial improvement was observed for sugar consumption. The R²_TS increased from 0.9638 to 0.9926, while SSD–TS decreased from 3.34 to 0.56. By directly linking total sugar (TS) consumption to biomass formation (b_X–TS = 0.40 h⁻¹), the coupled model more accurately represented substrate depletion, confirming that sugar utilization is primarily growth-driven rather than purely time-dependent.

An even more pronounced improvement was observed for nitrogen consumption. The R²_TN increased from 0.9512 to 0.9954, and SSD–TN decreased by nearly one order of magnitude (from 8.27 × 10⁻³ to 8.88 × 10⁻⁴). This indicates that total nitrogen (TN) depletion is strongly correlated with biomass production (b_X–TN = 0.44 h⁻¹), and that the coupling strategy effectively captures this relationship.

Overall, in CW medium, the coupled models improved the description of biomass, TAA, TS, and TN, maintained a similar fit for lactic acid, and slightly reduced the statistical performance for pH. Nevertheless, the coupled system provides a more coherent and biologically consistent representation of the fermentation process (Table 2 and Table 3).

In MRS broth, where the uncoupled models already provided excellent fits, the improvements introduced by coupling were less pronounced, although still relevant.

For biomass, R²_X decreased slightly from 0.9979 to 0.9938, and the SSD–X increased modestly. Nevertheless, the fit remained excellent, and the growth–pH interaction (b_pH–X = 0.51 h⁻¹) confirms that acidification influences growth even in a nutritionally rich medium.

For lactic acid (LA) synthesis, the coupled model slightly reduced R²_LA (0.9915 vs. 0.9981), although the fit remained very good. Because LA production in MRS follows a well-defined logistic profile, the uncoupled model was already sufficient (Figure 2A); thus, coupling mainly improved mechanistic interpretation rather than statistical performance.

In the case of the TAA model, a slight improvement was observed (R²_TAA = 0.9958 vs. 0.9856), accompanied by a marked reduction in SSD–TAA (0.12 vs. 42.07). The dual dependence on biomass and pH significantly enhanced model accuracy, highlighting the synergistic regulation of antibacterial compound synthesis (Figure 2A and 2B).

The coupled pH model slightly improved the fit (R²_pH = 0.9969 vs. 0.9949; lower SSD–pH), indicating that linking pH decline to LA production (b_LA–pH = 0.68 h⁻¹) is appropriate in MRS broth (Table 2 and Table 3).

For total sugars (TS) and total nitrogen (TN), the coupled models showed clear improvements. For TS, R²_TS increased from 0.9809 to 0.9966, while SSD–TS decreased from 1.28 to 0.18. For TN, R²_TN increased from 0.9563 to 0.9957; with a substantial reduction in SSD–TN. These results further confirm that nutrient depletion is tightly coupled to biomass formation and that the coupled models provide a more accurate description of substrate utilization kinetics.

Overall, across both media, the coupled models improved or maintained high R² values for most variables and significantly reduced SSD for biomass, TAA, TS, and TN. However, they slightly reduced statistical flexibility for LA (in CW) and pH (in CW) due to the imposed biological constraint (initial experimental LA value = predicted LA value).

The most notable improvements were observed for nutrient consumption and antibacterial activity, variables that strongly depend on biomass and environmental conditions.

To further demonstrate the effectiveness of the coupled models, kinetic data from fermentations conducted in media containing different initial concentrations of total nitrogen (TN) or total sugars (TS) were used for model fitting. Under these conditions, the time-course profiles of TS and TN differed substantially [31,32,33,34,35,36], allowing a more rigorous evaluation of the ability of the coupled models to describe the kinetics of the six culture variables.

Accordingly, both the uncoupled and coupled models were applied to describe the evolution of the six culture variables in DW and MPW media supplemented with different concentrations of glutamic acid, as well as in DW medium supplemented with varying concentrations of glucose.

3.2. Modeling of Batch Growth, Product Formation, Culture pH Drop, and Nutrient Consumption by Strain CECT 539 in DW and MPW Media Supplemented with Different Concentrations of Glutamic Acid

3.2.1. Uncoupled Models

The results of fitting the uncoupled models (1), (2), (3), (6), (7), and (8) to the fermentation kinetics of L. lactis CECT 539 in both unsupplemented and glutamic acid–supplemented DW and MPW media are presented in Figure 3A and 3B and Table 4 and Table 5, respectively.

In both media, biomass growth and the production of lactic acid and total antibacterial activity followed well-defined sigmoidal (S-shaped) logistic profiles, whereas culture pH exhibited a C-shaped profile (Figure 3A and 3B). The good agreement between the experimental time courses and the model predictions for these four variables (X, LA, TAA, and pH) supported the fitting of models (1), (2), (3), and (6), respectively (Table 4 and Table 5).

However, supplementation of DW and MPW media with glutamic acid resulted in more pronounced inverted S-shaped profiles for total sugars (TS) and total nitrogen (TN) compared with the corresponding unsupplemented media (Figure 3A and 3B). This behavior affected the fitting of the uncoupled TS and TN models in DW cultures, which showed the poorest fits, with R² values ranging from 0.9066 to 0.9545 for TS and from 0.9041 to 0.9479 for TN (Table 4).

In MPW cultures, the R² values obtained for the TS model (0.9614–0.9765) and the TN model (0.9645–0.9846) (Table 5) were higher than those observed in DW media (Table 4). This difference was likely due to the wider variation ranges of TS and TN in DW, DW–GA1, DW–GA2, DW–GA3, DW–GA4, and DW–GA5 media (Figure 3A) compared with the corresponding MPW media (Figure 3B).

In DW, DW–GA1, DW–GA2, DW–GA3, DW–GA4, and DW–GA5 media, TS concentrations ranged from 26.10 to 23.56, 26.16 to 22.70, 26.19 to 22.10, 25.87 to 21.64, 25.44 to 20.60, and 25.76 to 20.55 g/L, respectively. The corresponding TN concentrations ranged from 0.473 to 0.323, 0.570 to 0.334, 0.664 to 0.332, 0.758 to 0.335, 0.867 to 0.346, and 0.945 to 0.332 g/L (Figure 3A). In contrast, in MPW media, TS concentrations ranged from 5.32 to 3.05, 5.32 to 2.84, 5.32 to 2.41, 5.36 to 2.26, 5.31 to 2.09, and 5.31 to 1.68 g/L, respectively, while TN concentrations ranged from 0.668 to 0.492, 0.765 to 0.492, 0.857 to 0.490, 0.951 to 0.493, 1.043 to 0.481, and 1.137 to 0.468 g/L (Figure 3B).

3.2.2. Coupled Models

The results of fitting the coupled models (16)–(21) to the fermentation kinetics of L. lactis CECT 539 in both unsupplemented and glutamic acid–supplemented DW and MPW media are presented in Figure 4A and 4B and Table 6 and Table 7, respectively.

A comparison between the uncoupled (Table 4 and Table 5) and coupled models (Table 6 and Table 7) clearly shows that coupling the models substantially improved the overall goodness of fit, particularly for substrate (TS) and nitrogen (TN) consumption, while maintaining similar predictive performance for biomass (X), lactic acid (LA), total antibacterial activity (TAA), and pH.

For biomass in DW media, R²_X values (0.9945–0.9985) were above 0.9900 for all cultures with the uncoupled model (Table 4). When the coupled model was applied (Table 6), slight decreases in R² were observed specifically in DW–GA2 (from 0.9985 to 0.9955), DW–GA3 (from 0.9966 to 0.9899), and DW–GA4 (from 0.9966 to 0.9872). In DW and DW–GA5, the differences were minimal (DW increased from 0.9945 to 0.9984, whereas DW–GA5 changed slightly from 0.9945 to 0.9941) (Table 4 and Table 6).

In MPW media, R² values for X were generally comparable between approaches (Table 5 and Table 7). A clear improvement with the coupled model was observed in MPW–GA2 (from 0.9797 to 0.9968). In the remaining MPW cultures, differences were small: slight decreases occurred in MPW (0.9900 to 0.9920, actually a small improvement), MPW–GA1 (0.9925 to 0.9919), MPW–GA3 (0.9861 to 0.9913, improvement), MPW–GA4 (0.9892 to 0.9886), and MPW–GA5 (0.9938 to 0.9940, slight improvement). Overall, the changes were minor and confirm that biomass was already well described by the uncoupled logistic model.

For LA and TAA, the coupled models produced similar or slightly higher R² values compared with the uncoupled models, especially in DW cultures (Table 4, Table 5, Table 6 and Table 7). The improvement was more noticeable for TAA in DW medium, where R²_TAA increased to values up to 0.9991 (Table 4 and Table 6), indicating that linking TAA production to biomass growth and pH enhanced the model’s descriptive capacity.

The most remarkable improvement was observed for TS and TN. In DW cultures, the coupled models increased R² dramatically to values between 0.9982–0.9993 for TS and 0.9973–1.0000 for TN (Table 4 and Table 6). A similar trend was observed in MPW media, where R² values for TS increased from 0.9614–0.9765 (uncoupled) to 0.9938–0.9981 (coupled), and for TN from 0.9645–0.9846 to approximately 0.9942–0.9968 (Table 5 and Table 7). These results indicate that substrate and nitrogen consumption were more accurately described when directly linked to biomass formation (Figure 4A and 4B).

The SSD values confirm the trends observed for R². For X, LA, TAA, and pH, SSD values were of the same order of magnitude in both approaches, with slight improvements in some cases under the coupled formulation (Table 4 and Table 6).

However, for TS and TN in DW medium, the reduction in SSD was substantial. For example, SSD–TS decreased from values between 1.38–4.12 (uncoupled) to 1.86×10⁻²–7.63×10⁻² (coupled), representing a reduction of approximately two orders of magnitude. An even more pronounced decrease was observed for TN, where SSD values dropped from 10⁻³–10⁻² to 10⁻⁴–10⁻⁵ (Table 4 and Table 6).

Similar, though less dramatic, reductions were observed in MPW cultures (Table 5 and Table 7). This marked decrease demonstrates the superior predictive accuracy of the coupled models for nutrient consumption.

The estimated constants for biomass (K_X and a_X) followed similar increasing trends with glutamic acid concentration in both modeling strategies, indicating that coupling did not distort the biological interpretation of growth stimulation by glutamic acid. Likewise, K_LA and K_TAA increased consistently with supplementation in both approaches, confirming the positive effect of glutamic acid on metabolite production.

In the coupled models, the relatively stable values of the coupling coefficients (b_X–LA, b_X–TAA, b_LA–pH, b_X–TS), which link product formation, pH decrease, and nutrient consumption to biomass growth and/or lactic acid production, across different media and supplementation levels suggest consistent physiological relationships among growth, acidification, and substrate utilization. In particular, the strong coupling between biomass and TS/TN consumption (b_X–TS and b_X–TN) explains the substantial improvement in the fit obtained for these variables.

These results emphasize that the coupled modeling approach provides a more robust and integrated framework for describing the kinetics of growth, product formation, and nutrient consumption in batch cultures of L. lactis CECT 539, particularly when interactions among variables play a significant role.

3.3. Comparison between Uncoupled and Coupled Models in Fermentations of Glucose-Supplemented DW Media with Strain CECT 539

In the previous set of cultures, the fermentation media were supplemented with different concentrations of glutamic acid to generate distinct TN time-course profiles (Figure 3 and Figure 4) and to evaluate the ability of the coupled models to describe these variations.

In the present analysis, the focus was on modeling the time courses of the six culture variables in glucose-supplemented DW media (Figure 5A and 5B), where different TS consumption profiles were observed [36]. To assess the performance of the coupled models under these conditions, a comparison between the fitting results obtained with the uncoupled (Table 8) and coupled models (Table 9) was carried out.

The comparison reveals clear differences in model performance, particularly for substrate (TS) and nitrogen (TN) consumption (Figure 5A and 5B).

For biomass, both approaches provided very good fits (R² generally > 0.9700) (Table 8 and Table 9). However, in contrast to the glutamic acid experiments, the uncoupled model performed slightly better in some cultures. For example, in DW, R²_X decreased from 0.9984 (uncoupled) to 0.9824 (coupled), and SSD–X increased from 8.10×10⁻⁴ to 1.27×10⁻². A similar trend was observed in DW–G20.

In the case of lactic acid production, both approaches yielded high R² values (>0.9800 in most cases). The coupled model improved the fit in several cultures, notably in DW (R²_LA increased from 0.9926 to 0.9970, with SSD–LA reduced by half) and in DW–G20 (R²_LA increased from 0.9821 to 0.9972, with SSD–LA decreasing from 4.19×10⁻² to 5.27×10⁻³). In other cases (e.g., DW–G5), the uncoupled model showed slightly lower SSD–LA values. Overall, coupling LA production to biomass formation (b_X–LA) improved consistency to some extent but did not substantially alter the predictive capacity, as the uncoupled model already provided an adequate fit (Table 8 and Table 9).

For TAA, both models performed well (R² ≈ 0.9900). The coupled model improved R²_TAA values in most glucose-supplemented cultures (e.g., DW–G15: 0.9931 to 0.9964; DW–G20: 0.9868 to 0.9983), although SSD–TAA values were not always lower (e.g., DW–G5). The inclusion of coupling terms (b_X–TAA and b_pH–TAA) enabled a more mechanistic representation of TAA formation, particularly at higher glucose concentrations, where growth and acidification effects are more pronounced (Figure 5).

For pH, both approaches showed good agreement with the experimental data (R²_pH between 0.9600 and 0.9900). Differences between models were relatively small. In some cases (DW–G10 and DW–G15), the coupled model slightly improved R²_pH and reduced SSD–pH, whereas in others (DW and DW–G20) the uncoupled model performed marginally better. Thus, coupling pH decrease to lactic acid production (b_LA–pH) provided a biologically consistent formulation but did not systematically enhance the statistical indicators.

The most significant differences were observed for TS. In the uncoupled model, R²_TS values were highly variable and, in some cases, low, particularly in DW (R²_TS = 0.8194, SSD–TS = 36.16). Even in glucose-supplemented media, R²_TS values ranged from 0.9425 to 0.9758. In contrast, the coupled model dramatically improved the fit for all cultures, with R²_TS values between 0.9972 and 0.9995 and a substantial reduction in SSD–TS (e.g., from 36.16 to 9.46×10⁻² in DW). This improvement confirms that TS consumption is strongly growth-associated and is more accurately described when directly linked to biomass formation (b_X–TS).

A similar pattern was observed for TN. The uncoupled model yielded moderate R²_TN values (0.8766–0.9209), indicating limited descriptive capacity (Table 8). The coupled model increased R²_TN to 0.9955–0.9989 and reduced SSD–TN by approximately one to two orders of magnitude (e.g., from 6.21×10⁻³ to 1.87×10⁻⁴ in DW). This marked improvement demonstrates that nitrogen consumption is closely associated with biomass growth and cannot be adequately represented as an independent process.

Overall, in glucose-supplemented DW cultures, the uncoupled models were generally sufficient to describe biomass growth, lactic acid production, TAA formation, and pH evolution. However, they were less accurate in capturing the dynamics of substrate (TS) and nitrogen (TN) consumption, particularly under conditions where glucose supplementation altered TS profiles.

The coupled models provided a more integrated and physiologically consistent description of the fermentation process. Although improvements for X, LA, TAA, and pH were moderate and not systematic, the enhancement in the fitting of TS and TN was substantial, as reflected by the marked increase in R² values and the pronounced reduction in SSD.

Therefore, as observed in the glutamic acid experiments, the coupled modeling approach is particularly advantageous for describing nutrient consumption kinetics when strong interactions exist among growth, substrate utilization, and metabolite production.

3.4. Performance of Uncoupled and Coupled Models in Describing Fermentations in CW, DW, and MRS Media with Lactobacillus casei CECT 4043 and Pediococcus acidilactici NRRL B–5627

Experimental data from batch cultures of Lact. casei CECT 4043 [34,35] and Ped. acidilactici NRRL B–5627 [31,32] grown in concentrated whey (CW), diluted whey (DW), and MRS broth were used to further evaluate the ability of the coupled models to describe the fermentation kinetics of additional LAB strains.

The time courses for strain CECT 4043 are presented in Figure 6A and 6B, while the estimated parameters and the corresponding R² and SSD values are summarized in Table 10 (uncoupled models) and Table 11 (coupled models).

For strain CECT 4043, the coupled models consistently yielded higher R² values and lower SSD values than the uncoupled models for the six culture variables (X, LA, TAA, pH, TS, and TN). Although the uncoupled biomass model already provided high R²_X values (0.9676–0.9909), the coupled model further improved the fit in CW and MRS media (R²_X up to 0.9961), with a marked reduction in SSD–X (e.g., from 1.05 × 10⁻³ to 2.20 × 10⁻⁴ in CW).

For lactic acid production, both approaches performed well. However, the coupled models slightly increased R²_LA in CW and DW and substantially reduced SSD–LA in DW. In MRS medium, although R²_LA remained high (>0.9900), SSD–LA increased in the coupled model, suggesting that under this specific condition the additional coupling did not significantly improve the dispersion of residuals despite maintaining a strong overall correlation.

In the case of total antibacterial activity (TAA), improvement with the coupled formulation was particularly evident in CW and MRS media, where R²_TAA increased to ≥0.9972 and SSD–TAA decreased considerably (e.g., from 1.58 to 0.44 in MRS).

Regarding culture pH, both models described pH evolution accurately (R²_pH ≈ 0.9900). Differences in SSD–pH were small, indicating that pH dynamics were already well captured by the uncoupled formulation, although slight improvements were observed in CW and MRS media with the coupled approach.

The culture variables TS and TN showed the most pronounced improvement when modeled using the coupled models (Figure 6A and 6B; Table 10 and Table 11). In MRS medium, the R²_TS value increased dramatically (from 0.8823 to 0.9973), and SSD–TS decreased from 29.48 to 1.32. Similarly, TN fitting improved across all media, with R²_TN values approaching or exceeding 0.9900 and substantial reductions in SSD–TN.

A similar trend was observed for Ped. acidilactici NRRL B–5627 in the three culture media (Figure 7A and 7B; Table 10 and Table 11). For biomass production, R²_X values were already very high for both modeling approaches (>0.9800); however, the coupled models generally reduced SSD–X, indicating improved precision even when correlation was strong.

For lactic acid production, the coupled approach significantly improved the fit in CW and MRS media (R²_LA up to 0.9944 and 0.9936, respectively), accompanied by marked reductions in SSD–LA (e.g., from 1.25 × 10⁻² to 9.51 × 10⁻⁴ in CW).

Substantial improvements were also observed for TAA synthesis in CW and MRS media, with R²_TAA values increasing to ≥0.9960 and SSD–TAA decreasing considerably (e.g., from 4417.129 to 1182.13 in MRS). Although SSD–TAA values remained relatively high due to the magnitude of TAA measurements, the relative reduction reflects improved predictive accuracy.

Both pH models provided excellent fits (R² ≥ 0.9600), and improvements with coupling were moderate, reflecting the already robust performance of the uncoupled models (Table 10 and Table 11).

Modeling of TS showed one of the most striking improvements in MRS medium, where R²_TS increased from 0.6789 to 0.9867 and SSD–TS decreased from 23.89 to 0.80. Similar improvements were observed in CW and DW (Table 10 and Table 11). In addition, the coupled models substantially improved TN prediction, with R²_TN values reaching up to 0.9997 in CW and SSD–TN decreasing by several orders of magnitude.

Overall, while the uncoupled models already described biomass growth and pH evolution reasonably well, coupling the kinetic expressions to biologically meaningful interactions (e.g., X–LA, X–TS, X–TN, and LA–pH relationships) significantly enhanced the mechanistic consistency and predictive robustness of the models.

Therefore, the coupled modeling strategy provides a more reliable and biologically realistic framework for describing LAB fermentations in CW, DW, and MRS media.

3.5. Performance of Uncoupled and Coupled Models in Describing Fermentations in TGE Broth with Pediococcus acidilactici LB42–923, Lactococcus lactis subsp. lactis ATCC 11454, Leuconostoc carnosum Lm1, and Lactobacillus sakei LB 706

The performance of the uncoupled (models 1, 3, and 6) and coupled (models 16, 18, and 22) formulations was evaluated using batch fermentation data in TGE broth for strains LB42–923, ATCC 11454, Lm1, and LB 706 (Table 12 and Table 13; Figure 8). In this case, the comparison was based on the coefficients of determination (R²) and the sums of squared differences (SSD) for biomass (X), bacteriocin production (B or TAA), and pH.

For all four strains, the uncoupled models already provided very good fits, with R²_X values ranging from 0.9786 to 0.9939 (Table 12). However, the coupled models systematically improved biomass prediction. Specifically, R²_X increased from 0.9920 to 0.9964 for strain LB42–923, from 0.9786 to 0.9911 for strain ATCC 11454, from 0.9845 to 0.9975 for strain Lm1, and from 0.9939 to 0.9996 for strain LB 706. Correspondingly, SSD–X values decreased from 0.16 to 8.36 × 10⁻² (LB42–923), from 3.25 × 10⁻² to 1.25 × 10⁻² (ATCC 11454), from 6.56 × 10⁻² to 1.23 × 10⁻² (Lm1), and markedly from 3.28 × 10⁻² to 2.77 × 10⁻³ (LB 706).

Greater differences between the modeling approaches were observed for bacteriocin production. For strain LB42–923, R²_TAA increased from 0.9774 to 0.9894, while SSD–TAA decreased dramatically from 34.35 to 0.18. For strain ATCC 11454, R²_TAA increased from 0.9886 to 0.9999, with SSD–TAA decreasing from 4.89 to 4.59 × 10⁻², indicating an almost perfect fit with the coupled model. In the case of strain Lm1, R²_TAA improved from 0.9902 to 0.9974, and SSD–TAA decreased from 13.21 to 3.55. However, for strain LB 706, R²_TAA decreased slightly (from 0.9948 to 0.9821) and SSD–TAA increased (from 4.86 to 12.62) with the coupled formulation. This suggests that, for this strain, bacteriocin production may not be strongly influenced by the additional coupling terms, or that the uncoupled model already described the kinetics adequately.

The pH profiles were described with high accuracy by the uncoupled models, with R²_pH values ranging from 0.9945 to 0.9985 (Table 12). The coupled models produced only minor changes: slight increases in R²_pH for strains LB42–923 and ATCC 11454, and small decreases for strains Lm1 and LB 706 (table 13). Similarly, SSD–pH values exhibited only modest variations, without a consistent trend toward improvement.

These results indicate that pH dynamics were adequately captured by both modeling approaches and that explicitly coupling pH to biomass did not substantially enhance predictive performance in most cases.

5. Conclusions

The results of this study demonstrate that coupled modeling approaches effectively describe biomass growth, metabolite production (LA and TAA), pH evolution, and nutrient (TS and TN) consumption by LAB across different fermentation media. Although more mathematically restrictive than uncoupled models, coupled models provide a more integrated and physiologically consistent representation of LAB fermentation dynamics by linking biomass growth, pH variation, nutrient utilization, and metabolite formation. These stronger inter-variable relationships reduce arbitrary parameter independence and enhance physiological interpretability, particularly in complex media such as concentrated whey.

By incorporating biologically meaningful relationships among variables, the coupled approach provides mechanistic insights beyond simple statistical fitting and enhances predictive capability under varying media compositions and supplementation strategies. Furthermore, coupled models enable the representation of both the decreasing C-shaped and inverted S-shaped profiles commonly observed in nutrient consumption across most LAB cultures.

For practical applications such as process optimization and scale-up, coupled models offer a more robust, integrated, and reliable framework, especially when nutrient management is critical. By emphasizing key physiological relationships—biomass–substrate coupling, acidification–product interactions, and nutrient utilization—they support improved medium formulation, supplementation strategies, and strain selection.