# Mechanistic modeling of the dynamics of phage attack during milk acidification in the cheese-making process

Michèle Bou Habib, Emmanuel Bernuau, Benjamín-José Sánchez, Dominique

Swennen, Ahmad A. Zeidan, Cristian Trelea, Jannik Vindeloev

► To cite this version:

Michèle Bou Habib, Emmanuel Bernuau, Benjamín-José Sánchez, Dominique Swennen, Ahmad A. Zeidan, et al.. Mechanistic modeling of the dynamics of phage attack during milk acidification in the cheese-making process. Journal of Food Engineering, 2025, 387, pp.112329. 10.1016/j.jfoodeng.2024.112329. hal-04741360

HAL Id: hal-04741360

https://hal.science/hal-04741360v1

Submitted on 17 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
1 Highlights

2 Mechanistic modeling of the dynamics of phage attack during milk acidification in the cheese-
making process

4 Michèle Bou Habib,Emmanuel Bernuau,Benjamín José Sánchez,Dominique Swennen,Ahmad A. Zeidan,Ioan-Cristian

5 Trelea, Jannik Vindeloev

6 • Prediction of success/failure of milk acidification in the presence of phage

7 • Development and experimental validation of a dynamic mechanistic model

8 • Identification of key interactions to enhance phage attack interventions
# Mechanistic modeling of the dynamics of phage attack during milk acidification in the cheese-making process

Michèle Bou Habib<sup>a</sup>, Emmanuel Bernuau<sup>a,*</sup>, Benjamín José Sánchez<sup>b</sup>, Dominique Swennen<sup>a</sup>, Ahmad A. Zeidan<sup>b</sup>, Ioan-Cristian Trelea<sup>a</sup> and Jannik Vindelov<sup>b</sup>

## ARTICLE INFO

Keywords:

kinetic model

virulent phage

Lactococcus lactis

model validation

## ABSTRACT

Bacteriophage attacks represent a major threat in the dairy industry. Here, an unstructured mechanistic model predicting the dynamics of milk acidification in case of phage attack was developed and experimentally validated. Multiple acidification experiments were run with different combinations of initial phage titers and bacterial concentrations and the resulting pH dynamics were recorded. The model could successfully predict the success or failure of milk acidification. Using the model, important biological parameters were deduced from simple, low-cost acidification measurements. These parameters included bacteria's maximum growth and lysis rates, phages' burst size, etc. Sensitivity analysis helped identify biologically relevant aspects of phage-host interactions. Growth and lysis kinetics were shown to have the most important impacts. This knowledge can be used to develop easy routine strategies to fight phage attack in the dairy industry. The model can be used to raise awareness amongst cheese makers on the importance of cleaning to avoid food and material waste.

## 1. Introduction

Microorganisms play a crucial role in cheese production, both during the acidification and ripening steps (Marcó et al., 2012). During the former, lactic acid bacteria (LAB), such as *Streptococcus thermophilus*, *Lactococcus lactis*, and *Lactococcus cremoris*, are used to convert lactose to lactic acid, thus acidifying milk. This step largely impacts the final cheese composition and quality (Kongo, 2013). To control the fermentation and obtain high-quality end-products, an inoculation in the order of $10^7$ Colony Forming Units per mL (CFU/mL) of milk is typically required (Garneau and Moineau, 2011).

According to Moineau and Lévesque (2005), between 0.1% and 10% of industrial milk fermentations are negatively affected by bacteriophages (or phages), i.e., viruses that attack bacteria. In milk fermentation, phages slow down or even completely inhibit LAB growth, therefore preventing milk acidification. It is a common practice in dairies to use the same cheese vat repeatedly, with only limited cleaning or rinsing between fillings. The LAB culture is added during the filling of the cheese vat with milk to save time and better distribute the cells. After a set fermentation time, the content of the vat is emptied and processed further downstream Fox and McSweeney (2017). The vat is considered ready for another filling round after a quick rinse. During these successive fillings, phages may build up and cause fermentation failure. Whenever the fermentation is halted in this way, significant economic losses and environmental

<figure><img src="image_2.png" /></figure>

<footnote><sup>a</sup>Université Paris-Saclay, INRAE, AgroParisTech, UMR SayFood, Palaiseau, 91120, France</footnote>

<footnote><sup>b</sup>Novonesis, R&D, Hørsholm, 2970, Denmark</footnote>

<footnote>*Corresponding author</footnote>

<footnote>emmanuel.bernuau@agroparistech.fr (E. Bernuau)</footnote>
issues arise, as milk is discarded and the plant needs to be thoroughly cleaned (Brüssow, 2001). In less drastic cases, the production is slowed down, resulting in reduced cheese quality and longer production times (Ledeboer et al., 2002).

As shown in Figure 1, there are two types of phage attack cycles: the lytic and the lysogenic cycles. Virulent phages undergo exclusively lytic cycles, whereas temperate phages can undergo both cycles.

In the beginning of both cycles, the phage is adsorbed on the surface of the host bacterium. Afterward, the phage injects its genetic material into the bacterium. In the case of a lytic cycle, the phage's genome is replicated inside the bacterial cell, and after some time, known as the latency time, the infected cell is lysed and releases the newly created phages in a sudden burst. The released progeny is then ready to infect new cells, and the cycle is repeated. In the case of a lysogenic cycle, the phage's genome is integrated inside the bacteria's chromosome and forms what is known as a "prophage." The phage's genome stays and replicates within the bacteria's chromosome for multiple generations until the prophage enters the lytic cycle due to environmental stresses (Sinha et al., 2020). Most lactococcal phages found in failed industrial milk fermentations belong to either of three groups: the 936 and c2 phage groups, which are all virulent, and the P335 phage group, whose phages can be either temperate or virulent (Mahony et al., 2015).

Simulations are a valuable tool to gain deeper insight into the dynamics of phage attack (Santos et al., 2014). One approach is to use dynamic, unstructured mechanistic models, as they do not only simulate a process' temporal evolution but also help better understand the physical, chemical, and biological phenomena involved. On the one hand, this is accomplished by selecting the model structure that fits the experimental data best. On the other hand, the study of the model's sensitivity to parameter variations gives an indication of the most influential phenomena.

Campbell (1961) was the first to propose a model to describe the dynamics of phage attack. The model consisted of two differential equations describing the evolution of bacteria and phages. It included mechanisms for bacterial growth, phage adsorption, phage inactivation, and burst of progeny while considering the latency time between bacterial infection and lysis. This model was a starting point from which several other models have been derived by changing the expression of some terms or studying in more detail some mechanisms (Santos et al., 2014; Chaudhry et al., 2018; Beretta and Kuang, 1998; Bull et al., 2006; Levin et al., 1977). Key biological parameters such as adsorption rate, latent period, burst size, bacterial growth rate, and substrate uptake rate are common to all these models. Levin et al. (1977) proposed a generalized model to describe any number of species of susceptible bacteria. They modeled infected bacteria, growth-limiting substrates, and phages as well. Beretta and Kuang (1998) proposed a model where susceptible bacteria reproduce while infected bacteria are removed by lysis. A new parameter, lysis or death rate, was introduced. It is the inverse of the latent period and indicates the number of cells lysed per unit time. Finally, in some studies, models have been extended to include acquired bacterial resistance to phages (Santos et al., 2014; Chaudhry et al., 2018; Lenski and Levin, 1985; Cairns et al., 2009).
Among the few published phage attack models, most are in the field of phage therapy and disease epidemics; models applied to the fields of food science, wastewater, and bioremediation are still lacking. Additionally, only a few of the developed models have been validated experimentally. Among these experimentally validated models, such as Santos et al. (2014) and Levin et al. (1977), experiments were conducted under only one specific set of conditions and in only one replicate. Finally, the applicability to a range of practical situations and the robustness to variations in experimental conditions were not investigated.

The objective of this work was to develop and validate an unstructured mechanistic model to predict milk acidification by *Lactococcus lactis*, a key player in cheese fermentation, in the presence of a virulent phage.

In this paper, we start by explaining the experimental setup in Section 2. In Section 3, we show the experimental acidification results We then detail, in Section 4, the model development based on these experiments. Then, in Section 5, we focus on estimating key biologically meaningful parameters Finally, we conclude in Section 6.

## 2. Experimental Setup

### 2.1. Strains and milk

One *Lactococcus lactis* strain (strain A) and one *Ceduovirus Lactococcus* phage c2 (phage B), both obtained from the Chr. Hansen Culture Collection, were used to model phage-host dynamics. In this experimental setup, the bacterial strain was maintained as a frozen pellets at -20°C (or lower), while a phage stock was kept at 4°C after preparation of the lysate (Poulsen et al., 2019). Milk acidification experiments were carried out using B-milk, which is prepared using low-fat skim milk powder to a level of dry matter of 9.5% and subjected to heating at 99°C for 30 min, followed by cooling at 30°C. The milk was kept at 4°C until further use. B-milk was then supplemented with 50 mg of bromocresol purple salt and 50 mg of bromocresol green salt (both from Sigma Aldrich, St. Louis, Missouri, United States) as a pH color indicator, as previously described (Poulsen et al., 2019).

### 2.2. Experimental design and pH measurements

Acidification experiments were performed in a volume of 400µL using deep 96-well plates, and a 2D concentration-gradient of strain-to-phage was prepared using a Hamilton MicroLab platform (Hamilton, Bonaduz, Switzerland). The sensitive strain was supplied in frozen pellets to concentrations ranging from 10g/100L of dry weight (~ 10⁶ CFU/mL) to 80g/100L (~ 8.0 × 10⁶ CFU/mL), while the phage concentrations varied between 0 and 10⁷ Plaque Forming Units per mL (PFU/mL).

After inoculation of the sensitive strain and its phage, the plates were incubated for 20h on top of a flat-bed scanner (HP ScanJet G4010) at 30°C (mesophilic culture). pH-dependent changes in color were recorded every 4 min and converted into pH values using *pH Coprah* software (v.2.1.52; Videometer, Copenhagen, Denmark).
The data consisted of 2 separate runs, with 3 replicates for each run, for 48 acidification curves with 12 different initial phage ($P_0$) titers and 4 different initial susceptible bacteria ($X_{s0}$) concentrations in 3 replicates. The experiment is summarized in a conceptual diagram in Figure 2.

### 3. Experimental Acidification Data

Figure 3 shows an example of available data of runs Run_01 and Run_02 for an initial bacteria concentration of 40 g/100L and an initial phage titer of 0.

After inoculation, bacteria require time, known as lag time, to adapt to a new environment. Afterward, acidification starts. The initial pH of milk is 6.6. A steep decrease in pH indicates fast acidification for the first 4 hours where pH reached a value of 4.6. The acidification then slows down until it completely stops towards the end of the experiment. The minimum pH reached is around 4.2. This could be due either to substrate depletion or product inhibition. In this experiment, lactose is known to be in excess. Therefore, product, namely lactic acid, inhibition is assumed and should be considered while developing the model. For these initial conditions, results from both runs are in agreement. It also shows a pH difference in the order of 0.1 between the three samples from the same run.

Figure 4 shows the results of the 48 experiments. Each curve is similar to the one in Figure 3: time varies on the x-axis between 0 and 20 hours, and pH varies on the y-axis between a minimum of 4.2 and a maximum of 6.6. Each column represents an initial phage titer in PFU/mL, and each line represents an initial bacteria concentration in g/100L. In the first column, no phages are added to the milk. As the initial concentration of bacteria decreases, acidification is delayed, but the same final pH is reached.

Curves in Figure 4 can be divided into three distinct zones. The blue zone represents successful acidification, the final pH reached for all replicates is lower than 4.5. The gray zone represents the transition phase between success and failure. In this zone, the final pH is between 4.5 and 5.5. Finally, the red zone indicates the complete failure of acidification, all final pH values are higher than 5.5.

It can be seen in the transition zone that there are discrepancies between the final pH values and the curve shapes between replicates of the same repetition. For each set of initial conditions, the standard deviation between pH measurements at each time was calculated (i) separately for the 3 replicates of each run and (ii) including the 6 measurements from both runs. Then, the conditions were grouped based on the final outcome: success, transition, and failure, and the standard deviations were averaged by zone.

The results are represented in Table 1. The standard deviations are higher in the transition zone compared to the success and failure zones for both runs, indicating more variability in the transition zone. This can be due to various possible factors such as experimental variabilities on the initial conditions or lack of homogeneity within the samples. It suggests a high sensitivity of the phenomena on initial conditions. Modeling should be used to quantify the behaviors in the
transition zone. In some cases, a slight increase in pH is observed such as for the initial conditions $X_{s0} = 10$ g/100L and $P_0 = 1.6 \times 10^5$ PFU/mL. It is thought to be linked to the release of amino acids from lysed bacteria making the milk more alkaline.

## 4. Model Development

The model consists of 5 ordinary differential equations with 11 parameters ($\theta$). The susceptible ($X_s$) and infected ($X_i$) bacteria, lactose ($L$), and lactate concentrations ($A_H$), and phage titer ($P$) are grouped in a vector of state variables $X(t, \theta)$, which is a solution of the set of ordinary nonlinear differential equations:

$$ \frac{dX}{dt} = f(t, X, \theta), \qquad (1) $$

where $\theta \in \Theta$ is the vector of model parameters and $\Theta$ is the feasible parameter space. The pH is the only observable of the experiments, and is denoted as $y$, such that:

$$ y = g(X), \qquad (2) $$

where $g(X)$ is a vector describing the relationship between lactic acid and pH and is further detailed in equations (10) and (11) in Appendix A. The optimal set of parameters $\hat{\theta}$ was estimated using the least squares method. This was done using Python version 3.8 and the `least_squares` function from the `scipy optimize` subpackage, version 1.10.1. The function was given another function that computes the vector of residuals between the experimental and simulated calibration data. Other arguments were left to default, such as method= "trf" for the Trust Region Reflective algorithm (Branch et al. (1999)), the tolerance for termination by the change of the cost function ftol=1e-8, and tolerance for termination by the change of the independent variable xtol=1e-8.

Figure 5 is a graphical representation of the model. Table 2 shows the symbols of state variables and parameters. Detailed description of considered phenomena and model equations are given next.

### 4.1. Bacterial Concentrations

The growth of susceptible LAB (not yet infected by phages) was assumed to follow Monod-like growth kinetics (Monod, 1949). In milk fermentation, the substrate, i.e., lactose, limits bacterial growth in the Monod growth law, hence the presence of the half-saturation constant parameter $K_L$ (Equation 4). The Monod half-saturation constant $K_L$ is the substrate concentration at which the growth rate $\mu$ is half of the value of $\mu_{max} \frac{L}{K_L+L}$. Growth is also inhibited by the product, namely lactic acid. We adopted the model of Luong (1985) to represent this inhibition. In addition, when bacteria are first put in a particular environment, they need an initial period to adapt. We adopted a modified
form (García et al., 2017) of Baranyi's model (Baranyi and Roberts, 1994) to represent this lag term (Equation 4).

When phages are present, they attack bacteria by first being adsorbed on their surface, as shown in Figure 1. This was represented in Equation 3, where the adsorption term is a product of the concentration of susceptible bacteria, the phage titer, and the average adsorption rate $k_A$. Thus, the dynamics of bacterial growth and growth rate are represented by Equations 3 and 4, respectively:

$$ \frac{dX_s}{dt} = \mu X_s - k_A X_s P \quad (3) $$

$$ \mu = \mu_{\max} \frac{L}{L + K_L} \cdot \frac{a_0}{a_0 + (1 - a_0)e^{-\mu_{\max}t}} \left( 1 - \left( \frac{A_H}{A_{H_{\max}}} \right)^{\gamma} \right) \quad (4) $$

Susceptible bacteria ($X_s$) attacked by phages become infected bacteria ($X_i$) and subsequently die with a specific lysis rate $\lambda$. Equation 5 describes the dynamics of infected bacteria:

$$ \frac{dX_i}{dt} = k_A X_s P - \lambda X_i \quad (5) $$

## 4.2. Lactose metabolism to Lactic Acid

The dynamics of consumption of lactose and production of lactic acid are assumed to be proportional to the growth of bacteria. Yield coefficients of substrate consumption per biomass produced ($Y_{L/X}$) and product formation per substrate consumed ($Y_{A_H/L}$) are the proportionality factors. In the case of phage attack, infected bacteria are considered to compete with susceptible bacteria for substrate consumption (Santos et al., 2014). The dynamics of lactose utilization and lactic acid production are represented by Equations 6 and 7, respectively:

$$ \frac{dL}{dt} = -Y_{L/X} \mu (X_s + X_i) \quad (6) $$

$$ \frac{dA_H}{dt} = Y_{A_H/L} Y_{L/X} \mu (X_s + X_i) \quad (7) $$

## 4.3. Phage Infection

Phage multiplication is assumed to be proportional to the lysis of infected bacteria with $\beta$, the burst size, as the proportionality coefficient. Phages are removed from the medium when adsorbed. Phage deactivation is neglected because the experiments occur for a relatively short time. Equation 8 describes the evolution of phage titer:
$$ \frac{dP}{dt} = -\alpha k_A X_s P + \beta \lambda X_i \qquad (8) $$

## 4.4. Calibration and analysis of model subparts

Initially, model parameters were chosen based on the literature and expert knowledge. Then, the model was optimized by running a two-step optimization: first, an estimation of non-phage-related parameters, and then the estimation of the remaining parameters. The estimations were done by minimizing the sum of discrepancies between the experimental and simulated data, and the optimal set of parameters $\hat{\theta}$ was computed. The simulated pH was deduced from lactic acid concentration using the relation described in the Appendix A. Experimental data was divided into calibration data to estimate the parameters and validation data to independently assess model performance in a checkerboard pattern with respect to Figure 4 (if both the column and row numbers were even or odd, the initial condition was part of the validation data, otherwise it was part of the calibration data).

To begin with, a sub-model without phages was considered. Equations 5 and 8 were not taken into account and the phage adsorption term in Equation 3 was omitted. Among the 7 non-phage related parameters listed in Table 1, 2 parameters were fixed as follows. The yield coefficients $Y_{L/X}$ and $Y_{A_H/L}$ were highly correlated and could not be individually identified with the available data. $Y_{A_H/L}$ was fixed to a value of 0.8 (g/L) $(g/L)^{-1}$ that is commonly found in the literature ((Cachon and Divies, 1994), (Parente et al., 1994)). The substrate, lactose, was not limiting. The parameter $K_L$ was then fixed to 0.5 g/L (Zacharof and Lovitt (2013)). Once the remaining 5 non-phage-related parameters of the sub-model were estimated, they were fixed, and the 4 phage-related parameters listed in Table 1 were estimated.

The initial conditions of phage titer and bacterial concentration were set as shown in Figure 4. The milk had an initial lactose concentration of 46 g/L. The initial lactic acid concentration was considered to be negligible. No infected bacteria were present at the beginning of the experiment.

## 4.5. Sensitivity analysis

To better understand the role of each phenomenon on the acidification dynamics, the sensitivity of final pH with respect to model parameters was computed for different initial conditions. The relative sensitivity of pH with respect to a parameter p was computed using Equation 9 from Cacuci et al. (2005):

$$ \text{relative sensitivity} = \left( \frac{pH_{\text{up}} - pH_{\text{down}}}{pH} \right) / \left( \frac{pH_{\text{up}} - pH_{\text{down}}}{p} \right) \qquad (9) $$
where $p_{\text{down}}$ and $p_{\text{up}}$ are the values of the parameter decreased and increased by 5% respectively, and $pH_{\text{down}}$ and $pH_{\text{up}}$ are the values of final pH for $p_{\text{down}}$ and $p_{\text{up}}$, respectively. A variation of 5% was selected to minimize potential numerical issues and to ensure that the observed variation was due to the parameter change rather than numerical noise. A similar analysis was done with 1% change in parameters, and the results were consistent with those obtained using a 5% change.

# 5. Calibration, Validation, and Discussion

## 5.1. Parameter Estimation

The optimal set of parameters $\hat{\theta}$ was estimated. The relative standard deviation (RSD) of each parameter was calculated. Table 3 summarizes the values of these parameters, their units, and their RSD.

The RSD was less than 2% for all estimated parameters. In general, parameter values estimated were within the ranges of values found in the literature. However, the values of the parameters vary depending on the specific strain, environmental conditions, and other factors.

For *Lactococcus lactis*, $\mu_{\max}$ values in other studies range from 0.2 to 1.2 $h^{-1}$ (Parente et al., 1994; Boonmee et al., 2003; Zacharof and Lovitt, 2013; Chen et al., 2015). The value found using our model is within the same order of magnitude as values in the literature. The value of $\mu_{\max}$ depends on the specific strain, the culture medium, and the conditions in which the bacteria are grown. Parente et al. (1994) and Zacharof and Lovitt (2013) found a yield of L. lactis growth over lactose consumption ($Y_{L/X_s}$) ranging from 0.04 to 0.09 g.(g/100L)$^{-1}$, which is close to our value. Cachon and Divies (1994) found a maximum lactic acid concentration ($A_{H_{\max}}$) at uncontrolled pH of 9 g/L which is again compatible with our value. Lactoccal phages' burst size ($\beta$) values range between 10 and 250 particles per cell (Müller-Merbach et al., 2007), i.e. $10$ to $250 \times 10^5$ (PFU/mL).(g/100L)$^{-1}$. The cell lysis rate ($\lambda$) is comparable to the inverse latency time between cell infection and burst. On average, phages have a latency period of 40 to 60 minutes (Müller-Merbach et al., 2007). Therefore, $\lambda$ takes values between 1 and $1.5 \ h^{-1}$. The adsorption rate $k_A$ could not be compared to literature values due to unit incompatibility. The parameter $a_0$ can take values from 0 to 1 (García et al., 2017). The estimated value of 0.01 means that at the beginning of the fermentation, the growth is slow.

## 5.2. Model Validation

Figure 4 shows the model and the experimental data for pH for all 48 conditions. Globally, the model captured the major phenomena and mechanisms observed. It accurately delineated the three colored zones: acidification success, failure, and transition. It correctly predicted the shape of the curves in the blue (success) zone. There was high variability between the experimental acidifications in the gray (transition) zone. However, the model mostly fell between the
different experimental curves. Finally, the model predicted no acidification in the red (failure) zone. No mechanism to explain pH increase was included in the model, and no pH increase was predicted, as expected.

Besides pH, the model predicted the dynamics of the state variables: susceptible and infected LAB, lactose, lactic acid, and phage. As an example, Figure 6 shows the predicted dynamics of these variables, and susceptible bacteria's growth, for the case with an initial LAB concentration of 80 g/100L and an initial phage titer of $10^4$ PFU/mL. Susceptible LAB grew to reach a maximum value of 180 g/100L, and when phages took over, their concentration decreased to become almost null towards the end of the experiment. In parallel, the concentration of infected LAB increased mostly after the susceptible bacteria's peak. Initially, phage titer decreased quickly due to adsorption on susceptible bacteria. It then increased due to phage burst by cell lysis. When the susceptible bacteria's concentration became null, the phage titer remained constant since no phage death mechanism was considered. At the end of the experiment, around 12 g/L of lactose were consumed and 9.5 g/L of lactic acid were produced. Susceptible bacteria's growth initially increased slowly due to the latency factor. It then increased rapidly until it reached a peak at 2.5 hours, after which it decreased due to lactic acid inhibition. Its largest value (≈ 0.3 h⁻¹) never approached $\mu_{max} = 1.88$ h⁻¹ (Table 3), confirming that the significance of parameters (here, $\mu_{max}$) is model-dependent.

## 5.3. Model use to assess acidification in case of phage attack

The final pH reached is considered a good indicator of the success or failure of mild acidification. The initial pH of milk is typically around 6.6. In the absence of phages, the final pH reached is 4.3, as acidification stops due to acid inhibition. In the case of phage attack, milk acidification can still be considered successful when the final pH is between 4.3 and 5.0. If the final pH is between 5.0 and 6.0, this is the transition zone between normal and no acidification. Finally, a final pH above 6.0 indicates failure of milk acidification due to phage attack.

Several simulations of the developed model were run, with different combinations of initial LAB concentrations and phage titers in the ranges of [10,80] (g/100L) and $[10^4, 10^7]$ (PFU/mL), respectively, where the model was validated. Figure 7a represents the final pH simulated using the model as a function of the initial conditions of bacteria and phages.

This experiment is a simplified real-life milk acidification experiment, where only one LAB strain and one phage interact. In such a case, if the initial concentration of LAB in the milk is known and if the phage contamination can be assessed, this graph can be used to predict the success or failure of acidification. If the final pH falls in the blue zone, it is safe to continue the cheese-making process. If it falls in the transition zone, which has been seen to be highly sensitive to the variability in initial conditions, special attention should be given to ensure that milk acidification will occur. However, if the pH falls in the red zone, it is advisable to stop the production process and clean the platform.
It can be seen in Figure 4 that there are discrepancies between the different replicates having theoretically the same initial conditions. This is most obvious in the transition (gray) zone. These discrepancies might be justifiable if small variations in the initial conditions have a large impact on the change of pH over time. To assess this hypothesis, the calibrated model was used to study the sensitivity of the final pH to the initial phage titer. Figure 7b shows the sensitivity analysis results. The model is clearly most sensitive in the transition zone. For example, for the initial conditions $X_{s0} = 10$ g/100L and $P_0 = 8 \times 10^4$ PFU/mL, the relative sensitivity of final pH to $P_0$ is 0.125. In this case, the model predicts a final pH of 5. A change of +50% and -50% in $P_0$ (likely to happen experimentally) leads to a final pH of 5.6 and 4.3, respectively, which is congruent with the experimental values. This high sensitivity supports the explanation of the observed experimental variability in the transition zone as being due to the initial phage titer variability. However, this variation could also be due to other factors, such as differences in the interaction between phages and bacteria. In the transition zone, the infection might be subject to biological variability, allowing bacteria in some repetitions/runs to escape the phages while others do not.

## 5.4. Sensitivity analysis of model parameters

Sensitivity analyses were performed to understand the dynamics of phage attack and the phenomena that govern it. Sensitivities of final pH (after 20 hours) to model parameters were computed for five different initial conditions marked in Figure 4: one in the success zone (a), three in the transition zone (b1, b2, and b3), and one in the failure zone (c). Figure 8 shows the results grouped per parameter. A positive value of relative sensitivity indicates that an increase in a parameter results in a corresponding increase in the final pH, whereas a negative value implies that an increase in a parameter leads to a decrease in the final pH.

For condition (a) in the success zone, the final pH was mostly sensitive to $A_{Hmax}$. This result is to be expected as the reason behind acidification halting in this zone is acid inhibition and is not related to the presence of phages. The effect of acid inhibition was reduced in the transition and failure zones.

In the transition zone (for conditions (b1), (b2) and (b3)), four parameters ($\mu_{max}$, $\alpha$, $\beta$, and $\lambda$) have the highest impacts. An increased growth rate (higher $\mu_{max}$) led to more lactose being metabolized to lactic acid by susceptible bacteria, reducing pH. Parameter $\alpha$ represents the number of phages disappearing as they get adsorbed on a bacterial cell, whether they infect it or not. A higher $\alpha$ means more infectious phages are disappearing, slowing the infection down. Parameter $\beta$, burst size, represents the number of phages released after cell lysis. As this burst size increases, more bacteria-infecting phages are released, which accelerates infection and increases the final pH. The parameter $\lambda$ plays a role in two phenomena. On the one hand, an increase in the rate of cell lysis results in the faster elimination of infected cells that produce lactic acid, ultimately leading to a higher final pH. On the other hand, a higher cell lysis
rate leads to a quicker release of phages, targeting more susceptible bacteria and stopping lactic acid production which also results in a higher final pH.

For condition (c), i.e., in the failure zone, the sensitivities of the final pH to model parameters were very small or null. A higher μmax (an increased growth), or a smaller λ (a slower cell lysis), led to a lower final pH.

This analysis gives guidelines to shift the transition zone while expanding the success zone and reducing the failure zone. Based on this model, adsorption kinetics ($k_A$) have shown to be much less important than burst size ($β$) and lysis rate ($λ$) under selected initial conditions. In conclusion, phage-resistant LAB should be selected accordingly.

## 5.5. Model limitations and future developments

Although the kinetic model developed in this study was able to properly predict the outcome of a number of contamination scenarios, and provide insights into process parameters and strain selection strategies, it has important limitations that should be considered when interpreting the results. First, as any unstructured kinetic model built from ordinary differential equations, it assumes that the tank where the milk acidification occurs is well mixed, and no concentration gradients exist across it. This may not hold true in the case of our experimental setup, where sedimentation of bacteria has occasionally been observed. This is also not the case in cheese production, as mixing is only done at the beginning of the process, whereas most of the fermentation occurs in a solid matrix where there is inherently uneven nutrient distribution. Secondly, we assume Monod-like kinetics for growth, where the main variables controlling growth are the concentrations of the substrate (lactose) and the product (lactate). However, other variables that play an important role on determining the growth rate, such as metabolic bottlenecks, temperature, osmotic pressure, and other limiting nutrients (e.g., casein peptides), are not accounted for; this is especially relevant in large-scale fermentations where some of these parameters do not remain constant. Finally, the developed model represents the attack of a single type of phage on a single bacterial strain. This is normally not the case in cheese-making, where often multiple types of phages can infect a microbial culture containing a community of different strains. A model that would allow for multiple types of phages and bacteria should be expected to behave differently from the model developed in this study, as not only it would simulate how the different phages interact with the different bacteria, but also how different phages compete for the same bacterial strain, and how different bacteria interact among each other.

Going forward, it is possible to extend our method to quantify the impact of temperature, milk composition, and additives such as rennet, on the dynamics of phage-host interactions. Moreover, this setup has the potential to characterize new host-phage systems in an inexpensive and fast manner. Done systematically and at large scale, such characterization may be used to classify many phage-host pairs to provide new insights about general phage virulence. In addition, the use of non-invasive methods such as spectrophotometry, could also help quantify more state variables, which would provide the foundation for even more advanced models of phage-host dynamics. Combined with modern
molecular biology and bioinformatics tools, it may even be possible to identify biomarkers in the host and phage that predicts phage virulence. Such biomarkers could be used to select more resilient strains for the dairy industry. We expect the model developed in this study to play a significant role in designing strategies for mitigating phage attacks, at the strain, community, and application levels. At the strain level, sensitivity analysis can identify which aspects of phage-host interactions should be addressed to make more robust strains. At the microbial community level, once this model is extended to account for other strains, it should be possible to find a blend of strains that minimize the risk of failure on the culture level if one or more strains are hit by phages. Finally, at the application level, the model can be used to devise strategies that mitigate the risk of fermentation failure in successive fillings, such as increasing rinsing (i.e. reducing the amount of phage left in the vat), late inoculation (i.e. after re-filling the tank with fresh milk), or reducing fermentation time (to limit phage build-up). These model-driven strategies can result in significant saving for the cheese-making industry.

## 6. Conclusion

This study contributes to the development and validation of dynamic phage attack models. We used a high throughput method to systematically create many inhibition trials in microtiter plates of different combinations of a host and its phage. Our model was validated on a broad range of initial conditions, and the experimental dataset included many replicates. Previous models have not been validated with such an extensive dataset, making our approach unique. The results were used to estimate bacteria- and phage-related parameters, such as specific growth rate, adsorption rate, latency time, and burst size, by fitting our model predictions to milk acidification curves.

This setup offers a fast and cost-effective method to characterize new host-phage systems. By integrating this model with molecular biology and bioinformatics tools, it could guide the selection of more resilient strains in the dairy industry. The model also provides strategies to mitigate phage attacks that could result in significant savings for cheese-making by reducing the risk of fermentation failure.

## A. Appendix: Lactic Acid to pH calibration curve

In this study, pH is the only experimentally measured output. In order to link pH to lactic acid concentration, a state variable of the model, a direct ion equilibrium-based relation is not possible, because of the pH-dependent buffer effect of milk. Instead, a lactic acid titration in milk was generated and the data was fitted to a model proposed by Vereecken and Van Impe (2002):

$$Z = A_H - \frac{b_1}{2(b_1 - c_1)} \left[ (A_H + b_1) - \sqrt{(A_H + b_1)^2 - 4(b_1 - c_1)A_H} \right] \quad (10)$$
$$pH = \frac{1}{2c_2} \left[ A(b_2 - 2c_2)Z - \sqrt{A^2 b_2^2 Z^2 + 4Ac_2b_2^2Z} \right] + pH_{ini} \quad (11)$$

where Z is an intermediate variable and $A_H$ is the concentration of total dissociated and undissociated lactic acid in mmol/L. Graphical interpretation of parameters A, $b_1$, $b_2$, $c_1$, and $c_2$ can be found in Vereecken and Van Impe (2002). Table 4 shows the values of the estimated parameters from equations 10 and 11. Once the parameters were determined, they remained unchanged throughout the rest of the study. Squared average deviation was chosen as a measurement of the goodness of fit. It represents the square root of the average of squared differences between simulated and experimental data. Its value was 0.03 pH units. Figure 9 shows the fitting of the model equations to experimental data.

To acquire the data, increasing amounts of lactic acid were added to B-milk until a pH of approximately 4.0 was reached. Milk was placed at room temperature for 3 hours. Then, 100mL of milk was transferred to a beaker, where it was stirred with a magnetic stirrer for the duration of the experiment. The pH-meter was calibrated before use. Using a 10% solution of lactic acid, pH was stepwise decreased starting with 100 µL intervals of added acid. When an increasing buffering activity was observed, the interval of added acid was adjusted conveniently. The pH was measured and plotted in relationship to the amount of added acids.

## Declaration of competing interest

Benjamín José Sánchez, Ahmad A. Zeidan, and Jannik Vindelov were employed by Chr. Hansen A/S, a global supplier of bacterial cultures for the dairy industry, at the time of writing this manuscript. The authors' views presented in this manuscript, however, are solely based on scientific grounds and do not reflect the commercial interests of their employer.

## Data availability

Experimental Acidification data can be found here: https://doi.org/10.57745/W3ZDG2.

## Acknowledgements

This paper is based on the experimental data of acidification curves obtained by the method developed and optimized by Louise Gammeltoft, Mette Schultz, Bjørn Serritzlev, Josue Castro, and Jeorgos Trihaas from Chr. Hansen A/S, Hørsholm, Denmark. We are grateful for their expertise and assistance in this project. In addition we would like to thank Ditte Ellegaard Christiansen for preparing the phages, Maiken Bjerrum for preparing the milk titration, and Thomas Janzen for the selection of phages and discussions.
Funding

379

This work is part of the E-MUSE project which was funded by the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement N° 956126.

381

Credit author statement

382

**Michèle Bou Habib:** Conceptualization, Software, Formal Analysis, Methodology, Validation, Writing original draft. **Emmanuel Bernueau:** Conceptualization, Software, Formal analysis, Methodology, Review & Editing. **Ioan-Cristian Trelea:** Conceptualization, Formal analysis, Methodology, Review & Editing. **Benjamín-José Sánchez:** Conceptualization, Formal Analysis, Methodology, Review & Editing. **Dominique Swennen:** Conceptualization, Funding acquisition, Review & Editing. **Ahmad Zeidan:** Funding acquisition, Review & Editing. **Jannik Vindeloev:** Conceptualization, Data curation, Methodology, Review & Editing.

383

384

385

386

387

388

References

389

Baranyi, J., Roberts, T.A., 1994. A dynamic approach to predicting bacterial growth in food. International journal of food microbiology 23, 277-294.

390

Beretta, E., Kuang, Y., 1998. Modeling and analysis of a marine bacteriophage infection. Mathematical biosciences 149, 57-76.

391

Boonmee, M., Leksawasdi, N., Bridge, W., Rogers, P.L., 2003. Batch and continuous culture of lactococcus lactis nz133: experimental data and model development. Biochemical engineering journal 14, 127-135.

392

393

Branch, M.A., Coleman, T.F., Li, Y., 1999. A subspace, interior, and conjugate gradient method for large-scale bound-constrained minimization problems. SIAM Journal on Scientific Computing 21, 1-23.

394

395

Brüssow, H., 2001. Phages of dairy bacteria. Annual Reviews in Microbiology 55, 283-303.

396

Bull, J., Millstein, J., Orcutt, J., Wichman, H., 2006. Evolutionary feedback mediated through population density, illustrated with viruses in chemostats. The American Naturalist 167, E39-E51.

397

398

Cachon, R., Divies, C., 1994. Generalized model of the effect of ph on lactate fermentation and citrate bioconversion in lactococcus lactis ssp. lactis biovar. diacetylactis. Applied microbiology and biotechnology 41, 694-699.

399

400

Cacuci, D.G., Ionescu-Bujor, M., Navon, I.M., 2005. Sensitivity and uncertainty analysis, volume II: applications to large-scale systems. CRC press.

401

402

Cairns, B.J., Timms, A.R., Jansen, V.A., Connerton, I.F., Payne, R.J., 2009. Quantitative models of in vitro bacteriophage-host dynamics and their application to phage therapy. PLoS Pathogens 5, e1000253.

403

404

Campbell, A., 1961. Conditions for the existence of bacteriophage. Evolution, 153-165.

405

Chaudhry, W.N., Pleška, M., Shah, N.N., Weiss, H., McCall, I.C., Meyer, J.R., Gupta, A., Guet, C.C., Levin, B.R., 2018. Leaky resistance and the conditions for the existence of lytic bacteriophage. PLoS biology 16, e2005971.

406

407

Chen, J., Shen, J., Ingvar Hellgren, L., Ruhdal Jensen, P., Solem, C., 2015. Adaptation of lactococcus lactis to high growth temperature leads to a dramatic increase in acidification rate. Scientific reports 5, 14199.

408

409

Fox, P.F., McSweeney, P.L., 2017. Cheese: an overview. Cheese, 5-21.

410
411 García, M.R., Cabo, M.L., Herrera, J.R., Ramilo-Fernández, G., Alonso, A.A., Balsa-Canto, E., 2017. Smart sensor to predict retail fresh fish quality under ice storage. Journal of food engineering 197, 87–97.

412 Garneau, J.E., Moineau, S., 2011. Bacteriophages of lactic acid bacteria and their impact on milk fermentations. Microbial cell factories 10, 1–10.

413 Kongo, J.M., 2013. Lactic acid bacteria as starter-cultures for cheese processing: past, present and future developments. Lactic acid bacteria-R & D for food, health and livestock purposes, 1–22.

414 Ledeboer, A., Bezemer, S., De Haard, J., Schaffers, I., Verrips, C., Van Vliet, C., Düsterhöft, E.M., Zoon, P., Moineau, S., Frenken, L., 2002. Preventing phage lysis of lactococcus lactis in cheese production using a neutralizing heavy-chain antibody fragment from llama. Journal of dairy science 85, 1376–1382.

415 Lenski, R.E., Levin, B.R., 1985. Constraints on the coevolution of bacteria and virulent phage: a model, some experiments, and predictions for natural communities. The American Naturalist 125, 585–602.

416 Levin, B.R., Stewart, F.M., Chao, L., 1977. Resource-limited growth, competition, and predation: a model and experimental studies with bacteria and bacteriophage. The American Naturalist 111, 3–24.

417 Luong, J., 1985. Kinetics of ethanol inhibition in alcohol fermentation. Biotechnology and bioengineering 27, 280–285.

418 Mahony, J., Tremblay, D.M., Labrie, S.J., Moineau, S., van Sinderen, D., 2015. Investigating the requirement for calcium during lactococcal phage infection. International Journal of Food Microbiology 201, 47–51.

418 Marcó, M.B., Moineau, S., Quiberoni, A., 2012. Bacteriophages and dairy fermentations. Bacteriophage 2, 149–158.

419 Moineau, S., Lévesque, C., 2005. Control of bacteriophages in industrial fermentations. Bacteriophages: biology and applications, 285–296.

420 Monod, J., 1949. The growth of bacterial cultures. Annual review of microbiology 3, 371–394.

421 Müller-Merbach, M., Kohler, K., Hinrichs, J., 2007. Environmental factors for phage-induced fermentation problems: replication and adsorption of the lactococcus lactis phage p008 as influenced by temperature and ph. Food microbiology 24, 695–702.

422 Parente, E., Ricciardi, A., Addario, G., 1994. Influence of ph on growth and bacteriocin production by lactococcus lactis subsp. lactis 14onwc during batch fermentation. Applied Microbiology and Biotechnology 41, 388–394.

423 Poulsen, V.K., Derkx, P., Oregaard, G., 2019. High-throughput screening for texturing lactococcus strains. FEMS microbiology letters 366, fnz001.

424 Santos, S.B., Carvalho, C., Azeredo, J., Ferreira, E.C., 2014. Population dynamics of a salmonella lytic phage and its host: implications of the host bacterial growth rate in modelling. PloS one 9, e102507.

425 Sinha, S., Grewal, R.K., Roy, S., 2020. Modeling phage-bacteria dynamics, in: Immunoinformatics. Springer, pp. 309–327.

426 Vereecken, K.M., Van Impe, J.F., 2002. Analysis and practical implementation of a model for combined growth and metabolite production of lactic acid bacteria. International Journal of Food Microbiology 73, 239–250.

427 Zacharof, M.P., Lovitt, R.W., 2013. Modelling and simulation of cell growth dynamics, substrate consumption, and lactic acid production kinetics of lactococcus lactis. Biotechnology and bioprocess engineering 18, 52–64.
Predicting milk acidification

<figure><img src="image_3.png" /><figcaption>**Figure 1:** The lytic and lysogenic cycles of phage attack, after Sinha et al. (2020). Created with BioRender.com.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2:** Conceptual diagram of the experiment. Acidification experiments were performed in a volume of 400L using deep 96-well plates, and a 2D concentration gradient of strain-to-phage. After inoculation of the sensitive strain and its phage, the plates were incubated for 20h on top of a flat-bed scanner at 30°C. pH-dependent changes in color were recorded every 4 min and converted into pH. Created with BioRender.com.</figcaption></figure>
Predicting milk acidification

<figure><img src="image_5.png" /><figcaption>**Figure 3:** Acidification curves for $P_0$ = 0 PFU/mL and $X_{50}$ = 40 g/100L. Experiments were done in two separate runs 'Run_01' and 'Run_02', each with three repetitions.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4:** Acidification curves for all 48 conditions of phage titer and LAB concentration. Each curve has the same scale as Figure 3. The x-axis represents time, ranging from 0 to 20 hours, and the y-axis represents pH, ranging from 4 to 6.6. The black curve is the result of the model simulation after parameter fitting (section 4.4). a, b1, b2, b3, and c are chosen conditions for the sensitivity analysis in 5.4. Background: blue: success, gray: transition zone, red: failure.</figcaption></figure>

<table><thead><tr><th></th><th>Success zone</th><th>Transition zone</th><th>Failure zone</th></tr></thead><tbody><tr><td>Run_01</td><td>0.030</td><td>0.051</td><td>0.009</td></tr><tr><td>Run_02</td><td>0.038</td><td>0.055</td><td>0.013</td></tr><tr><td>Both Runs</td><td>0.065</td><td>0.186</td><td>0.039</td></tr></tbody></table>

**Table 1**

Values of the average standard deviations per zone. For each set of initial conditions, the standard deviation between pH measurements at each time was calculated (i) separately for the 3 replicates of each run and (ii) including the 6 measurements from both runs. Then, the conditions were grouped based on the final outcome: success, transition, and failure, and the standard deviations were averaged by zone.
<figure><img src="image_7.png" /><figcaption>**Figure 5:** Graphical representation of the phage attack model. Symbols are all explained in Table 2. Created with BioRender.com.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 6:** Model dynamics predictions for $X_{s0} = 80$ g/100L and $P_0 = 10^4$ PFU/mL</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>(a) Simulated final pH as a function of initial conditions</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) Relative sensitivity on final pH on initial phage titer as a function of initial conditions</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>**Figure 8:** Final pH sensitivity to model parameters for 5 different initial conditions marked in Figure 4</figcaption></figure>
<figure><img src="image_12.png" /><figcaption>**Figure 9:** Conversion of total lactic acid to pH. The titration was done by adding 10% solution of lactic acid to standardized B-milk.</figcaption></figure>

<table><thead><tr><th>Category</th><th>Symbol</th><th>Definition</th><th>Unit</th></tr></thead><tbody><tr><td></td><td>t</td><td>Time</td><td>h</td></tr><tr><td rowspan="5">State Variables</td><td>A<sub>H</sub></td><td>Lactic Acid Concentration</td><td>g/L</td></tr><tr><td>L</td><td>Lactose Concentration</td><td>g/L</td></tr><tr><td>X<sub>s</sub></td><td>Susceptible Bacteria Concentration</td><td>g/100L</td></tr><tr><td>X<sub>i</sub></td><td>Infected Bacteria Concentration</td><td>g/100L</td></tr><tr><td>P</td><td>Phage Titer</td><td>PFU/mL</td></tr><tr><td rowspan="7">Non phage-related parameters</td><td>μ<sub>max</sub></td><td>Theoretical Maximal Growth Rate</td><td>h<sup>-1</sup></td></tr><tr><td>Y<sub>A<sub>H</sub>/L</sub></td><td>Yield of Product over Substrate</td><td>(g/L)·(g/L)<sup>-1</sup></td></tr><tr><td>Y<sub>L/X</sub></td><td>Yield of Substrate Consumption over Biomass Growth</td><td>(g/L)·(g/100L)<sup>-1</sup></td></tr><tr><td>γ</td><td>Shape Factor</td><td>-</td></tr><tr><td>K<sub>L</sub></td><td>Half-saturation Constant</td><td>g/L</td></tr><tr><td>A<sub>Hmax</sub></td><td>Maximal Acid Concentration</td><td>g/L</td></tr><tr><td>a<sub>0</sub></td><td>Lag Parameter</td><td>-</td></tr><tr><td rowspan="4">Phage-related parameters</td><td>α</td><td>Phage to Bacteria Ratio</td><td>(PFU/mL)·(g/100L)<sup>-1</sup></td></tr><tr><td>k<sub>λ</sub></td><td>Adsorption Rate</td><td>(PFU/mL)<sup>-1</sup>·h<sup>-1</sup></td></tr><tr><td>β</td><td>Burst Size</td><td>(PFU/mL)·(g/100L)<sup>-1</sup></td></tr><tr><td>λ</td><td>Cell Lysis Rate</td><td>h<sup>-1</sup></td></tr></tbody></table>

**Table 2**

List of symbols for the phage attack model. State variables were initially set to experimental design initial conditions. Non-phage related parameters were estimated from acidification data with no phages, except for $K_L$ and $Y_{A_H/L}$ which were fixed. Phage-related parameters were estimated from the rest of the acidification data.
<table><caption>**Table 3**</caption><thead><tr><th>Symbol</th><th>Definition</th><th>Value</th><th>Relative Standard Deviation (%)</th><th>Unit</th></tr></thead><tbody><tr><td>μ<sub>max</sub></td><td>Theoretical Maximal Growth Rate</td><td>1.88</td><td>0.42</td><td>h<sup>-1</sup></td></tr><tr><td>Y<sub>AH/L</sub></td><td>Yield of Product over Substrate</td><td>0.80</td><td>Fixed according to literature value</td><td>g.(g/100L)<sup>-1</sup></td></tr><tr><td>Y<sub>L/X</sub></td><td>Yield of Substrate</td><td>0.10</td><td>1.27</td><td>g.(g/100L)<sup>-1</sup></td></tr><tr><td>γ</td><td>Consumption over Biomass Growth</td><td>0.24</td><td>0.58</td><td>-</td></tr><tr><td>K<sub>L</sub></td><td>Shape Factor</td><td>0.50</td><td>as substrate is not limiting</td><td>g/L</td></tr><tr><td>A<sub>Hmax</sub></td><td>Half-saturation Constant</td><td>9.58</td><td>0.03</td><td>g/L</td></tr><tr><td>a<sub>0</sub></td><td>Lag Parameter</td><td>0.01</td><td>1.53</td><td>-</td></tr><tr><td>α</td><td>Phage to Bacteria Ratio</td><td>8.79 ×10<sup>5</sup></td><td>0.18</td><td>(PFU/mL).(g/100L)<sup>-1</sup></td></tr><tr><td>k<sub>A</sub></td><td>Adsorption Rate</td><td>2.73 ×10<sup>-6</sup></td><td>0.02</td><td>(PFU/mL)<sup>-1</sup>·h<sup>-1</sup></td></tr><tr><td>β</td><td>Burst Size</td><td>18.90 ×10<sup>5</sup></td><td>0.35</td><td>(PFU/mL).(g/100L)<sup>-1</sup></td></tr><tr><td>λ</td><td>Cell Lysis Rate</td><td>0.92</td><td>0.26</td><td>h<sup>-1</sup></td></tr></tbody></table>

Parameter values and standard deviations after estimation. Parameters $K_L$ and $Y_{AH/L}$ were fixed.

<table><caption>**Table 4**</caption><thead><tr><th>Parameter</th><th>Estimated Value</th><th>Relative Standard Deviation (%)</th><th>Unit</th></tr></thead><tbody><tr><td>A</td><td>0.01</td><td>1.22</td><td>(mmol/L)<sup>-1</sup></td></tr><tr><td>b<sub>1</sub></td><td>8.12</td><td>1.34</td><td>mmol/L</td></tr><tr><td>c<sub>1</sub></td><td>1.07</td><td>7.08</td><td>mmol/L</td></tr><tr><td>b<sub>2</sub></td><td>1.59</td><td>1.93</td><td>unitless</td></tr><tr><td>c<sub>2</sub></td><td>0.29</td><td>6.08</td><td>unitless</td></tr></tbody></table>

Estimation of parameters for Lactic Acid to pH conversion