Evaluation and Mathematical Analysis of a Four-Dimensional Lotka–Volterra-like Equation Designed to Describe the Batch Nisin Production System

Nisin, an antibacterial compound produced by Lactococcus lactis strains, has been approved by the US Food and Drug Administration to be used as a safe food additive to control the growth of undesirable pathogenic bacteria. Nisin is commonly described as a pH-dependent primary metabolite since its production depends on growth and culture pH evolution. However, the relationships between bacteriocin synthesis (BT), biomass production (X), culture pH, and the consumption of the limiting nutrient (total nitrogen: TN) have not been described until now. Therefore, this study aims to develop a competitive four-dimensional Lotka–Volterra-like Equation (predator-prey system) to describe these complex relationships in three series of batch fermentations with L. lactis CECT 539 in diluted whey (DW)-based media. The developed four-dimensional predator-prey system accurately described each individual culture, providing a good description of the relationships between pH, TN, X, and BT, higher values for R2 and F-ratios, lower values (<10%) for the mean relative percentage deviation modulus, with bias and accuracy factor values approximately equal to one. The mathematical analysis of the developed equation showed the existence of one asymptotically stable equilibrium point, and the phase’s diagram obtained did not show the closed elliptic trajectories observed in biological predator-prey systems.


Introduction
Nisin, a bacteriocin produced by Lactococcus lactis strains, has a wide antibacterial activity against food spoilage and pathogenic bacteria. For this reason, this biomolecule has been recognized by the US Food and Drug Administration as a natural and safe biopreservative in food products, being allowed in the USA and several European Union countries. The advantages of using nisin in foods include the reduction in both the thermal treatment and addition of chemicals to food products and an increase in their shelf life [1].
For high nisin production at low cost, it is necessary to know the relationship between the main culture variables, which could be elucidated with the use of appropriate mathematical models. This could also allow proper monitoring and control of these bioprocesses [2]. Different mathematical models have been commonly used to describe the kinetics of growth (e.g., Verhulst, Gompertz, Richards, Bertalanffy, Weibull, and Monod) and (iii) Biomass (X(t)) could be considered as one predator that grows logistically competing with nisin for the nitrogen source and depending on the culture pH and the TN source concentration [5], as follows: (iv) Nisin (BT(t)) could be considered as the second predator that is produced by the biomass but competes with it for the nitrogen source and depends on the culture pH and the TN source concentration [5], as follows:

Modeling the Batch Nisin Production System in Different Series of Fermentations in DW Media Using a Global Set of Model Parameters
The capability of the four-dimensional predator-prey system (1)- (4) to describe the batch nisin production system was first assessed by adjusting a unique global set of equation parameters (general equation) to the entire set of experimental data in different series of batch cultures. The fermentations were performed in culture media (Tables 1 and 2) prepared with diluted whey (DW) containing the following different initials: (i) concentrations of glucose (DW-G series) [16], (ii) concentrations of total sugars and phosphorous (DW-TS-TP series) [17], or (iii) concentrations of MRS broth nutrients (DW-MRS series) [18]. Table 1. Initial concentrations (mean ± standard deviations) of total sugars (TS), nitrogen (TN), phosphorous (TP), and proteins (Pr) in culture media prepared with deproteinized diluted whey (DW) and concentrated mussel-processing wastes (CMPW). The data (symbols) corresponding to the four series of batch cultures [16][17][18] and the corresponding predictions (dashed lines) of the developed four-dimensional predator-prey system (1)-(4) are shown in Figures 1-4. The values for the constants and the statistical analysis of each equation in each series of cultures are shown in Table 3. Table 3. Parameter values (as estimates ± confidence intervals) calculated with the global set of model parameters of the four-dimensional predator-prey system (1)- (4) to describe the batch nisin production system in the different series of fermentations.  The parameter value is considered statistically significant if its corresponding p-value is lower than 0.05.  Table 3. Solid lines were obtained by adjusting the four-dimensional predator-prey system to the experimental data corresponding to each individual culture (see parameter values in Table 4). Reproduced with permission from Costas et al. [16], Appl. Microbiol. Biotechnol.; published by Springer Nature, 2016.  Table 3. Solid lines were obtained by adjusting the four-dimensional predator-prey system to the experimental data corresponding to each individual culture (see parameter values in Table 4). Reproduced with permission from Costas et al. [16], Appl. Microbiol. Biotechnol.; published by Springer Nature, 2016. The parameter value is considered statistically significant if its corresponding p-value is lower than 0.05.

Figure 2.
Experimental data (symbols) of culture pH, TN consumption, and X and Nis production by L. lactis CECT 539 in the first seven batch cultures of the experimental matrix (Table 2), corresponding to the fermentation series DW-TS-TP. Dashed lines drawn through the experimental data are predictions of the global four-dimensional predator-prey systems (1)-(4) obtained with the parameters shown in Table 3. Solid lines were obtained by adjusting the four-dimensional predator-prey system (1)-(4) to the experimental data corresponding to each individual culture (see parameter values in Table 5). Reproduced with permission from Costas et al. [17], 3-Biotech; published by Springer Nature, 2018. Table 5. Statistically significant (p < 0.05) parameter values (as estimates ± confidence intervals) calculated with the four-dimensional predator-prey system (1)-(4) for each individual culture corresponding to the four factorial points and the four axial points of the experimental matrix (Table 2)   The parameter value is considered statistically significant if its corresponding p-value is lower than 0.05.
Mathematics 2022, 10, x FOR PEER REVIEW 11 of 35 Figure 3. Experimental data (symbols) of culture pH, TN consumption, and X and Nis production by L. lactis CECT 539 in the last six batch cultures of the experimental matrix (Table 2) and in the optimum conditions (OC), corresponding to the fermentation series DW-TS-TP. Dashed lines drawn through the experimental data are predictions of the global four-dimensional predator-prey system (1)-(4) obtained with the parameters shown in Table 3. Solid lines were obtained by adjusting the four-dimensional predator-prey system (1)-(4) to the experimental data corresponding to each individual culture (see parameter values in Tables 5 and 6). Reproduced with permission from Costas et al. [17], 3-Biotech; published by Springer Nature, 2018.  (Table 2) and in the optimum conditions (OC), corresponding to the fermentation series DW-TS-TP. Dashed lines drawn through the experimental data are predictions of the global four-dimensional predator-prey system (1)-(4) obtained with the parameters shown in Table 3. Solid lines were obtained by adjusting the fourdimensional predator-prey system (1)-(4) to the experimental data corresponding to each individual culture (see parameter values in Tables 5 and 6). Reproduced with permission from Costas et al. [17], 3-Biotech; published by Springer Nature, 2018. Table 6. Statistically significant (p < 0.05) parameter values (as estimates ± confidence intervals) calculated with the four-dimensional predator-prey system (1)-(4) for each individual culture corresponding to the five center points (TS = 37.0 g/L, TP = 0.43 g/L) of the experimental matrix (Table 2) and to the optimum conditions (TS = 22.6 g/L, TP = 0.46 g/L) of the series of fermentation DW-TS-TP.    The parameter value is considered statistically significant if its corresponding p-value is lower than 0.05.  Table 3. Solid lines were obtained by adjusting the four-dimensional predator-prey system (1)-(4) to the experimental data corresponding to each individual culture (see parameter values in Table 7).  Table 3. Solid lines were obtained by adjusting the four-dimensional predator-prey system (1)-(4) to the experimental data corresponding to each individual culture (see parameter values in Table 7).  The parameter value is considered statistically significant if its corresponding p-value is lower than 0.05.

Center Points
When the global Equations (1)-(4) were set to describe the time course of the culture pH, TN, X, and BT in the DW-G, DW-TS-TP, and DW-MRS series of cultures, the results obtained were not satisfactory. Thus, although in some cases, statistically significant values (p < 0.0001) for both the parameters and global pH equations were obtained, the values of R pH 2 and F-ratio were relatively low, the RPDM values were almost always higher than 10, and both the B f and A f values were generally far from one (Table 3). In addition, the pH, TN, X, and BT trajectories predicted by the global Equations (1)-(4) for the three series of cultures showed a clear deviation from the experimental pH, TN, X, and BT data (dashed lines in Figures 1-4).
These observations suggest that the four-dimensional predator-prey system could not be used as a general equation to describe the nisin production system in the DW-G, DW-TS-TP, and DW-MRS series of batch fermentations.
These unsatisfactory results could be related to the different initial compositions of the media used in each series of cultures: DW media supplemented with different initial concentrations of glucose (Figure 1), TS and TP (Figures 2 and 3), and MRS broth nutrients ( Figure 4). Therefore, it could be considered that, in each series of cultures, each fermentation was independent of the other ones since the fermentation substrates used were different. So that the different initial media composition in the following three series of batch fermentations (DW-G, DW-TS-TP, and DW-MRS) modulated the time-course of the culture variables: pH and the concentrations of total nitrogen, biomass, and nisin (Figures 1-4). In this way, it is well known that the initial culture conditions affect the evolution of these culture variables (pH, TN, X, and BT) in different ways [8,[16][17][18]21]. For example, the pH drop depends on the presence and interaction between some compounds (salts, organic acids, proteins, and free amino acids) with buffering capacity in the culture medium [22] and organic acid production by growing cells [16][17][18]. TN consumption during fermentation depends on the initial medium composition, mainly the type and concentration of the nitrogen source [5,[16][17][18]21] and culture pH [18,20]. In fact, the consumption of TN [18] or amino acids [20] in L. lactis strains was maximal when the culture pH reached values between 5.8 and 6.5, and decreased abruptly for high and low pH values.
On the other hand, biomass production depends on different factors, including the initial medium composition (concentration and type of nutrients, mainly carbon, nitrogen, and phosphorous sources), initial and final pH values in the cultures, pH evolution, and production of inhibitory compounds [5,[16][17][18]. Nisin synthesis depends not only on the time course of biomass concentration, but also on (i) the amount of biomass produced, (ii) the initial concentration and type of nutrient (carbon, nitrogen, and phosphorous sources), and (iii) initial and final pH value, pH evolution, and pH drop generated in the cultures [5,8,12,[16][17][18]. So that the specific effects of these factors on the response variables (culture pH, TN consumption, biomass, and nisin production) could be non-synchronous, producing a different change in the time course of the latter variables and, consequently, in their relationships.
For example, the buffering capacity (BC), which is a measure of the resistance of the culture medium to pH changes, affects biomass and nisin synthesis differently. On the one hand, the increase in BC favors biomass production since the cultures remain longer within the optimum pH range (between 5.8 and 6.5) for nutrient consumption for L. lactis CECT 539 [18,20]. On the other hand, these high pH values inhibit bacteriocin synthesis, which was higher at an optimum pH value of 4.90 in DW medium. The latter was due to the need for a low pH value to favor the maturation of the nisin molecule [5,18]. The value of this optimum final pH for nisin production depends on the producer strain and composition of the culture medium [8,12,23].
In addition, it has been observed that higher pH drops (rpH) enhance nisin production [12,18] before the cultures reached an inappropriate pH for survival and cell growth of L. lactis [24]. Thus, in the series of fermentations DW-G, DW-TS-TP, and DW-MRS ( Figures 1-4), it can be observed that the evolution of culture pH, biomass production, and nisin synthesis in each culture was different.
For these reasons, it is very difficult to develop a general four-dimensional predatorprey system to explain the variations in the time courses of the four variables (culture pH, TN, X, and BT) for each or all series of cultures. In addition, with the use of a general four-dimensional equation, the effect of different initial culture conditions on the evolution of the four dependent variables could not be explained, leading to a misinterpretation of the kinetics of the cultures.
To solve this problem, we first fitted the four-dimensional predator-prey system (1)-(4) to each individual culture of each series of fermentation to accurately determine how the values of the different parameters change with changes in the initial culture conditions (concentrations of glucose, TP and TS, and MRS broth nutrients). Afterward, we intend to modify the four-dimensional predator-prey system (1)-(4) (when this was possible) by including a term for the specific effect of the initial culture conditions on the evolution of pH, TN, X, and BT.

Modeling the Batch Nisin Production System in Individual Cultures Corresponding to Each Series of Fermentations
When the four-dimensional predator-prey system (1)-(4) was used to describe the relationships between the four response variables (pH, TN, X, and BT) in each individual culture, both equations and the values of the parameters were statistically significant (p < 0.050), with R 2 and F-values considerably higher, and B f and A f values~1 (Tables 4-7).
In addition, the predictions of the four-dimensional predator-prey system (1)-(4) for each response variable (solid lines in Figures 1-4) were in perfect agreement with the corresponding experimental data. This indicates that the developed four-dimensional predator-prey system (1)-(4) is consistent and robust enough to accurately describe the trend observed in the experimental data of culture pH, TN, X, and BT.
The results obtained for each series of fermentations are discussed below. Table 4 shows the parameter values as well as the statistical analysis obtained when the four-dimensional predator-prey system (1)-(4) was fitted to the experimental data of cultures pH, TN, X, and BT in the DW-G cultures.

Series of Fermentation DW-G
In this case, all values of the parameters in Equations (1)-(4) were significant (p < 0.05) and considerably higher values for R pH 2 (between 0.9968 and 0.9980), R TN 2 (between 0.9987 and 0.9998), R X 2 (between 0.9990 and 1.0000), and R BT 2 (between 0.9992 and 1.0000) were obtained. In addition, the values of B f and A f calculated for Equations (1)-(4) were~1 and the RPDM values were considerably lower than 10% (Table 4). Therefore, it could be considered that the use of the four-dimensional predator-prey system (1)-(4) accurately described the trend observed for the culture pH, TN, X, and BT in the DW-G cultures.
Regarding the equation parameters, it can be noted that the values of a and b in Equation (1) had a negative sign and their absolute values decreased (Table 4). This is because the culture pH drop (the difference between the initial and final pH values) decreased and the final pH increased with the increase in the initial concentration of glucose [G 0 ] (from 0 to 25 g/L) in the media (Figure 1). In addition, the value of c decreased from 0.039 to 0.002 with the increase in [G 0 ] due to the inhibition that increasing glucose concentration produced on the growth of L. lactis CECT 539 and, consequently, on lactic acid production [16], causing a gradual reduction in the pH drop in the culture media ( Figure 1).
In the case of Equation (2), negative values for d, e, f, and g, and positive values for h were obtained ( Table 4). The decrease in the absolute values obtained for d and e could be related to the reduction in the TN consumption rates and the increase in the final TN values observed with the increase in [G 0 ] (Figure 1). Similarly, the values of f, g, and h were almost similar for all cultures because the TN consumption decreased with the increase in [G 0 ], in parallel with the reduction in biomass production, nisin synthesis, and pH drop (Figure 1).
Equation (3) also provides an accurate description of biomass production in each culture. In this case, the values of i decreased because of the inhibition that the increasing initial glucose concentrations produced on the growth rate of L. lactis (Figure 1). On the other hand, the value of j did not vary, indicating that the reduction in the growth rate was proportional to the reduction in the maximum biomass concentration produced in the different cultures (Figure 1).
Similarly, a constant value for k was obtained, indicating that the competition between biomass production and nisin synthesis for the nitrogen source was very similar in the different glucose-supplemented cultures. As observed in Figure 1, the reduction in growth caused by the increase in [G 0 ] was proportional to that observed in nisin synthesis because this bacteriocin was produced in this series of cultures as a pH-dependent primary metabolite [16].
As expected, the values of l and m decreased (Table 4), due to the decrease in the pH drops and TN consumption caused by the reduction in biomass production ( Figure 1).
The detailed analysis of the results obtained for Equation (4) showed a decrease in the values of n in agreement with the reduction in nisin production with the increase in [G 0 ] (Figure 1). The constant value obtained for the o constant could be explained by a proportional decrease in the nisin production rate and maximum nisin levels produced by L. lactis CECT 539. The constant value obtained for p is in perfect agreement with the constant value obtained for k in Equation (3), indicating again that the competition between biomass production and nisin synthesis for the nitrogen source was very similar in the different glucose-supplemented cultures.
The decrease in the values obtained for q could be explained by the fact that the reduction in pH drops with the increase in [G 0 ] negatively affected the synthesis of nisin. In the unsupplemented culture ([G 0 ] = 0), the final pH value was 4.73, which is in perfect agreement with the optimum final pH between 4.78 and 4.90 observed in L. lactis cultures in whey [8]. In contrast, a constant value was obtained for r (Table 4), indicating proportional TN consumption for nisin synthesis.
These results corroborate the affirmation that nisin was produced by L. lactis CECT 539 as a pH-dependent primary metabolite [16] and indicate that the TN consumption was proportional to the production of bacteriocin and biomass.
As discussed above, each culture variable (pH, TN, X, and BT) evolved, describing a similar profile in the different fermentations of the DW-G series ( Figure 1). However, the rates of culture pH (rpH(t)) and TN (rTN(t)) decreased, and biomass (rX(t)) and nisin (rBT(t)) production in the glucose-supplemented cultures did not exhibit a clear trend compared with the respective rates in the culture in the unsupplemented culture ( Figure 5).
From the detailed observation of Figure 5, it can be noted that the highest rpH(t), rX(t), rBT(t)), and rTN(t) were obtained in the unsupplemented DW substrate during the first 4, 8, 9, and 12 h of fermentation, respectively. However, after these times, the four rates in the unsupplemented DW medium were lower than the corresponding rates calculated in the glucose-supplemented cultures.
For this reason, the four-dimensional Equations (1)-(4) could not be modified by including a term for explaining the inhibitory effect of the increase in [G 0 ] on the rates rpH(t), rTN(t), rX(t), and rBT(t)). compared with the respective rates in the culture in the unsupplemented culture ( Figure  5). From the detailed observation of Figure 5, it can be noted that the highest rpH(t), rX(t), rBT(t)), and rTN(t) were obtained in the unsupplemented DW substrate during the first 4, 8, 9, and 12 h of fermentation, respectively. However, after these times, the four rates in the unsupplemented DW medium were lower than the corresponding rates calculated in the glucose-supplemented cultures.
For this reason, the four-dimensional Equations (1)-(4) could not be modified by including a term for explaining the inhibitory effect of the increase in [G0] on the rates rpH(t), rTN(t), rX(t), and rBT(t)). Tables 5 and 6 show the results obtained when the four-dimensional Equations (1)-(4) was fitted to each individual culture of the series of fermentation DW-TS-TP. The predictions of Equations (1)-(4) are shown as solid lines in Figures 2 and 3. As observed before for the series of fermentation DW-G, the DW-TS-TP cultures were satisfactorily described using this modeling procedure (Tables 5 and 6).

Series of fermentation DW-TS-TP
The values of a, b, and c in Equation (1), as expected, depended on the initial chemical composition of the media (mainly the initial concentrations of TS and TP). So that the highest a value was obtained in the culture performed at the optimum conditions (TS = 22.6 g/L, TP = 0.46 g/L, Table 6), in which the highest pH drop (difference between the initial and final pH value) was generated (Figures 2 and 3). The calculated values for b and c varied between 0.007 and 0.013 and between 0.033 and 0.057, respectively (Tables 5  and 6), which were dependent on the growth of L. lactis in the different cultures ( Figures  2 and 3).
The highest values for d and f in Equation (2) were obtained at the optimum conditions (Table 6), in which the highest amounts of TN and biomass were consumed (0.230 g/L) and produced (0.716 g/L), respectively (Figures 2 and 3). This suggests that the total  Tables 5 and 6 show the results obtained when the four-dimensional Equations (1)-(4) was fitted to each individual culture of the series of fermentation DW-TS-TP. The predictions of Equations (1)-(4) are shown as solid lines in Figures 2 and 3. As observed before for the series of fermentation DW-G, the DW-TS-TP cultures were satisfactorily described using this modeling procedure (Tables 5 and 6).

Series of Fermentation DW-TS-TP
The values of a, b, and c in Equation (1), as expected, depended on the initial chemical composition of the media (mainly the initial concentrations of TS and TP). So that the highest a value was obtained in the culture performed at the optimum conditions (TS = 22.6 g/L, TP = 0.46 g/L, Table 6), in which the highest pH drop (difference between the initial and final pH value) was generated (Figures 2 and 3). The calculated values for b and c varied between 0.007 and 0.013 and between 0.033 and 0.057, respectively (Tables 5 and 6), which were dependent on the growth of L. lactis in the different cultures (Figures 2 and 3).
The highest values for d and f in Equation (2) were obtained at the optimum conditions (Table 6), in which the highest amounts of TN and biomass were consumed (0.230 g/L) and produced (0.716 g/L), respectively (Figures 2 and 3). This suggests that the total nitrogen source consumption depended on the initial composition of the fermentation medium, as indicated before [2,5,[16][17][18]. The value of e depended on the initial media composition, but the values of g and h were constant, indicating a proportional consumption of TN for nisin production and a similar effect of pH on TN assimilation.
The maximum growth rate was observed in the culture performed at TS = 22.60 g/L, TP = 0.46 g/L, and in accordance with this, the highest value for i in Equation (3) was obtained (Tables 5 and 6). In addition, the values of j and l varied as a function of the initial TS and TP concentrations in the different culture media, but the value of the competition coefficient k was constant in the different cultures, indicating a proportional efficiency of TN utilization for biomass production and nisin synthesis. The constant value calculated for m suggests that the nitrogen source consumption was directly correlated with biomass production.
When Equation (4) was fitted to the experimental data of nisin synthesis, the highest n value was obtained in the culture performed at the optimum conditions (Tables 5 and 6) due to the highest bacteriocin production rate observed in this culture (Figures 2 and 3). The coefficient o did not show a significant variation since the relationship between the nisin synthesis rates and the maximum bacteriocin levels produced in the cultures was almost constant. Additionally, the coefficients p, q, and r exhibited constant values, suggesting constancy in the competition between biomass and nisin production for the TN source and in the effect of pH and TN consumption on bacteriocin synthesis.
However, in the case of nisin production, the values of RPDM corresponding to some cultures were higher than 10% (Tables 5 and 6), due to the lack of fit between the experimental and calculated values observed during the first 7 h of fermentation (Figures 2 and 3). This was probably because nisin production was quantified by a photometric bioassay using an indicator strain [16] and, consequently, the experimental error in determining nisin titers could be greater than that of the analytical methods used in pH, total nitrogen, and biomass measurements. So that in nisin determination, the differences between the experimental and predicted values during the first 7 h of fermentation were low (Figures 2 and 3), but the experimental nisin data in this interval, used as the denominator in Equation (8), were also considerably low, increasing the RPDM value (Tables 5 and 6).
In the series of cultures DW-TS-TP, the effects of the initial TS and TP concentrations on both the growth and bacteriocin production were described by empirical quadratic equations [17]. The inclusion of terms for explaining these effects in the four-dimensional predator-prey system (1)-(4) could contribute to obtaining a general equation for describing the evolution of the four response variables (pH, TN, X, and BT). However, this approach has several drawbacks since too large equations could be obtained, and information about the true relationship between the dependent variables (pH(t), TN(t), X(t), BT(t)), and the own essence of the LV equation would be lost.
In fact, the rates of culture pH drop (rpH), total nitrogen consumption (rTN), and biomass (rX) and nisin (rBT) production in the different experiments (1-14) did not show a clear dependence on changes in initial TS and TP concentrations ( Figure 6).

Series of Fermentation DW-MRS
In this series of fermentation, the absolute values of the constant a increased when the DW medium was supplemented with MRS nutrients from 0 to 50% (fermentations

Series of Fermentation DW-MRS
In this series of fermentation, the absolute values of the constant a increased when the DW medium was supplemented with MRS nutrients from 0 to 50% (fermentations DW-MR0%, DW-MRS25%, and DW-MRS50%) since the pH drop increased slightly in these cultures from 2.27 to 2.36 (Table 7, Figure 4). The absolute values of b decreased slightly; meanwhile, the constant c decreased from 0.103 to 0.022 because of the increase in the buffering capacity of media supplemented with increasing MRS nutrient concentrations from 0 to 50% [18]. This counteracted the reduction in pH values due to acid organic production by the nisin-producing strain.
However, the values of a and b had a positive value in the following cultures (fermentations DW-MR75%, DW-MRS100%, and DW-MRS125%), and c increased in DW-MRS100% and DW-MRS125% fermentations. The sign change observed for a and b from negative in the first three fermentations (DW-MR0%, DW-MRS25%, and DW-MRS50%) to positive in the latter three cultures (DW-MR75%, DW-MRS100%, and DW-MRS125%) could be related to the change in the trajectories described by the culture pH. These pH trajectories evolved from convex curves in the first three cultures to inverted S-curves in the three latter cultures ( Figure 4).
Thus, the values of a, b, and c increased in the fermentation DW-MRS100% compared to the fermentation DW-MRS75% due to the increase in the growth of L. lactis. However, in the fermentation DW-MRS125%, the increase in the buffering capacity of the supplemented media counteracted the effect of lactic acid production by the growing strain [18]. For this reason, the values of a, b, and c decreased slightly (Table 7).
Regarding Equation (2), it can be noted that the values of d and e (with negative signs) increased and decreased, respectively, from the fermentation DW-MRS0% to DW−75% (Table 7). However, the sign of both constants becomes positive in the following fermentations due to the change in the curvature of the TN trajectories from convex curves (first four cultures) to inverted S-curves (fermentations DW-MRS100% and DW-125%), as observed before for the culture pH curves (Figure 4).
In fermentations DW-MRS0% to DW-75%, the values of f and h increased with nutrient supplementation since the addition of MRS nutrients led to an increase in the TN consumption, growth, and pH gradient. The latter probably affected the TN consumption rate since the consumption of this nutrient depends on the culture pH, as explained above [18,20]. The constant e varied as a function of the value of d and the lowest TN concentration reached in the cultures, and g was approximately constant (−0.009 ± 0.002) in the four fermentations (Table 7). In comparison with fermentation DW-MRS100%, fermentation DW-125% provided higher values of d and e, in agreement with the increase in TN consumption in the latter culture. However, the value of f decreased while those of g and h were constant in both cultures.
Equation (3) provided increasing values for i due to the increase in the growth rates caused by the increase in MRS nutrient supplementation; meanwhile, the value of j depended on the values of i and the maximum biomass level reached in each culture. The values of k (0.007 ± 0.000) and l (0.087 ± 0.002) were almost constant, and m showed an increasing trend.
In Equation (4), the values of n increased with the increase in nutrient supplementation, indicating a stimulation in nisin production, and the values of the constant o varied depending on the values of n and the maximum nisin titers reached in the cultures. In addition, the values of p and r did not vary, but the constant q increased (Table 7), probably because of the changes in the trajectories described by the culture pH that affected the evolution of nisin production in the cultures (Figure 4), as commented above [5,8,12].
In this series of cultures, the increase in the MRS nutrients added into the DW medium affected both the evolution of the culture variables and the final concentrations of biomass and nisin obtained (Figure 4), as well as the rates rpH, rTN, rX, and rBT (Figure 7). The rates of culture pH decrease exhibited a transition from exponential decay-shaped curves (in the DW25, DW50, and DW75 media) to bell-shaped curves (in the DW100 and DW125 media); meanwhile, the rTN, rX, and rBT profiles showed bell-shaped curves [25]. However, the rpH, rTN, rX, and rBT profiles did not show an appreciable relationship (linear, quadratic, sigmoidal, etc.) with the initial MRS nutrient concentration ([Nut] 0 ) in the medium.
because of the changes in the trajectories described by the culture pH that affected the evolution of nisin production in the cultures (Figure 4), as commented above [5,8,12].
In this series of cultures, the increase in the MRS nutrients added into the DW medium affected both the evolution of the culture variables and the final concentrations of biomass and nisin obtained (Figure 4), as well as the rates rpH, rTN, rX, and rBT ( Figure  7). The rates of culture pH decrease exhibited a transition from exponential decay-shaped curves (in the DW25, DW50, and DW75 media) to bell-shaped curves (in the DW100 and DW125 media); meanwhile, the rTN, rX, and rBT profiles showed bell-shaped curves [25]. However, the rpH, rTN, rX, and rBT profiles did not show an appreciable relationship (linear, quadratic, sigmoidal, etc.) with the initial MRS nutrient concentration ([Nut]0) in the medium.  Therefore, in this case, it is also difficult to develop a general four-dimensional predator-prey system describing the evolution of pH, TN, X, and BT as a function of the initial concentrations of MRS nutrients.

Mathematical Analysis of the Four-Dimensional Lotka-Volterra Equation
After demonstrating the feasibility of the designed four-dimensional Lotka-Volterra equation to describe the batch nisin production system in different batch cultures, the following step was focused on the mathematical study of the equation by determining its approximate solutions and analyzing its equilibrium points and trajectories around the stable equilibrium points.

Generalized Four-Dimensional Lotka-Volterra Equation
Given two column vectors of R n , x, and y, their component-by-component product can be defined by the following: . .
x n y n      Given x ∈ R n and A ∈ M n×n (R) we define the following: where A i is the i-th row of the matrix A.
The generalized Lotka-Volterra equations are given by the system of differential equations given by the following: where b is a column vector and A is a square matrix.

Lemma 1.
Suppose A is nonsingular. The Jacobian of F(x) = x·(b − Ax)satisfies the following:

Proof. Note that
which proves the Lemma 1.
As a consequence of the above lemma and the linearization theorem of Liapunov and Poincaré, the following theorem follows: Theorem 1. If A is nonsingular, the generalized Lotka-Volterra equations have a unique equilibrium point with all their entries different from zero, given by x = A −1 b. If, in addition, the matrix A· x has all its eigenvalues with a positive real part, then x is asymptotically stable [26]. Note that, since J F ( x) = −A· x that corresponds to the situation in which all the eigenvalues of the Jacobian have a negative real part.
As observed in Figure 9B, the four culture variables (pH, TN, X, and BT) gradually reach a plateau phase and stabilize at the above-mentioned stable equilibrium point. In fact, the coordinates of the latter point are in perfect agreement with the experimental values of the culture variables (pH = 4.730, TN = 0.319 g/L, X = 0.482 g/L, BT = 22.905 BU/mL) obtained after 17 h of incubation, from which the culture variables stabilized (Figure 1).
The phase graphs ( Figure 10) show that there is an interaction between the four culture variables over time. Thus, the relationships between TN vs. pH and BT vs. X were directly proportional because TN concentration decreased with the decrease in pH, and BT increased with the growth of L. lactis since this bacteriocin was produced in the logarithmic phase of growth as a primary metabolite (Figure 1).  As observed in Figure 9B, the four culture variables (pH, TN, X, and BT) graduall reach a plateau phase and stabilize at the above-mentioned stable equilibrium point. I fact, the coordinates of the latter point are in perfect agreement with the experimenta values of the culture variables (pH = 4.730, TN = 0.319 g/L, X = 0.482 g/L, BT = 22.90 BU/mL) obtained after 17 h of incubation, from which the culture variables stabilized (Fig  ure 1).
The phase graphs ( Figure 10) show that there is an interaction between the four cu  directly proportional because TN concentration decreased with the decrease in pH, and BT increased with the growth of L. lactis since this bacteriocin was produced in the logarithmic phase of growth as a primary metabolite (Figure 1). In contrast, the relationships X vs. pH, BT vs. pH, X vs. TN, and BT vs. TN were inversely proportional because the biomass and nisin synthesis increased with the decrease in pH and increase in TN consumption, which reduced the concentration of the nitrogen source (Figure 1).
On the other hand, closed elliptic orbits previously observed in classic predator-prey models were not observed when the four-dimensional LV-like equation was used to describe the batch nisin system, because these microbial relationships are not periodic, as those of the populations of predators and prey in biological systems [27]. In contrast, the relationships X vs. pH, BT vs. pH, X vs. TN, and BT vs. TN were inversely proportional because the biomass and nisin synthesis increased with the decrease in pH and increase in TN consumption, which reduced the concentration of the nitrogen source ( Figure 1).
On the other hand, closed elliptic orbits previously observed in classic predatorprey models were not observed when the four-dimensional LV-like equation was used to describe the batch nisin system, because these microbial relationships are not periodic, as those of the populations of predators and prey in biological systems [27].

Study of the Parameter Values for Which There Is an Asymptotically Stable Solution
In what follows, we will establish, by means of a Monte Carlo study [28,29], which conditions must be satisfied by the parameters defining the matrix A and the vector b so that the Lotka-Volterra equation has a single asymptotically stable equilibrium point with all its positive (and non-zero) coordinates.
If we denote by {p i } 18 i=1 the generic parameters (i.e., a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r), and by {p i0 } 18 i=1 the current parameters, the aim of our study is to look for positive values, R 1 and R 2 , as large as possible, so that if . . ,18 Then x = A −1 b satisfies that all its entries are positive and the matrix A· x has all its eigenvalues with a positive real part.
For this purpose, we fixed R 1 , and the step h = 0.001, and we determined the first value of a natural k such that if R 2 =k·h, and 10.000 random values were taken for each one of the intervals has any of its entries nonpositive or the matrix A· x has some eigenvalue with a nonpositive real part.
We observed that if R 1 > 0.12, there are no asymptotically stable equilibrium points with all their non-zero coordinates. The left and right parts of Figure 11 show the curve R 1 → R 2 and the curves R 1 → 1 − R 1 and R 1 → 1 + R 2 , respectively. As can be observed in these figures, the range of the parameters' validity is approximately equal to 0.12|p i0 |, i = 1,2, . . . ,18. In what follows, we will establish, by means of a Monte Carlo study [28,29], which conditions must be satisfied by the parameters defining the matrix A and the vector b so that the Lotka-Volterra equation has a single asymptotically stable equilibrium point with all its positive (and non-zero) coordinates.
If we denote by the generic parameters (i. e., a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p,q,r), and by the current parameters, the aim of our study is to look for positive values, R1 and R2, as large as possible, so that if Then = satisfies that all its entries are positive and the matrix • has all its eigenvalues with a positive real part.
For this purpose, we fixed R1, and the step h = 0.001, and we determined the first value of a natural k such that if R2 =k·h, and 10.000 random values were taken for each one of the intervals (1 − ), < 0), then = has any of its entries nonpositive or the matrix • has some eigenvalue with a nonpositive real part.
We observed that if R1 > 0.12, there are no asymptotically stable equilibrium points with all their non-zero coordinates. The left and right parts of Figure 11 show the curve R1⟶ R2 and the curves R1⟶1−R1 and R1⟶1+R2, respectively. As can be observed in these figures, the range of the parameters' validity is approximately equal to 0.12| |, i = 1,2,…,18.

Microorganisms, Culture Media, and Inoculum Preparation
In this work, Lactococcus lactis CECT 539 and Carnobacterium piscicola CECT 4020 were used as the nisin-producing strain and target bacterium (in the nisin activity bioassay), respectively. Both strains were obtained from the Spanish Type Culture Collection (CECT) and cultured at 30 °C in MRS (de Man, Rogosa and Sharpe) agar slants or broth.
Diluted whey (DW) and concentrated mussel-processing waste (CMPW) were used to prepare the different culture media (Table 1). Sterilization (121 °C/15 min) of these substrates led to the precipitation of a protein fraction that interfered with biomass measurements. For this reason, the precipitated material was removed by acidification of the DW and CMPW substrates to pH 4.5 with 5 N HCl, heating (121 °C/15 min) and centrifugation (12,000 × g for 15 min) [16][17][18]. Figure 11. Plots of the curves R 1 → R 2 (left part), and R 1 → 1 − R 1 and R 1 → 1+ R 2 (right part).

Microorganisms, Culture Media, and Inoculum Preparation
In this work, Lactococcus lactis CECT 539 and Carnobacterium piscicola CECT 4020 were used as the nisin-producing strain and target bacterium (in the nisin activity bioassay), respectively. Both strains were obtained from the Spanish Type Culture Collection (CECT) and cultured at 30 • C in MRS (de Man, Rogosa and Sharpe) agar slants or broth.
Diluted whey (DW) and concentrated mussel-processing waste (CMPW) were used to prepare the different culture media (Table 1). Sterilization (121 • C/15 min) of these substrates led to the precipitation of a protein fraction that interfered with biomass measurements. For this reason, the precipitated material was removed by acidification of the DW and CMPW substrates to pH 4.5 with 5 N HCl, heating (121 • C/15 min) and centrifugation (12,000 × g for 15 min) [16][17][18].
Given that the nisin-producing strain is not an amylolytic bacterium, the glycogen contained in the CMPW was enzymatically hydrolyzed to produce a glucose-containing substrate as described in Costas et al. [17].
To prepare the different fermentation substrates, DW medium was supplemented with the following nutrients: (i) different amounts of glucose to obtain 5, 10, 15, 20, and 25 g glucose/L of medium (series of fermentation DW-G) [16], (ii) different volumes of CMPW medium (101.33 g glucose/L) and amounts of KH 2 PO 4 to obtain initial total sugars and phosphorous concentrations between 22.61 and 51.35 g/L, and 0.24 and 0.63 g/L, respectively (series of fermentation DW-TS-TP) [17], and (iii) MRS broth nutrients (except glucose and Tween 80) at 25, 50, 75, 100, and 125% (w/v) of their standard concentrations in the complex substrate to produce the DW25, DW50, DW75, DW100, and DW125 media (series of fermentation DW-MRS) [18]. In these three series of fermentation, control cultures in unsupplemented DW medium were performed to obtain data for the comparisons [16][17][18]. Tables 1 and 2 show the mean compositions of the resulting culture media used in this work.
To prepare the preculture, cells of L. lactis CECT 539 from MRS agar slants were used to inoculate sterile MRS broth (10 mL) and incubated at 30 • C for 12 h with shaking at 200 rpm. After that, 50 mL of inoculum medium (which was, in each case, similar to the corresponding fermentation substrate) were inoculated with 1 mL of the preculture and subsequently incubated for 12 h at 30 • C with shaking at 200 rpm. An appropriate volume of the latter culture was used to inoculate the corresponding fermentation medium used in the different batch cultures to give an initial viable cell count of 1.5 × 10 9 colony-forming units/mL [16][17][18].

Batch Cultures
The different experimental data used in this research were collected from previous batch cultures of L. lactis CECT 539 [16][17][18]. The fermentations were conducted in duplicate in 250 mL Erlenmeyer flasks that contained 50 mL of the corresponding DW-based medium. After inoculation, the flasks were incubated at 30 •

Analytical Methods
The corresponding analytical methods used to measure culture pH and concentrations of total nitrogen, biomass, and nisin were previously described in Costas et al. [16].

Statistical Significance of the Parameters and Equation
Before being used to fit the four-dimensional predator-prey system, the experimental data of the culture pH and remaining concentrations of total nitrogen, biomass, and nisin [16][17][18] were smoothed using the following logistic equations (5-7): For the culture pH (Q(t) = pH(t)) and total nitrogen (Q(t) = TN(t)) decrease, we modified the logistic decline equation presented by Goudar et al. [28] as follows: For biomass (X(t)) production, the logistic equation presented by Goudar et al. [30] was used by considering that the death cell rate was zero, as follows: Being D = X max and E = X max −X o X o . For nisin (BT(t)) synthesis, we modified the logistic decline equation [30] as follows: Being G = BT max .
In this work, the values of the different constants were first obtained by numerical integration of the four differential equations (dpH(t)/dt, dTN(t)/dt, dX(t)/dt, and dBT(t)/dt) of the competitive predator-prey system (1)-(4), minimizing the sum of quadratic differences between equation-predicted and experimental values, with the non-linear least squares (quasi-Newton) method included in the Solver tool of Microsoft Excel 2016 spreadsheet. Then, the values and statistical significance (p < 0.05) of both the constants and the four-dimensional Lotka-Volterra-like equation were corroborated with an accurate statistical analysis using the statistics software SigmaPlot for Windows version 12.0 (Systat Software, Inc., San Jose, CA, USA, 2012) and the Regression (non-lineal) module of the software package IBM SPSS Statistics 20.0 for Windows (Release 20.0.1; SPSS Inc., Chicago, IL, USA, 2021). The significance of the different constants in each differential equation was evaluated by using Student's t-test and considering the corresponding p values [31] for each constant. Thus, the values of the constants with a high t-value and a p value lower than 0.05 were considered statistically significant. The significance of the different constants in each differential equation was evaluated by using Student's t-test and considering the corresponding p values for each constant. The convergence of the parameters was checked by using the Levenberg-Marquardt method [32,33]  The global consistency of the four differential equations was verified by using the overall Fisher's F-test (α = 0.05) and considering the corresponding p values for each equation. Thus, a high F-value and a p value lower than 0.05 means that the differential equation was statistically significant.
The goodness-of-fit of the four differential equations was also checked using the mean relative percentage deviation modulus (RPDM) values [5] where s is the number of experimental data, Yexp i is the experimental value and Ypred i is the value predicted by the equation. Values of R 2 ≥ 0.95, RPDM < 10% [5], and B f and A f close to 1 [34] indicate that the corresponding equation was accurately fitted to the experimental data.

Mathematical Analysis
The mathematical analysis of the four-dimensional Equation (Equations (1)-(4)) was performed with the MATLAB Runtime R2021a (MathWorks Inc., Natick, MA, USA) using the values of the constants corresponding to the unsupplemented culture of the DW-G series.

Conclusions
The main contribution of this paper is the development and mathematical analysis, for the first time, of a four-dimensional predator-prey system for an accurate description of the batch nisin production system by L. lactis CECT 539 in different series of fermentation in DW media supplemented with different concentrations of glucose (DW-G cultures), total sugars and phosphorous (DW-TS-TP cultures), or MRS broth nutrients (DW-MRS cultures).
The results obtained in this paper demonstrated that the microbiological bacteriocin production system could be explained by considering the biological approach used to describe the relationships between species that compete for the same nutrient sources. Thus, further knowledge is provided about the relationship between the main culture variables (culture pH, total nitrogen consumption, and the synthesis of biomass and bacteriocin) involved in nisin production, which is usually difficult to explain.
The mathematical novelty of this paper relies on the determination, for the first time, of the equilibrium points and trajectories around the stable equilibrium point obtained from the mathematical analysis of the four-dimensional model designed in this paper. Thus, the existence of an asymptotically stable equilibrium point with biological sense (pH = 4.665, TN = 0.301 g/L, X = 0.477 g/L, BT = 23.812 BU/mL), to which the culture evolved to reach the stationary state, was proved. The existence of this equilibrium point is in perfect agreement with the experimental point (pH Funding: Financial support from Xunta de Galicia (ED431C2018/42-GRC) and the European Regional Development Fund (ERDF) is gratefully acknowledged.
Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Data are available in the manuscript.

Conflicts of Interest:
The authors declare no conflict of interest.

Time (h) pH(t)
Culture pH value over the time pH 0 Initial culture pH value pH f Final culture pH a Intrinsic pH drop rate (h −1 ) b Quotient between the intrinsic pH drop rate and the theoretical minimum pH value for growth (h -1 ) c Constant that represents the effect of biomass production on pH time course (L/g/h) TN(t) Total nitrogen concentration (g/L) over the time TN 0 Initial total nitrogen concentration (g/L) TN f Final total nitrogen concentration (g/L) d Intrinsic TN consumption rate (h −1 ) e Quotient between the intrinsic TN consumption rate and the theoretical maximum TN concentration that biomass can consume (L/g/h) f Intrinsic TN consumption rate (L/g/h) for biomass production g Intrinsic TN consumption rate (mL/BU/h) for nisin production h Constant that represents the effect of pH time course on TN consumption (h −1 ) X(t) Biomass concentration (g/L) over the time X 0 Initial biomass concentration (g/L) X max Maximum biomass concentration (g/L) i Intrinsic growth rate (h −1 ) j Quotient between the intrinsic growth rate and the theoretical maximum biomass concentration that system can support (L/g/h) k Efficiency of TN utilization (mL/BU/h) to be channeled into cells of L. lactis rather than nisin (competition coefficient) l Constant that represents the effect of pH time course on the growth (h −1 ) m Constant that represents the effect of TN time course on the growth (L/g/h) BT(t) Nisin concentration (BU/mL) over the time BT max Maximum nisin concentration (BU/mL) n Intrinsic nisin production rate (h −1 ) o Quotient between the intrinsic nisin production rate and the theoretical maximum nisin concentration that biomass can produce (h −1 ) p Efficiency of TN utilization (mL/BU/h) to be channeled into nisin rather than cells of L. lactis (competition coefficient) q Constant that represents the effect of pH time course on nisin synthesis (h −1 ) r Constant that represents the effect of TN time course on nisin synthesis (L/g/h) A, B, C 0 , C 1 , C 2 Constants in Equation (5)  Predicted values by the corresponding equation for culture pH and concentration of total nitrogen, biomass, and nisin.

R 2
Correlation coefficient