# Dynamic modeling of milk acidification: an empirical approach

Carlos Eduardo Robles-Rodríguez, Ewa Szymańska, Thom Huppertz, Leyla Özkan

► To cite this version:

Carlos Eduardo Robles-Rodríguez, Ewa Szymańska, Thom Huppertz, Leyla Özkan. Dynamic modeling of milk acidification: an empirical approach. Food and Bioproducts Processing, 2021, 128, pp.41-51. 10.1016/j.fbp.2021.04.010. hal-03328584

HAL Id: hal-03328584

https://hal.inrae.fr/hal-03328584

Submitted on 30 Aug 2021

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>

Distributed under a Creative Commons Attribution 4.0 International License
# Dynamic modeling of milk acidification: an empirical approach

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

Carlos Eduardo Robles-Rodríguez<sup>a,*</sup>, Ewa Szymańska<sup>b</sup>, Thom Huppertz<sup>b,c</sup>, Leyla Özkan<sup>a</sup>

## ARTICLE INFO

### Article history:

Received 29 October 2020

Received in revised form 24 March 2021

Accepted 18 April 2021

Available online 22 April 2021

### Keywords:

Acid casein

Empirical approach

Modeling

Milk acidification

pH

Yoghurt

## ABSTRACT

Acidification is an important process in the production of several dairy products and has an impact in taste, viscosity, and shelf life of the products. This paper deals with the development of an empirical mathematical model to describe milk acidification in terms of pH dynamics. The model is built upon experiments of pasteurized skimmed milk acidification for acid-induced casein precipitation and extended to predict the partitioning of calcium and phosphates in the whey and the dynamics of pH during fermentation of yoghurt manufacture. Furthermore, the model can be used directly or in optimization problems to provide recommendations about product design.

© 2021 The Author(s). Published by Elsevier B.V. on behalf of Institution of Chemical Engineers. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).

## 1. Introduction

Acidified milk products are some of the oldest and most popular foodstuffs. The popularity of fermented milks and yoghurts is, at least partly, due to various health claims and therapeutic benefits that have been attributed to consumption of these products (de Oliveira, 2014). One important process in the production of these fermented dairy products is acidification.

Acidification is an important operation in food industry, as it affects the taste, viscosity, and shelf life of dairy products and non-dairy products (Mudgil and Barak, 2019). In fermented milk manufacture, the pH of the product typically decreases to values in the pH range 4.0-4.8, compared to a natural milk pH of 6.6-6.8. A major application of acidification in dairy product manufacture is in the production of acid casein, the base material from which caseinates can be prepared. Caseins and

caseinates are used as ingredients in a wide variety of food and nonfood products (Sarode et al., 2015). Major differences between the acidification of milk for acid casein manufacture, acidified dairy products and the fermented dairy products are in the type of acid that is used and the time-scale of acidification. For acid casein manufacture, mineral acids are added rapidly for near-instantaneous acidification; whereas for fermented milk products, acidification occurs slowly, over a period of hours, due to the production of lactic acid by lactic acid bacteria.

When the pH of milk is decreased, the different acido-basic groups of milk's constituents (organic and inorganic phosphate, citrate, carboxylic acid residues) become increasingly more protonated, which has two major effects. The first is that charge neutralization of the κ-casein segments protruding from the casein micelle surface leads to loss of solvency and therewith reduced steric and electrostatic stabilization and concomitant aggregation of casein micelles (Lucey, 2017). Furthermore, the colloidal calcium phosphate (CCP), which was originally in the casein micelles, is dissolved and this

<footnote><sup>a</sup> Control Systems Group, Department of Electrical Engineering, Eindhoven University of Technology, De Zaale, 5612</footnote>

<footnote>AJ Eindhoven, The Netherlands</footnote>

<footnote><sup>b</sup> FrieslandCampina, Amersfoort, The Netherlands</footnote>

<footnote><sup>c</sup> Wageningen University and Research, Wageningen, The Netherlands</footnote>

<footnote>* Corresponding author. Present address: TBI, Université de Toulouse, CNRS, INRAE, INSA, Toulouse, France.</footnote>

<footnote>E-mail address: roblesro@insa-toulouse.fr (C.E. Robles-Rodríguez). https://doi.org/10.1016/j.fbp.2021.04.010</footnote>
alters the micelle's internal structure. The extent of CCP solubilization depends on both pH and temperature (Holt et al., 1981).

During fermentation of milk with lactic acid bacteria, these processes occur slowly, over a time scale of several hours, and a large proportion of the CCP has solubilized at the point when coagulation of casein micelles commences, e.g., from pH ~ 5.0 in the case of unheated milk. During acid casein manufacture, however, the very rapid acidification can lead to several processes occurring in very short succession or even simultaneously, i.e., mixing of acid and milk; acid-induced coagulation, solubilization of CCP, aggregate growth and aggregate break-up. It has been demonstrated that mixing of acid and milk is a rate-controlling step in the formation and growth of the precipitate particles (Jablonka et al., 1988).

In order to gain understanding, mathematical models have been proposed for salt specification in milk and the influence of acidification thereof. For instance, the works of Mekmene et al. (2009) and Holt (2004) have focused on the detailed description of salt equilibria in milk. These models are accurate in steady state calculation and they require a lot of calculations about many ionic species. Hofland et al. (2003) studied the dynamics of milk acidification with a model to describe the pH profiles considering the electro-neutrality condition, where information about the particle size and the coagulation time are needed for the calculation of the model. Although the model is accurate for the analysis of the pH data, it is data-dependent because it needs to specify the coagulation time. Other works have focused on surface methodologies and statistical models to relate temperature and rheological parameters (Kristo et al., 2003).

The aim of this work is two-fold: to develop a dynamic model to describe the pH dynamics during acidification of milk and to test several uses of the model for experimental and product design. To this end, we propose an empirical model to describe milk acidification in a simple manner capturing the effect of measured variables (i.e. temperature, acidification times, and incubation times). The empirical model is built upon experimental data for the acidification of pasteurized skimmed milk. Additionally, the model is extended to predict calcium, phosphates and protein distribution during milk acidification. A second extension of the model is introduced and applied to literature data of yoghurt acidification where it is shown that the model structure is robust, but a new parameter estimation is required to adapt the model to the different data set.

# 2. Materials and methods

## 2.1. Experimental set-up

The precipitation experiments were carried out in a 200 mL jacketed glass beaker. The pH was monitored using a Mettler Toledo DL50 Graphix autotitrator (Mettler Toledo, Switzerland) operating in equivalence/end point titration method. The titrator was equipped with a DL7X rotor (Mettler Toledo, Switzerland) with a three bladed propeller stirring rod (Mettler Toledo, Switzerland), a DG115-SC pH electrode (Mettler Toledo, Switzerland) and a DT120 Pt100 temperature sensor (Mettler Toledo, Switzerland).

The titration was controlled by the LabX titration V2.6 software (Mettler Toledo, Switzerland) where the flow of acid was fixed to be constant. The beaker temperature was con-

<table><caption>Table 1 – Conditions of the experiments of pasteurized skimmed milk with a protein content of 3.47% w/w performed at a constant stirring speed of 500 rpm.</caption><thead><tr><td>Experiment</td><td>t<sub>A</sub> (min)</td><td>T (°C)</td><td>pH<sub>0</sub></td></tr></thead><tbody><tr><td>1</td><td>0</td><td>44</td><td>6.42</td></tr><tr><td>2</td><td>2.5</td><td>44</td><td>6.45</td></tr><tr><td>3</td><td>10</td><td>44</td><td>6.46</td></tr><tr><td>4</td><td>30</td><td>44</td><td>6.44</td></tr><tr><td>5</td><td>0</td><td>20</td><td>6.62</td></tr><tr><td>6</td><td>2.5</td><td>20</td><td>6.67</td></tr><tr><td>7</td><td>10</td><td>20</td><td>6.66</td></tr><tr><td>8</td><td>30</td><td>20</td><td>6.63</td></tr><tr><td>9</td><td>0</td><td>5</td><td>6.74</td></tr><tr><td>10</td><td>2.5</td><td>5</td><td>6.80</td></tr><tr><td>11</td><td>10</td><td>5</td><td>6.78</td></tr><tr><td>12</td><td>30</td><td>5</td><td>6.76</td></tr></tbody></table>

trolled using a Julambo F33 water bath (Julambo Labortechnik, Seelbach, Germany) operating at medium flow. The starting temperature was determined with an IKA ETS-D4 Fuzzy temperature controller connected to an IKA-RH KT/C stirring and heating plate (IKA, Germany).

## 2.2. Experiments of milk acidification

A series of experiments was carried out consisting of the acidification of pasteurized skimmed milk at different acidification rates and temperatures T. The experiments investigate the dynamics of the precipitation of casein in pasteurized skimmed milk regarding pH measurements. Before acidification 150 ± 2 mL of pasteurized skimmed milk was equilibrated at the required temperature (5, 20 and 44 °C) in the jacketed glass beaker. Subsequently, the samples were acidified with a solution of 0.5 M sulfuric acid. A total volume of acid (A<sub>T</sub>) of 7.92 mL of H<sub>2</sub>SO<sub>4</sub> was found to be required to lower the pH of the 150 mL of the pasteurized skimmed milk to 4.6. The acid was added over an acidification time (t<sub>A</sub>) of 0, 2.5, 10, or 30 min. The acidification time t<sub>A</sub> = 0 means that all acid has been added very fast (a short pulse addition). The acidification rate was defined as Q<sub>A</sub> = A<sub>T</sub>/t<sub>A</sub>. All the experiments were performed at the same stirring rate of 500 rpm. The pH measurements had an accuracy of ±0.02 and were recorded every 10 s.

The experimental data representing the titration curves of the milk are displayed in Fig. 1, where the replicates of the curves are also presented. The mean values of the data sets were calculated in order to use only one set per acidification time and temperature. In total, 12 data sets were gathered whose characteristics and initial conditions are reported in Table 1.

## 2.3. Modeling pH dynamics in milk acidification

An empirical model to capture the dynamics of pH is developed based on the measured variables, such as acidification times t<sub>A</sub> (min), temperature T (K), the acid addition rate Q<sub>A</sub> (L<sub>acid</sub>/min) and the volume of milk V<sub>milk</sub> (L<sub>milk</sub>). For the sake of simplicity, let us define A as,

$$A = \frac{\text{Cumulativevolumeofaddedacid}}{\text{Volumeofmilk}} = \sum \frac{Q_A}{V_{\text{milk}}}$$

which represents the cumulative added acid per volume of milk (L<sub>acid</sub>/L<sub>milk</sub>) and can be considered dimensionless.
<figure><img src="image_3.png" /><figcaption>Fig. 1 – Experimental data for milk acidification at different temperature and acidification times. The numbers in the parenthesis indicate the replicate number.</figcaption></figure>

The modeling approach consists in finding a function φ for pH depending on the measured variables and mixing in a structured manner, where independent functions f(·) interact. Therefore, the model can be expressed as:

$$pH = \phi(A, t_A, T) = f(A) \cdot f(T) \cdot f(t_A) \cdot f(\text{mixing}) \quad (1)$$

From the data displayed in Fig. 1, two stages can be distinguished: (i) acid addition stage where the acid is being added to the milk ($t \le t_A$); and (ii) diffusion stage ($t > t_A$), where the acid that has not been used solubilizes in the milk until pH settles around an equilibrium pH value thanks to the action of mixing.

In this work, we make use of an interesting property observed in the data. For instance, when the dynamics of pH are displayed against A, a time invariant relationship is found. This invariant relationship is displayed in Fig. 2. It is observed that the variation of pH in terms of A is nonlinear and it follows a polynomial function. It is important to note that this invariant behavior holds for large acidification times and lower temperatures (i.e., 5 and 20 °C). There is, however, an inflection point at pH values around 5.1. According to Hofland et al. (2003), around this pH value small particles would rapidly grow and the buffering groups of the protein and residual micellar calcium phosphate become less accessible for the added acid. This implies that at pH values around 5.1 caseins have precipitated and the micellar calcium phosphate is completely solubilized (Mekmene et al., 2010). At this moment, the trajectories of pH (Fig. 2 can be considered to vary according to a quadratic relationship as a function of the cumulative added acid per volume of milk A. This empirical relationship is used to build $f_1(A, T)$, which incorporates the variation of temperature as an exponential function expressed similarly to an Arrhenius type of equation. Temperature is included since it has a big effect on milk acidification, especially on casein aggregation and milk coagulation. For instance, at low temperature (<10 °C) coagulation does not occur (Lucey and

Singh, 1997). This is also observed in Fig. 2, where there is not an evident change at pH 5.1 at lower temperatures (5 and 20 °C) indicating that no coagulation is taking place and only small particles of casein are formed.

Below pH 5.1, the slope of the curves in Fig. 2(a) is steeper. This is probably related to the fact that almost all caseins have precipitated, the aggregation proceeds to form bigger particles, and the acid is being diffused in the serum. This diffusion acts similarly to a buffer (Hofland et al., 2003). In this case, we propose to represent the pH dynamics below 5.1 by another quadratic function that also depends on temperature $f_2(A, T)$. However, we consider that the pH variations also depend on mixing and the acidification times $t_A$. This is observed in Fig. 1a where pH decreases abruptly (undershoots) when acidification times of $t_A = 0, 2.5$ min are used at 44 °C. For the sake of clarification, let us define the pH undershoot as the difference between the minimum pH and the equilibrium pH value at the end of the experiment. The pH undershoot is neither visible when the acidification occurred over $t_A = 30$ min nor when lower temperatures ($T = 5, 20$ °C) are explored. Therefore, the parameters related to temperature must have a different magnitude and sign to indicate that at lower temperatures the effect of acid addition below pH 5.1 is smaller with temperature. Additionally, two corrections terms denoted by $\beta$ are included in the model to take into account for the changes related to the different acidification times, especially for the shorter times where the reaction/precipitation is faster than mixing. Furthermore, we consider a term for mixing depending on the mixing time which is defined as,

$$t_m = \frac{V_T}{r_c d_s^3 N_i}$$

where $V_T$ is the volume of the beaker (m³), $d_s$ is the diameter of the impeller (m), $N_i$ is the stirring speed (rpm) and $r_c$ is a constant to relate the diameter of the precipitate and micelles (Hofland et al., 2003). In this case, $r_c$ will be fixed to 1 and considered to be absorbed by other constants in the model.
<figure><img src="image_4.png" /><figcaption>Fig. 2 – pH as a function of added acid per volume of milk (A) at different temperatures and acidification times.</figcaption></figure>

The model for the two parts of the acid addition stage can be expressed as:

$$ \text{pH(A, T, t_A) - pH}_0 = \begin{cases} f_1(A, T), & \text{if pH} \ge 5.1 \\ f_2(A, T, t_A), & \text{if pH} < 5.1 \end{cases} \qquad (2) $$

where pH(A) denotes a function of pH depending on the acid addition and pH₀ is the initial pH. The empirical expressions of f₁(A, T) and f₂(A, T, t_A) correspond to

$$ f_1(A, T) = (\theta_1 A^2 + \theta_2 A)e^{\left(\frac{\lambda_1}{T}\right)} \qquad (3a) $$

$$ f_2(A, T, t_A) = t_m \left(\frac{1}{t_A}\right)^{\beta_1} (\theta_3 A^2 + \theta_4 A) e^{\left(\frac{\lambda_2}{T}\right) \left(\frac{1}{t_A}\right)^{\beta_2}} \qquad (3b) $$

The parameters θ₁–₄ are constants to deal with the quadratic representations of the change of pH with respect to A with units expressed as pH value/min and pH value/min² for θ₁,₂ and θ₃,₄ respectively. The terms λ₁–₂ are the temperature related constants (K). The parameters β₁ and β₂ are dimensionless constants that allow to generalize the model for different acidification rates and temperatures. In order to consider t_A = 0, we have assumed that the acidification time is shorter than the response of the pH probe, which is lower than 9 s (Hofland et al., 2003). Hence, for modeling purposes, we have set t_A = 3 s (0.05 min) when referring to t_A = 0 min.

The second stage of the experiments considers a pure diffusion process involving the increase of pH to values close to the isoelectric point where the pH values achieve equilibrium. This implies that most of the reactive groups are accessed within the time frame of the experiments. From the experiments (Fig. 1), it is observed that the pH decreases and then stabilizes due to the action of mixing and buffering. In order to find a function that describes the changes of pH, we have

<figure><img src="image_5.png" /><figcaption>Fig. 3 – pH variations after acid addition.</figcaption></figure>

plotted $\frac{dpH}{dt}$ with respect to time in Fig. 3. It can be observed that the trend follows an exponential decay that depends on t_A. Furthermore, the same type of function for temperature is also added to keep consistency with the previous stage of the model. The model for the second stage is then expressed as:

$$ \frac{dpH}{dt} = f_{\text{diff}} = \theta_5 t_m \left(\frac{1}{t_A}\right)^{\beta_3} e^{(\theta_6 (t-t_A))} e^{\left(\frac{\lambda_2}{T}\right)} \qquad (4) $$

where θ₅–₆ are constants with dimensions pH value/min² and 1/min, respectively. λ₃ is a constant to deal with temperature changes in (K). The parameter β₃ is introduced to manage the differences between acidification times.

From Eqs. (2) and (4), it is possible to obtain the dynamics of pH and H⁺ (H⁺ = 10⁻pH) by applying the chain rule of
differentiation. Hence, the complete model can be expressed as:

$$ \frac{dH^+}{dt} = \begin{cases} -\frac{\partial pH}{\partial A} \frac{\partial A}{\partial t} H^+ \ln(10), & \text{if } Q > 0 \\ -f_{\text{diff}} H^+ \ln(10), & \text{if } Q = 0 \end{cases} \qquad (5) $$

$$ \frac{dA}{dt} = Q \qquad (6) $$

where Q is the constant acid addition rate ($L_{acid}/L_{milk}/min$), whose cumulative sum is A. The model comprises two differential equations and 12 parameters to describe the change of pH at different temperatures and acidification times.

It is important to note that this modeling approach opens the possibility to propose different types of functions $f(\cdot)$ for each variable defined in Eq. (1).

# 3. Results

## 3.1. Calibration and validation of the model for milk acidification

The model has been calibrated using the experimental data described in Section 2.2, where two stages have been characterized.

The mean values of the data presented in Fig. 1 have been used to identify the model parameters regardless of the differences in the replicates at lower $t_A$. Parameter identification has been performed dividing the data into acid addition stage above pH 5.1, acid addition stage below pH 5.1, and diffusion stage. First, the parameters $\theta_1$, $\theta_2$ and $\lambda_1$ have been identified with nine data sets comprising experiments 2–4, 6–8 and 10–12 as defined in Table 1 with the values of pH $\ge$ 5.1. The data set for $t_A = 0$ is not used here since the second value of pH was already around pH of 4.6, making difficult to infer the dynamic behavior between the initial pH and this second point. A similar approach is used for the determination of the parameters $\theta_3$, $\theta_4$ and $\lambda_2$, which, in this case, are identified per $t_A$. This implies that three sets of parameters are available, whose differences are attributed to $t_A$ and mixing. However, as the mixing was constant along the experiments, only changes in $t_A$ were considered. In order to find a unique set of parameters for all $t_A$ and T in the model, we have proposed to add the correction factors $\beta_1$ and $\beta_2$ to take into account the differences regarding $t_A$. The values for $\beta_{1-2}$ are first inferred from the variations of the parameters $\theta_3$, $\theta_4$ and $\lambda_2$. Then, a re-calibration has been made to find an optimal value for the parameters involved in $f_2(A, T, t_A)$. It is worth noting that estimating the five parameters at the same time may hamper identifiability. This was addressed by setting small bounds for the parameters $\beta_{1-2}$. Finally, the twelve data sets for values after acid addition were used to identify the parameters of Eq. (4). Parameter identification have been performed by a global gradient free method called pattern search in MATLAB searching to minimize the sum of the squared errors of the experimental data with respect to the model. Sensitivity analysis has been performed to assess the influence of each parameter on the measured output (pH) as well as the possible interaction (correlation) between these effects. The local sensitivity was studied according to the first-order sensitivities for Eq. (5) as:

$$ \dot{S} = \frac{\partial f(H^+)}{\partial H^+} S + \frac{\partial f(H^+)}{\partial P} \qquad (7) $$

<table><caption>Table 2 – Parameters of the dynamic model.</caption><thead><tr><th></th><th>Value</th><th>Value</th></tr></thead><tbody><tr><td>θ<sub>1</sub> [pH/min]</td><td>157.23 ± 24.68</td><td>439.31 ± 43.34</td></tr><tr><td>θ<sub>2</sub> [pH/min]</td><td>-15.64 ± 2.18</td><td>-485.42 ± 70.54</td></tr><tr><td>θ<sub>3</sub> [pH/min<sup>2</sup>]</td><td>3.01e4 ± 4.05</td><td>-3511.34 ± 159.97</td></tr><tr><td>θ<sub>4</sub> [pH/min<sup>2</sup>]</td><td>-6.1e3 ± 1.93</td><td>0.273 ± 55.44</td></tr><tr><td>θ<sub>5</sub> [pH/min<sup>2</sup>]</td><td>1200 ± 0.03</td><td>0.15 ± 94.42</td></tr><tr><td>θ<sub>6</sub> [1/min]</td><td>-0.0066 ± 0.015</td><td>0.125 ± 4.47</td></tr></tbody></table>

where $f(H^+)$ denotes the right hand side of Eq. (5), P the 12 parameters of the model, and S the sensitivity matrix. As pH, or in this case $H^+$ is the only output, the sensitivities for $H^+$ with respect to the parameters P were calculated as

$$ \frac{\partial H^+}{\partial P} = S \qquad (8) $$

The sensitivity functions are evaluated at the identified parameter values P and used to compute the Fisher information matrix (FIM) which provides information about the confidence intervals and the correlation of the parameters.

The identified parameters together with their confidence intervals are reported in Table 2. The correlation matrix of the parameters is also shown in Table 3. The correlations were calculated only for the parameters involved in each partition of the model: 3 parameters for part one, and 5 parameters for part two of the acidification stage; and 4 parameters corresponding to the diffusion stage. Therefore there is no correlation between the parameters that correspond to other stage. It can be seen in Table 3 that some parameters are correlated, especially for the terms of the polynomials and the parameters for temperature $\lambda_{1-3}$, because they are multiplying each other. Furthermore, the parameters with larger confidence interval indicate that other parameters must be fixed. This is the case, for instance, of $\beta_1$ and $\beta_2$ whose values are determining the values of $\theta_3$, $\theta_4$ and $\lambda_2$.

The comparison of the model predictions with the experimental data is presented in Fig. 4 as solid lines for the different experimental conditions. It can be observed that when acidifying for a longer time, the pH gradually decreases resulting in slower precipitation. When longer acidification times ($t_A$) are applied, the curd spends longer time at higher pH values where precipitation is slow, resulting in finer particles (Hofland et al., 2003). Rapid acidification of milk $t_A = 0$ results in an undershoot, which can be explained by the rapid casein micelle coagulation forming large aggregates and a less gradual pH decrease (Bringe and Kinsella, 1990). The pH undershoot is not visible at low temperatures (5°C and 20°C). For instance, the whey clearly separates from the curd at 44°C, whereas at 20°C less separation occurs and at 5°C there is no separation.

The pH dynamics are followed quite accurately by the model simulation for all the data points. The performance of the model is also evaluated by means of the root mean squared errors (RMSE), whose values are reported in Table 4. The calibration errors are around 0.07 in pH units as an average of the experiments. However, we observe that the model predictions present lower errors for $t_A = 30$ min due to the presence of both parts of the acid addition stage and the diffusion stage in these data sets. It is also important to note that the model holds for all the data sets, and thus there is a trade-off between the fitting accuracy of all the data sets. The model was developed for a constant stirring speed, however, it contains a term for mixing. The model has been evaluated against an independent
<table><caption>Table 3 – Correlation matrix of the parameters.</caption><thead><tr><th></th><th>θ<sub>1</sub></th><th>θ<sub>2</sub></th><th>λ<sub>1</sub></th><th>θ<sub>3</sub></th><th>θ<sub>4</sub></th><th>λ<sub>2</sub></th><th>β<sub>1</sub></th><th>β<sub>2</sub></th><th>θ<sub>5</sub></th><th>θ<sub>6</sub></th><th>λ<sub>3</sub></th><th>β<sub>3</sub></th></tr></thead><tbody><tr><th>θ<sub>1</sub></th><td>1</td><td>-0.158</td><td>0.501</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>θ<sub>2</sub></th><td></td><td>1</td><td>-0.926</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>λ<sub>1</sub></th><td></td><td></td><td>1</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><th>θ<sub>3</sub></th><td></td><td></td><td></td><td>1</td><td>0.925</td><td>-0.035</td><td>-0.037</td><td>0.020</td><td></td><td></td><td></td><td></td></tr><tr><th>θ<sub>4</sub></th><td></td><td></td><td></td><td></td><td>1</td><td>-0.408</td><td>0.105</td><td>0.237</td><td></td><td></td><td></td><td></td></tr><tr><th>λ<sub>2</sub></th><td></td><td></td><td></td><td></td><td></td><td>1</td><td>-0.263</td><td>-0.494</td><td></td><td></td><td></td><td></td></tr><tr><th>β<sub>1</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>0.961</td><td></td><td></td><td></td><td></td></tr><tr><th>β<sub>2</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td></td><td></td><td></td><td></td></tr><tr><th>θ<sub>5</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>-0.873</td><td>-0.909</td><td>-0.070</td></tr><tr><th>θ<sub>6</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>0.681</td><td>-0.045</td></tr><tr><th>λ<sub>3</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td><td>0.048</td></tr><tr><th>β<sub>3</sub></th><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td>1</td></tr></tbody></table>

<figure><img src="image_6.png" /><figcaption>Fig. 4 – Model performance against experimental data. Dots represent experimental data, solid lines are model outputs for the different acidification times t<sub>A</sub>.</figcaption></figure>

<table><caption>Table 4 – Performance evaluation of the model in the calibration and validation data sets.</caption><thead><tr><th>t<sub>A</sub> (min)</th><th>RMSE</th><th>Std dev.</th><th>RMSE validation</th><th>Std dev.</th></tr></thead><tbody><tr><td>0</td><td>0.082</td><td>0.057</td><td>–</td><td>–</td></tr><tr><td>1</td><td>–</td><td>–</td><td>0.091</td><td>0.072</td></tr><tr><td>2.5</td><td>0.037</td><td>0.027</td><td>–</td><td>–</td></tr><tr><td>6</td><td>–</td><td>–</td><td>0.092</td><td>0.063</td></tr><tr><td>10</td><td>0.040</td><td>0.022</td><td>0.093</td><td>0.086</td></tr><tr><td>16</td><td>–</td><td>–</td><td>0.040</td><td>0.032</td></tr><tr><td>30</td><td>0.049</td><td>0.034</td><td>–</td><td>–</td></tr></tbody></table>

dent data set performed at higher mixing conditions, whose results are displayed in Fig. 5. These results highlight that the model can be used at different mixing conditions. The initial conditions and RMSE values for this set are reported in Table 5.

In order to assess the validity of the model, we have compared the model against literature data (Hofland et al., 2003)

reported for similar stirring speeds and a temperature of 40 °C. The results of the model with this experimental data set are shown in Fig. 6. It is observed that the model follows quite accurately the dynamic behavior of the pH. However, the change in the mixing conditions (i.e., volume and impeller diameter) make the model to predict higher values for short
<figure><img src="image_7.png" /><figcaption>Fig. 5 – Model performance against experimental data carried out at 44 °C and stirring speed of 600. Dots represent experimental data, solid lines are model outputs.</figcaption></figure>

<table><caption>Table 5 – Initial pH and performance evaluation of the model with the set with higher mixing (800 rpm).</caption><thead><tr><td>t<sub>A</sub> (min)</td><td>RMSE calibration</td><td>Std dev.</td><td>pH<sub>0</sub></td></tr></thead><tbody><tr><td>0</td><td>0.12</td><td>0.012</td><td>6.39</td></tr><tr><td>2.5</td><td>0.022</td><td>0.016</td><td>6.40</td></tr><tr><td>5</td><td>0.048</td><td>0.021</td><td>6.39</td></tr><tr><td>10</td><td>0.028</td><td>0.015</td><td>6.42</td></tr><tr><td>30</td><td>0.028</td><td>0.026</td><td>6.39</td></tr></tbody></table>

<figure><img src="image_8.png" /><figcaption>Fig. 6 – Validation of the model. Data from Hofland et al. (2003). Stirring rate 500 rpm, temperature 40 °C.</figcaption></figure>

acidification times $t_A$. Additionally, the model is not able to predict a smoother trajectory for $t_A = 10$ min. This can be explained by the fact that Hofland et al. (2003) used a different type of impeller which is not considered in the model. Furthermore, the experiments for validation were carried out in a 1 L reactor for the experiments, whereas the lab experiments for calibration presented in Section 2.2 were performed in flasks.

Similar to the experiments, the validation of the model has also been assessed by the RMSE (Table 4). The RMSE values for the different data sets are low leading to the conclusion that the model can be used to predict the dynamics of pH at different $t_A$ and temperatures. However, it is the interest of this

<figure><img src="image_9.png" /><figcaption>Fig. 7 – Model prediction for protein yield against data of protein yield evolution during acidification with respect to pH (Hofland et al., 1999).</figcaption></figure>

paper to use this model to gain insights about the dynamics of other components and to use it for product design.

## 3.2. Modeling the partitioning of protein and salts

A way to measure the impact of the acidification is the measurement of the calcium and phosphates in the whey. In this work, we propose to extend the model presented in Section 2.3 in order to track the partitioning of calcium and phosphates. To this end, we have taken the data from the work of Hofland et al. (1999) who have reported the changes of protein precipitation yield, the calcium concentration in the whey and the phosphates concentration in the whey with respect to pH when skimmed milk was acidified with the addition of 0.5 M of $H_2SO_4$. It is worth noting that the reported values from Hofland et al. (1999) have been obtained at similar conditions of the experiments presented in this paper. As the reported values were expressed with respect to pH, it is possible to make a connection or at least to infer the dynamics of the reported data. In the next subsections, we present the development of the proposed empirical models for the dynamics of the three minerals mentioned above. The choice of the model structure has been performed following the trend of the figures, which present a sigmoidal behavior. Therefore, we propose the use of logistic equations.

### 3.2.1. Partitioning of the protein

Hofland et al. (1999) reported steady state data for the yield of protein precipitation as shown in Fig. 7. The protein yield was calculated from the residual protein concentration in the whey assuming that the whey proteins did not co-precipitate and the final value was considered as a complete precipitation. Therefore, the final values are always 1. Nonetheless, the aim in this paper is to propose a model to describe the dynamics of the proteins (Pr) depending on pH and temperature. To this end, let us consider a regression of the data following a logistic function as:

$$Y_{Pr} = \frac{\gamma_1 e^{(\lambda_{Pr}/T)}}{e^{(\lambda_{Pr}/T)} + e^{(-\gamma_2 H^+ + \gamma_3)}} \quad (9)$$

where $Y_{Pr}$ is the protein yield, $T$ represents the actual temperature, $\lambda_{Pr}$ holds for a constant to deal with temperature changes (K) and $H^+$ are the hydrogen protons. The terms $\gamma_1-3$
<figure><img src="image_10.png" /><figcaption>Fig. 8 – Model prediction for calcium concentration data from Hofland et al. (1999).</figcaption></figure>

are constant parameters, where $γ_1$ is the maximum attained value which in this case is equal to 1.

The results of the proposed regression in Eq. (9) are displayed in Fig. 7, where a good agreement between the model output and the experimental data can be appreciated. The logistic function considers that $γ_1$ is the maximum value that can be achieved. In this case, it corresponds to 1 and thus the dynamics will converge toward this value.

The accuracy of the model is confirmed by the determination coefficient (R²) of 0.972 with standard deviation of error of 0.062. It is worth noting that as the values were obtained from 1 backwards, the predictions of the model do not start at zero, and therefore the model has to extrapolate to obtain those values.

### 3.2.2. Partitioning of calcium and phosphates in the whey

Similar to the protein content, data for the concentrations of calcium and phosphates in the whey has been obtained from Hofland et al. (1999). In this case, the same type of logistic model is proposed to also fit the titers of calcium in the whey:

$$Ca_{whey} = \frac{\alpha_1 e^{(\lambda_{Ca}/T)}}{1 + e^{(-\alpha_2 H^+ + \alpha_3)}} \quad (10)$$

where $Ca_{whey}$ stands for the calcium concentration in the whey and $α_{1–3}$ are the parameters of the model. $λ_{Ca}$ is the term to deal with temperature variations. The structure is different to the presented in Eq. (9), where the exponential term in the denominator was erased to consider variations in the maximum concentration of calcium. The values of the parameters $α$ are reported in Table 6. The fit of the model against experimental data is shown in Fig. 8, where the determination coefficient (R²) has been found to be 0.987 and standard deviation of 41.11.

Regarding phosphates, a different structure is proposed due to the fact that the steady state value changes with temperature. The mathematical representation for the $Phos_{whey}$ is as follows:

$$Phos_{whey} = \frac{\delta_1 e^{(\lambda_{Phos}/T)}}{1 + e^{(-\delta_2 H^+ + \delta_3)}} \quad (11)$$

where $δ_{1–3}$ are the parameters of the model, with $λ_{Phos}$ being a constant to deal with temperature. The parameters are reported in Table 6. It is important to mention that the phosphates have been studied in terms of phosphorus con-

<figure><img src="image_11.png" /><figcaption>Fig. 9 – Model prediction of phosphorus concentrations in the casein phosphate data from Hofland et al. (1999).</figcaption></figure>

centrations, as in Hofland et al. (1999). The performance of the model against the experimental data is presented in Fig. 9, where a good agreement between the data and the model is observed. A value of 0.989 for the determination coefficient (R²) and a value of 21.61 of standard deviation have been found.

The models for calcium and phosphates share the same structure and their parameters are in similar ranges, which highlights their intrinsic connection. Furthermore, the exponential terms of the logistic function confirm the fact that pH can be used to approximate their dynamic behavior. Nonetheless, the results are only for the calibration procedure, and more data is needed for validation.

### 3.2.3. Dynamics of the partition components in the precipitation

The idea of using the relationships of protein yield, calcium and phosphate with respect to pH is to exploit their dependency on pH in order to approximate their dynamics. This can be done by the use of the chain rule with the pH dependent expressions. Hence, the dynamics of the protein yield $Y_{Pr}$, calcium $Ca_{whey}$ and phosphates $Phos_{whey}$ can be approximated as:

$$\frac{dY_{Pr}}{dt} = \frac{\partial Y_{Pr}}{\partial pH} \frac{\partial pH}{\partial t} - \frac{Y_{Pr}}{V} \frac{dV}{dt} \quad (12)$$

$$\frac{dCa_{whey}}{dt} = \frac{\partial Ca}{\partial pH} \frac{\partial pH}{\partial t} - \frac{Ca}{V} \frac{dV}{dt} \quad (13)$$

$$\frac{dPhos_{whey}}{dt} = \frac{\partial Phos}{\partial pH} \frac{\partial pH}{\partial t} - \frac{Phos}{V} \frac{dV}{dt} \quad (14)$$

where V is the volume of the tank/reactor and the term containing the derivative $dV/dt$ is added to account for the dilution. Considering that most of the experiments were performed in fed-batch, the dilution term $dV/dt$ appears. However, this could be adjusted depending on the operation regime.

The model presented described by Eqs. (5) and (6) has been coupled to Eqs. (12)–(14) to compute the dynamic behavior of the protein, calcium and phosphates, whose results are displayed in Fig. 10 for the experiments at 44 °C considering different acidification rates. The initial conditions for the model were 0.01, 400 ppm and 400 ppm for protein yield, calcium and phosphorus, respectively.

In Fig. 10, it is observed that the precipitation always arrives at steady state before the $t_A$ which is in line with the
<table><caption>Table 6 – Parameters of the models for protein yield, calcium and phosphates.</caption><thead><tr><th>Protein</th><th>Value</th><th>Calcium</th><th>Value</th><th>Phosphate</th><th>Value</th></tr></thead><tbody><tr><td>γ<sub>1</sub></td><td>1</td><td>α<sub>1</sub></td><td>435.5 ± 366.5</td><td>δ<sub>1</sub></td><td>201.38 ± 73.2</td></tr><tr><td>γ<sub>2</sub></td><td>6.12e5 ± 1.02e5</td><td>α<sub>2</sub></td><td>5.1e5 ± 7.69e4</td><td>δ<sub>2</sub></td><td>5.13e5 ± 6.5e4</td></tr><tr><td>γ<sub>3</sub></td><td>-32.4 ± 8.4</td><td>α<sub>3</sub></td><td>0.73 ± 0.14</td><td>δ<sub>3</sub></td><td>-0.29 ± 0.07</td></tr><tr><td>λ<sub>Pr</sub> [K]</td><td>-1.12e4 ± 2.7e3</td><td>λ<sub>Ca</sub> [K]</td><td>355.1 ± 256</td><td>λ<sub>Phos</sub> [K]</td><td>385.3 ± 110.9</td></tr></tbody></table>

<figure><img src="image_12.png" /><figcaption>(a) Evolution of the precipitation of protein yield.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>(b) Calcium concentration in the whey.</figcaption></figure>

<figure><img src="image_14.png" /><figcaption>(c) Phosphorus concentration in the whey.</figcaption></figure>

<caption>Fig. 10 – Prediction of the dynamics of the partition components in the precipitation at 44 °C using different acidification times.</caption>

assumption that almost all caseins have precipitated at pH 5.1 (Mekmene et al., 2010; Hofland et al., 2003). These times indeed correspond to the times where pH 5.1 is reached.

In general, the steady state values are achieved faster than the stabilization of pH, which indicates that the diffusion is very important in the pH, whereas it does not have a strong impact in the release of components. The proposed empirical model can be used with different initial conditions. Nonetheless, we should mention that this extension is reported only for calibration results, which means that further validation and possibly re-estimation of parameters is required.

### 3.3. Yoghurt acidification

Yoghurt is widely consumed as a functional food due to its good taste and nutritional properties (rich in potassium, calcium, protein and vitamin B) and excellent vehicle to deliver probiotics to consumers (de Oliveira, 2014). Yoghurt is manufactured from milk fermentation with the use of thermophilic lactic acid bacteria, such as *Streptococcus thermophilus* and *Lactobacillus delbrueckii* sp. bulgaricus (Kristo et al., 2003). Fermentation of lactose to lactic acid by these bacteria gradually causes the pH of the milk to drop from about 6.7 to values around 4–5. This pH decrease results in a gel structure as a consequence of the coagulation of the milk proteins, due to the lactic acid secreted by defined species of bacteria cultures. Additionally, yoghurt can include some derivatives

of milk such as skimmed milk powders, whey concentrates, caseinates or creams (Sfakianakis and Tzia, 2014).

The production of yoghurt involves gradual acidification at 40–45 °C where the thermophilic bacteria produces the lactic acid which acidifies the milk and causes the precipitation of casein micelles. The phenomena occurring in acidification of yoghurt is, however, comparable to the acidification of milk. Both encounter several stages when pH drops from 6.7 to 4.6 (Lucey and Singh, 2003): (1) decrease in the net negative charge on the micelle due to the diffusion of protons; (2) neutralization of the micelles surfaces but electrostatic destabilization of micelles causing the colloidal calcium phosphates to be completely dissolved when pH 5; and (3) destabilization of micelles, gelation/coagulation occurs at pH 4.9.

Several factors might affect the acidification of yoghurt such as temperature, incubation time, addition of emulsifiers and fruit flavors, agitation, type of milk, etc. (de Oliveira, 2014), which have an impact on the texture, flavor and quality of the final yoghurt. For instance, texture of yoghurt is the result of both acid aggregation of casein micelles and production of exopolysaccharide (EPS), which is very important for the quality of the product. Hence, controlling the production of yoghurt is quite a challenge. Notwithstanding, one characteristic that can be simply monitored is the pH. When a pH value is specified, by means of the final product quality profile, then online control of the fermentation process can be carried out by monitoring pH (De Brabandere and De Baerdemaeker,
<figure><img src="image_15.png" /><figcaption>Fig. 11 – Model output for profiles of pH at 40 °C taken from Aguirre-Ezkauriatza et al. (2008). The symbols correspond to the samples taken at different locations where the shapes are as indicated in the paper. The solid line is the output of the model.</figcaption></figure>

1999). This entails several advantages such as low investment and functional cost, easy implementation, and absence of any fluctuations of yoghurt coagulation. In this context, it is our intention to extend the model proposed for milk acidification to capture the dynamics of pH in yoghurt fermentation.

### 3.3.1. Model extension for yoghurt acidification

In order to extend the model developed in Section 2.3 to yoghurt acidification, which is a notably slower process, we have used two data sets from literature. The first set is taken from the work of Aguirre-Ezkauriatza et al. (2008), where they have studied the effect of mixing in yoghurt fermentation in two types of reactors. In this case, the selected data comprises controlled experiments performed under constant temperature whose samples were taken at three different locations. The reader is referred to the paper of Aguirre-Ezkauriatza et al. (2008) for more details. In these experiments the pH was recorded for 5 h until the pH reached a value of 4.5. The dynamics of pH are displayed in Fig. 11. The second data set was taken from Soukoulis et al. (2007) where the authors monitored viscosity and pH in the industrial yoghurt manufacture. The data set used here corresponds to skimmed milk (to be in accordance with the model presented in Section 2.3), where the acidification of yoghurt was recorded for 6.5 h, at which point it had reached plateau pH of 4.5–4.6. The dynamics of pH are shown in Fig. 11.

The pH profiles for yoghurt acidification from both experimental data sets are described by three stages: (i) lag phase (slow decrease of pH), (ii) logarithmic stage (rapid decrease of pH), and (iii) a slowdown stage where pH settles into an equilibrium point. On the other hand, the experiments of milk acidification, which have shown only two stages: (1) acid addition and (2) diffusion. In the case of yoghurt, however, acid is produced in situ by the lactic acid bacteria, which are distributed homogeneously throughout the product. Therefore, the equation that we have considered for the diffusion stage (Eq. (4)) is not required for yoghurt. Another important difference concerns the acid addition rate, which is constant and fixed for the milk acidification experiments presented with sulfuric acid, but for yoghurt is defined by the bacterial conversion of lactose into lactic acid. Hence, it is needed to provide

<figure><img src="image_16.png" /><figcaption>Fig. 12 – Model output for profiles of pH at 40 °C taken from Soukoulis et al. (2007). The symbols are data points whilst the solid line is the model output.</figcaption></figure>

further information to the model to mimic the dynamics of pH.

In order to introduce a simple adaptation approach, we propose to modify the profile of acid addition Q ($L_{acid}/L_{milk}/min$). Several formulations were tested in order to mimic the acid formation by microorganisms, which normally follows a sigmoid function. The best representation is a logistic equation of the form $Q = \frac{c_1}{1+\exp(-c_2t+c_3)}$. However, this function has three new parameters and it is time dependent which makes difficult to extrapolate for different conditions. A simple time-independent relation has been chosen based on pH assuming that pH cannot be higher than the initial pH ($pH_0$). The acid addition is formulated as,

$$Q = \frac{\left(1 - \frac{pH}{pH_0 + c_1}\right)}{t_A} \qquad (15)$$

where $c_1 = 0.01$ corresponding to a small constant to consider that the initial acid addition rate is not zero.

Applying the acid addition profile Q and integrating the model provided in Eqs. (5) and (6), the model output follows the dynamics of pH above 5.1. However, when the model was used with the parameters from Table 2, the results showed a sharp decrease in the pH below 5.1. This suggested that the parameters $\theta_3$ and $\theta_4$ in Eq. (3b) generate an almost linear behavior. Hence, a correction of these parameters was needed to capture the dynamics of the pH. In order to keep the same model structure and the update of parameters as low as possible, we have decided to re-estimate uniquely the parameter $\beta_1$, which affect both $\theta_3$ and $\theta_4$ and it is a sensitive parameter according to the sensitivity analysis presented in last section. This re-estimation resulted in a parameter value of 0.42 for yoghurt, instead of the 0.273 obtained for milk. The larger value suggests a smoother acidification. For the re-estimation of this parameter and the calculation of $c_1$ only the data from Aguirre-Ezkauriatza et al. (2008) is employed. The data from Soukoulis et al. (2007) is considered as a validation set. The main difference between these sets is the initial pH.

The performance of the model with the re-estimated parameter is displayed in Fig. 11 and Fig. 12 for the experiments from Aguirre-Ezkauriatza et al. (2008) and Soukoulis et al. (2007), respectively.

The figures show a good fit of the model predictions with the experimental data with RMSE of 0.297 and 0.325 for the
calibration and validation sets, respectively. The standard
deviation of the calibration was found to be 0.119, whereas
it was 0.052 for the validation set. It should be mentioned that
the only changes were the acid addition rate in Eq. (15) and the
value of the parameter $\beta_1$. Therefore, we can conclude that the
model structure holds for both milk and yoghurt.

4. Conclusions and future work

This work has presented an empirical dynamic model to char-
acterize the dynamics of pH during acidification of milk since
pH is an important indicator of how precipitation proceeds.
The model has been calibrated with experimental data and
validated with an independent data set from literature. Addi-
tionally, the model has been used to connect the description
of steady state behavior of partitioning of some components
during precipitation, such as calcium and phosphates con-
centration with pH and estimate their dynamic behavior.
Although the model can approximate the dynamics of these
components, further experimental work is needed to validate
this extension of the model before its use for experiment
design. The model predicts pH values from 4 to 7 in tem-
peratures between 0 and 50 °C. Higher temperatures may
induce different phenomena that is not captured by the model.
Although different mixing conditions have been tested, the
use of more sophisticated models, e.g., CFD, is encouraged to
explore the impact of mixing in the precipitation of casein.

Even if the model has been developed for milk acidifica-
tion during casein production, we have demonstrated that it
could be adaptable to describe the dynamics of acidification in
other dairy products, e.g., yoghurt. Nonetheless, it is important
to remark that sometimes re-estimation of some parameters
is necessary. In the case of acidification of yogurt, only one
parameter has been updated. This implies that the structure of
the model is robust, but more data sets are required to provide
parameter values that can consider other cases. Re-estimation
of the model parameters is thus necessary if different con-
ditions are used, especially for higher initial pH values and
higher temperatures.

The empirical model can be used in optimization problems
to define the total amount of acid to be added or to calculate
incubation times. Other examples of its utilization could be
understanding digestion at different ages, where the uptake
of nutrients in the stomach is dependent on the acid used
to precipitate the nutrients. Notwithstanding, these examples
would require further experimental validation to assess the
performance of the model.

Declaration of competing interest

The authors report no declarations of interest.

Acknowledgments

This work is part of the INSPEC project (Integrating Sen-
sor Based Process Monitoring and Advanced Process Control)
managed by the Institute for Sustainable Process Technology
(ISPT) and co-funded by TKI-E & I with the supplementary
grant ‘TKI-Toeslag’ for Topconsortia for Knowledge and Inno-
vation (TKI’s) of the Ministry of Economic Affairs and Climate
Policy.

**Appendix A. Supplementary data**

Supplementary data associated with this article can be
found, in the online version, at https://doi.org/10.1016/
j.fbp.2021.04.010.

**References**

Aguirre-Ezkauriatza, E.J., Galarza-González, M.G., Uribe-Bujanda, A.I., Ríos-Licea, M., López-Pacheco, F., Hernández-Brenes, C.M., Alvarez, M.M., 2008. Effect of mixing during fermentation in yogurt manufacturing. J. Dairy Sci. 91 (12), 4454–4465.

Bringe, N.A., Kinsella, J.E., 1990. Acidic coagulation of casein micelles: mechanisms inferred from spectrophotometric studies. J. Dairy Res. 57 (3), 365–375.

De Brabandere, A.G., De Baerdemaeker, J.G., 1999. Effects of process conditions on the pH development during yogurt fermentation. J. Food Eng. 41 (3), 221–227.

de Oliveira, M.N., 2014. Fermented milks: fermented milks and yogurt. In: Encyclopedia of Food Microbiology, second edition. Elsevier Inc., pp. 908–922.

Hofland, G.W., Berkhoff, M.R., Witkamp, G.J., Van der Wielen, L.A., 2003. Dynamics of isoelectric precipitation of casein using sulfuric acid. AIChE J. 49 (8), 2211–2223.

Hofland, G.W., Van Es, M., Van Der Wielen, L.A., Witkamp, G.-J.J., van Es, M., van der Wielen, L.A.M., Witkamp, G.-J.J., 1999. Isoelectric precipitation of casein using high-pressure CO₂. Ind. Eng. Chem. Res. 38 (12), 4919–4927.

Holt, C., 2004. An equilibrium thermodynamic model of the sequestration of calcium phosphate by casein micelles and its application to the calculation of the partition of salts in milk. Eur. Biophys. J. 33 (5), 421–434.

Holt, C., Dalgleish, D.G., Jenness, R., 1981. Calculation of the ion equilibria in milk diffusate and comparison with experiment. Anal. Biochem. 113 (1), 154–163.

Jablonka, M.S., Munro, P.A., Duffy, G.G., 1988. Use of light scattering techniques to study the kinetics of precipitation of mineral acid casein from skim milk. J. Dairy Res. 55 (2), 179–188.

Kristo, E., Biliaderis, C.G., Tzanetakis, N., 2003. Modelling of the acidification process and rheological properties of milk fermented with a yogurt starter culture using response surface methodology. Food Chem. 83 (3), 437–446.

Lucey, J.A., 2017. Formation, structural properties, and rheology of acid-coagulated milk gels. Cheese: Chemistry, Physics and Microbiology, vol. 1., fourth edition Elsevier Inc., pp. 179–197.

Lucey, J.A., Singh, H., 1997. Formation and physical properties of acid milk gels: a review. Food Res. Int. 30 (7), 529–542.

Lucey, J.A., Singh, H., 2003. Acid coagulation of milk. In: Advanced Dairy Chemistry—1. Proteins. Springer US, pp. 1001–1025.

Mekmene, O., Le Graët, Y., Gaucheron, F., 2009. A model for predicting salt equilibria in milk and mineral-enriched milks. Food Chem. 116 (1), 233–239.

Mekmene, O., Le Graët, Y., Gaucheron, F., 2010. Theoretical model for calculating ionic equilibria in milk as a function of pH: comparison to experiment. J. Agric. Food Chem. 58 (7), 4440–4447.

Mudgil, D., Barak, S., 2019. Dairy-based functional beverages. In: Milk-Based Beverages: Volume 9: The Science of Beverages. Elsevier, pp. 67–93.

Sarode, A.R., Sawale, P.D., Khedkar, C.D., Kalyankar, S.D., Pawshe, R.D., 2015. Casein and caseinate: methods of manufacture. In: Encyclopedia of Food and Health., pp. 676–682.

Sfakianakis, P., Tzia, C., 2014. Conventional and innovative processing of milk for yogurt manufacture; development of texture and flavor: a review. Foods 3 (1), 176–193.

Soukoulis, C., Panagiotidis, P., Koureli, R., Tzia, C., 2007. Industrial yogurt manufacture: monitoring of fermentation process and improvement of final product quality. J. Dairy Sci. 90 (6), 2641–2654.