Food Research International 161 (2022) 111745 Contents lists available at ScienceDirect Food Research International journal homepage: www.elsevier.com/locate/foodres Developing an optical backscatter method for determining casein micelle particle size in heated milk Heather Taterka, Anna Zamora *, Manuel Castillo * Centre d’Innovació, Recerca i Transferència en Tecnologia dels Aliments (CIRTTA), Departament de Ciència Animal i dels Aliments, Facultat de Veterinària, Universitat Autònoma de Barcelona, 08193 Bellaterra, Cerdanyola del Vallès, Barcelona, Spain A R T I C L E I N F O A B S T R A C T Keywords: A plethora of different factors, such as heat treatment, pH, soluble calcium and phosphate concentrations, Casein micelle size colloidal calcium phosphate, ionic strength, redox potential, etc., affect functionally of critical milk components Whey protein denaturation such as casein micelles, fat globules and whey proteins. These physicochemical changes induce fat- or protein- Milk protein interactions that would be associated to changes in particle size that might be revealed using light In-line Process monitoring backscatter measurements. We hypothesized that inline, simple, low-cost light backscatter measurements might Optical sensor have the potential to provide functionally related information, representing an interesting opportunity for Near infrared spectroscopy process control. Casein micelle particle size and near infrared light backscatter spectra were measured in milks heat treated at 80 and 90 ◦ C and pH 6.3, 6.7 and 7.1 in order to obtain prediction models for estimating changes in casein micelle particle size during milk heat treatment. Light intensity was measured over a spectral range of 200–1100 nm using a simple optical backscatter sensor and was implemented into models for particle size predictions as a function of heat treatment temperature and pH. Models which included an exponential factor containing a ratio of two specific wavebands were found to improve R2 when compared to single wavelength models. The best model exhibited an R2 of 0.993 and SEP of 2.36 nm. The developed prediction models show promise for in-line monitoring of whey protein denaturation and casein micelle particle size. 1. Introduction et al. (1993) that NIR spectroscopy could be used to successfully quan­ tify protein, fat and total casein in goatś milk. NIR spectroscopy has also In-, on– and at-line monitoring of liquid milk products has imple­ been used to develop prediction equations for the determination of a mentation potential in the dairy industry for process control and time number of milk attributes with in-line application potential. For and cost savings. To date, several applications have been developed to example, for the differentiation of different heat treatments in milk that end, but it is still plenty of room for innovative process analytical (Downey et al., 1990), detection of adulteration by the addition of whey technologies that could provide useful, real-time information for process powder to milk powders (Giangiacomo et al., 1991) or addition of NaCl control optimization and improvement. In that sense, the monitoring of and skim milk powder to milk (Pedretti et al., 1993). Sørensen and various milk product manufacturing steps such as determination of milk Jepsen (1997) used NIR spectroscopy to detect cheese defects as a result gelation mechanisms and changes during the cheese making process has of Clostridium tyrobutyricum. been investigated using light scatter techniques (Castillo et al., 2005c; Development of prediction models has also been accomplished using Castillo et al., 2006; Fagan et al., 2007). Near infrared (NIR) spectros­ NIR and light scatter techniques to predict cheese making characteristics copy has been widely used at- or even in-line in milk to determine ab­ at real time. Castillo, et al. (2000) used NIR light backscatter for sorption in the infrared region (780–2500 nm) of bonds and chemical developing simple prediction models using only optical parameters for groups to quantify various milk components (García, 2004). Robert et al. the determination of cutting and clotting time. Further developments (1987) was able to establish specific wavelengths corresponding to have allowed the in-line transformation of the NIR light backscatter changes in fat, protein and lactose content in milk, however noted that response during milk coagulation into real time and continuous pre­ some interference was observed as a result of large water absorption and diction of the gel elastic modulus during cheese manufacturing (Arango light scattering of fat particles. As well, it was observed by Diaz-Carrillo and Castillo, 2018) and pH during yogurt fermentation (Arango et al., * Corresponding authors. E-mail addresses: anna.zamora@uab.cat (A. Zamora), manuel.castillo@uab.es (M. Castillo). https://doi.org/10.1016/j.foodres.2022.111745 Received 9 March 2022; Received in revised form 19 July 2022; Accepted 25 July 2022 Available online 29 July 2022 0963-9969/© 2022 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by- nc-nd/4.0/). H. Taterka et al. Food Research International 161 (2022) 111745 2020). Light sidescatter and transmission were also used to estimate at 21 ◦ C using 0.5 M HCl or 0.5 M NaOH. re-adjusted milk samples were whey fat concentration, where sidescatter produced models with a allowed to equilibrate for 2 h in a dark place, before final pH reading and higher R2 (greater than0.95) than corresponding transmission models minor re-adjustments. (less than0.5) (Castillo et al., 2005b). An in-line syneresis sensor was A stainless-steel plate (6 mm slot) containing 80 mL of milk was developed by Castillo, et al. (2007), which was applied by Fagan, et al. placed into a water bath (OvanTherm, Suministros Grupo Esper, S.L., (2008) to model cheese manufacturing indices such as whey fat, curd Badalona, Spain) controlled using an OvanTherm TC00 unit (resolution yield and curd moisture content using light backscatter sensor technol­ 0.1, stability ± 0.1 ◦ C) (Suministros Grupo Esper, S.L.) at either 80 or ogy achieving R2 values of 0.90 or more for the prediction models. In 90 ◦ C for 10 min. Samples took less than 30 s to reach the target tem­ particular, their work was mainly aimed toward the determination of perature and heat treatment was stopped by removing the heat plate and curd moisture content as a function of time in order to improve final rapidly placing it in a 0 ◦ C ice-water bath for 3 min. Milk was then moisture content during cheese making. Initial models were determined transferred to a test tube, refrigerated at 4 ◦ C for no more than 2 days, from a total of 40 parameters; however successful models utilized a and re-equilibrated to 21 ◦ C before analysis. combination of the parameters: temperature, percentage of protein, milk fat, and milk solids as well as milk fat protein ratio and light backscatter intensity ratios, which contain coagulation and syneresis information. 2.3. Casein micelle particle size Other studies have investigated the relation between light scatter and whey protein denaturation and the subsequent attachment of whey Particle size z-average was measured at 20 ± 0.5 ◦ C using dynamic proteins (WPs) to the surface of the casein micelle (Lamb, Payne, Xiong, light scattering set to 90◦ and refractive index set to 1.471 with a Mal­ and Castillo, 2013), however the strong pH effect of this mechanism has vern Zetasizer 4 (Malvern Instruments ltd., Malvern, Worcs., UK). not been thoroughly investigated for optical sensor development. Samples were suspended in Ca/imidazole buffer (20 mM-imidazole, 5 The authors hypothesize that particle size may be used as an indirect mM-CaCl2, 30 mM-NaCl, pH 7.0) at a concentration of 1:2 in order to measure of whey protein denaturation in heat treated milk. Particularly suspend the casein micelles and allow stability during measurements because the increase of particle size has been shown to be a good indi­ (Anema, 1997; Anema & Li, 2003; Anema et al., 2004). Measurements cator of the extent of binding of denatured whey proteins to the surface were 3 min each and an average of 3 measurements was taken for each of the casein micelle (Anema and Li, 2003, Taterka and Castillo, 2015). representative sampling. However, this mechanism is affected to a large degree by milk pH. The maximum attachment level of denaturated whey protein to the casein 2.4. Light backscatter measurements micelle surface has been found to occur at pH ~ 6.3, whereas the for­ mation of soluble whey protein aggregates is the preferential mechanism Full detail on the optical system used is given in Taterka and Castillo at higher pH (maximum at pH 7.1) (Vasbinder and de Kruif, 2003; (2015). Briefly, it consisted of a light scatter probe (Reflectronics Inc, Donato and Guyomarc’h, 2009; Kethireddipalli et al., 2010; Taterka and Lexington, KY, USA) connected through two 600-µm UV-NIR optical Castillo, 2015). Moreover, milk pH also affects the light backsctter signal fibers to a tungsten halogen light source (LS-1; Ocean Optics, Inc. as it has been exhibited that bound/soluble aggregate formation is Dunedin, FL, USA) and a high-resolution spectrometer (Model HR4000; highly dependent on pH and corresponds significantly to the intensity of Ocean Optics, Inc.). Temperature control was achieved by circulating light backscatter signal (Taterka and Castillo, 2015). Another important water at 25 ◦ C inside double jacketed walls below the sample well. factor is the temperature of heat treatment, which significantly increases Spectra Suite software (Ocean Optics, Inc.) was used to collect intensity the amount of denatured whey protein, resulting in bound and/or sol­ spectra in the range 200–1100 nm. For single wavelength prediction uble serum aggregates, depending on pH. It should be noted that, whey models, light backscatter maximum intensity (570 nm) was recorded protein denaturation plays a quite significant role in milk and milk and used for further analysis. In addition, an alternative technique of products from a functional point of view. For instance, it is well-known waveband ratios for modeling proposed by Lamb, Payne, Xiong, and that intensive whey protein denaturation impairs milk coagulation Castillo (2013) was employed (see section 3.2.1 for further detail). during cheese manufacturing while, contrarily, it improves yogurt texture and whey retention. Additionally, intense thermal denaturation is required for thermal stabilization of milk prior to sterilization of 2.5. Statistical analysis unsweetened condensed milk. Since the previous models obtained by Lamb et al. (2013) did not take into account the pH dependence of the Statistical analysis and interpretation of data was accomplished binding reaction, the objective of this study was to look further into using “Statistical Analysis System” (SAS, version 9.2, SAS Institute Inc., modeling the light scatter signal with respect to the changes in particle Cary, NC, USA). Single wavelength models (models 1–5) were developed size in heat treated milk, while taking into consideration the effect of using the PROC REG function. The development of models 1–5 is more milk pH. thoroughly described in Section 3.1. Models 6–10 utilized a technique of waveband ratios using a 2. Materials and methods grouping of two wavebands portions of the spectra, determined by SAS code from Lamb et al. (2013). Wavebands were defined as 35 nm por­ 2.1. Experimental design tions of the spectra, thus all grouping of waveband ratios yielded 27 distinct combinations. Prior studies on the behavior of particle size, it The experiment consisted of a 3 × 2 factorial design with three pH has been considered that in a nonlinear regression model the response values (6.3, 6.7 and 7.1) and two heat treatments (80 and 90 ◦ C for 10 variable (particle size) is an exponential function as a function of pH min). Each treatment was replicated three times. (Taterka and Castillo, 2018). The NLIN function was used to form nonlinear regression models, assuming that the response variable fol­ 2.2. Sample preparation lows an exponential model. From this, waveband ratios were developed with the objective of modeling the radius of the casein micelle versus pH Low-heat skim milk powder (low-heat, spray-dried skim milk pow­ and temperature. Models were assessed based on the determination der; pH = 6.5, solubility = 99 %, WPNI ≥ 7 mg g− 1, 800 cfu g− 1) was coefficients (R2) and the standard errors of prediction (SEP) for the re­ supplied by Chr. Hansen SL (Barcelona, Spain). Milk was reconstituted gressions between the predicted and actual values. Supplementary in­ with Milli-Q water to a final solids content of 12 % (w/w). Reconstituted formation regarding the development of wavelength ratios can be found milk, initially at pH 6.55 (±0.015), was pH-adjusted to 6.3, 6.7 and 7.1 in Section 3.2.1. 2 H. Taterka et al. Food Research International 161 (2022) 111745 3. Results and discussion however those models which did not show a large difference in R2 be­ tween quadratic and linear models were simplified into the linear model 3.1. Single wavelength prediction models form. Also, in some “integrated” models, temperature (T) was used for the model as the R2 value increased markedly in some models when A model for particle size as a function of light backscatter intensity at including the T variable. Models in which the addition of T did not in­ 570 nm (single wavelength model) was developed. Fig. 1 shows the crease the R2 value to a large extent were simplified into temperature trends for both light backscatter at 570 nm and z-average particle size as combined models in order to maintain degrees of freedom (DF) for the a function of pH and temperature. Models were “separated” between model. Table 1 shows the various models that were developed and their heat treatment temperatures (80 and 90 ◦ C), temperature “integrated” summary statistics. Temperature “separated” models and temperature or temperature ¨combined¨. Temperature “integrated” models were “integrated” models were developed to predict light backscatter developed as a function of temperature, whereas temperature combined maximum intensity as a function of pH and temperature (model 1 and 2) used all data points but did not include temperature in the prediction (Fig. 2) and particle size z-average as a function of pH and temperature equation. Two types of models were obtained: linear models with the (model 3 and 4) (Table 1) (Fig. 3). Model 1 utilizes a simple linear model form (y = ax + b) and quadratic models in the form (y = ax2 + bx + c). As using only pH as a dependent variable, yet still results in high R2 values a first modeling approach, both light backscatter and particle size were for both 80 and 90 ◦ C models (R2 = 0.979 and 0.988, respectively) modeled as a function of pH and/or T as follows: (Table 1) (Fig. 2A). In the case of the temperature “integrated” model (model 2) (Fig. 2B), all temperature data points were considered (DF = LB = β0 + β1 pH (1) 15) and a high R2 value of 0.984 was obtained. For the particle size temperature “separated” models (model 3) (Fig. 3A) R2 values of 0.881 LB = β0 + β1 pH + β2 T (2) and 0.992 were obtained for 80 and 90 ◦ C models, respectively, whereas PS = β0 + β1 pH + β2 pH 2 (3) the temperature “integrated” model (model 4) (Fig. 3B) had a slightly lower R2 value of 0.791. In general, it can be seen that light backscatter PS = β0 + β1 pH + β2 pH 2 + β3 T (4) showed a consistent linearly decreasing trend with pH (Fig. 1A), whereas changes in particle size tended to show a more quadratic where LB is the light backscatter intensity (bits) measured at the response as a function of pH, where 90 ◦ C samples showed greater maximum intensity wavelength encountered within the wavelength changes with pH than 80 ◦ C samples (Fig. 1B). range between 200 and 1100 nm (i.e., 570 nm), PS is the particle size z- Trends for changes in particle size and light backscatter maximum average (nm), T is the heat treatment temperature (◦ C), and β0-3 are intensity were investigated and prediction models were developed with regression coefficients. high R2 values as a function of milk pH (Models 1–4). As the technique to Light backscatter and particle size experimental data was fit to both measure particle size can be time consuming and costly, the light linear and quadratic mathematical models to estimate regression co­ backscatter method could be a good alternative to particle size efficients and summary statistics using the PROC REG function in Sta­ tistical Analysis System (SAS version 9.2, SAS Institute Inc., Cary, NC, USA). In all cases, higher R2 values were observed in quadratic models, Table 1 Predictive models 1 to 4 for light-backscatter intensity (LB) and particle size z- average (PS). Model Prediction Temperature DF Regression R2 SEP equation err coef. 1 LB = β0 + 80 ◦ C 7 β0 = -5.11 • 0.979 320 β1 pH 104 β1 = 5.95 • 103 90 ◦ C 7 β0 = -5.07 • 0.988 235 104 β1 = 5.82 • 103 2 LB = β0 + Integrated 15 β0 = 4.68 • 0.984 273 β1 pH + β2 T 104 (80 and β1 = -5.88 • 90 ◦ C) 103 β2 = 47.7 3 PS = β0 + 80 ◦ C 6 β0 = 1.97 • 0.881 3.8 β1 pH + β2 pH2 103 β1 = -526 β2 = 37.4 90 ◦ C 6 β0 = 7.87 • 0.992 3.2 103 β1 = -2.23 • 103 β2 = 161 4 PS = β0 + Integrated 14 β0 = 4.78 • 0.791 11.7 β1 pH + 103 β2 pH2 + β3 T (80 and β1 = -1.38 • 90 ◦ C) 103 β2 = 99.1 β3 = 1.60 N = 18; DF err, degrees of freedom for error; coef.: coefficients; β0-3, prediction Fig. 1. (A) Light backscatter maximum intensity (LB) as a function of pH at 80 coefficients; LB, light backscatter intensity (bits); PS, particle size z-average and 90 ◦ C (N = 18); (B) Particle size z-average as a function of pH at 80 and (nm); R2, determination coefficient; SEP, standard error of prediction for the 90 ◦ C (N = 18). model (bits for LB and nm for PS). 3 H. Taterka et al. Food Research International 161 (2022) 111745 Fig. 2. (A) Model 1: Light backscatter maximum intensity modeled as a func­ tion of pH at 80 and 90 ◦ C (N = 18); (B) Model 2: Light backscatter maximum Fig. 3. (A) Model 3: Particle size z-average modeled as a function of pH at 80 intensity modeled as a function of pH and temperature (N = 18). and 90 ◦ C (N = 18); (B) Model 4: Particle size z-average modeled as a function of pH and temperature (N = 18). measurements in milk as it is inexpensive and nondestructive, and can easily be implemented for in-line application. Initial work by Anema and Table 2 Li (2003) presented a strong correlation of denatured whey proteins Predictive model 5 for particle size z-average (PS). attached to the casein micelle and changes in particle size. As well, previous work has shown good correlation between both the light Model Prediction Temperature DF Regression R2 SEP equation err coef. backscatter maximum intensity (r = 0.77) (Taterka and Castillo, 2015) and bound whey protein (r = 0.70) (Taterka and Castillo, 2018) to 5 PS = β0 + 80 ◦ C 6 β0 = 219 0.847 4.4 changes in casein micelle particle size. A model of particle size as a β1 LB + β2 LB2 β1 = -0.0221 function of light backscatter intensity would enable in-line monitoriza­ β2 = 1.18 • tion of the extent of changes in denaturation in heat treated milk. Since 10-6 in-line process control shows promise toward optimization of the 90 ◦ C 6 β0 = 572 0.992 3.1 manufacturing process of milk and milk products such as cheese and β1 = -0.0899 yogurt, a simple, non-invasive, quick and inexpensive technique β2 = 4.41 • 10-6 providing real time information about the binding of denatured whey Combined 15 β0 = 479 0.825 10.3 proteins to the casein micelle and/or soluble aggregate formation would (80 and β1 = -0.0725 be highly advantageous to the dairy industry. Note that different pro­ 90 ◦ C) cessing factors inducing particle interactions in milk are relevant for β2 = 3.59 • 10-6 milk functionality and would be associated to changes in particle size that might be real-time revealed using light backscatter measurements. N = 18; DF err, degrees of freedom for error; coef.: coefficients; β0-3, prediction And a relatively simple equation for the prediction of particle size as a coefficients; LB, light backscatter intensity (bits); PS, particle size z-average function of the light backscatter maximum, was developed: (nm); R2, determination coefficient; SEP, standard error of prediction for the model (bits for LB and nm for PS). PS = β0 + β1 LB + β2 LB2 (5) 2 Model 5 (Table 2) used a quadratic form and yielded high R values which has a higher R2 and lower SEP value than Model 1 temperature in the case of 80 and 90 ◦ C models (R2 = 0.847 and 0.992, respectively) separated at 80 ◦ C). However, all temperature “integrated” models (Fig. 4A) and temperature “combined” models (R2 = 0.825) (Fig. 4B). As maintained an R2 value higher than or equal to 0.791 and all 80 and it can be observed in Table 1 and 2, temperature separate models had 90 ◦ C models with R2 greater than 0.847 (Tables 1 and 2). higher R2 and lower standard error of prediction (SEP) in nearly all Other authors have looked at various characteristics of milk using models (except in the case of the temperature “integrated” Model 2, NIR (780–2500 nm) or mid-infrared (MIR) (2500–15000 nm) 4 H. Taterka et al. Food Research International 161 (2022) 111745 In addition, a technique of ratios has been used to decrease the number of variables used in the model by combining certain parameters, such as in the case of Fagan et al. (2008), which incorporated a non-optical milk- fat protein ratio into the model for curd moisture content. Therefore, in order to pursue a higher level of accuracy for predictive models, a method of waveband ratios was also tested. Note that Fagan et al. (2008) tested Partial Least Squares (PLS) regression analysis for use in optical sensor development and found little improvement in models when compared to individual wavelength and waveband ratio models. In addition, complex model development techniques are inherently more complicated and expensive as they require sensors that use multiple wavelength analysis and, as a result, it was opted out of using any more than two wavelength combinations for the development of models in the present study. Initially, wavebands of 15, 25 and 35 nm portions of the light backscatter spectra (from 200 to 1100 nm) were investigated in pre­ liminary models, however no major differences were observed. As a result, wavebands were defined as an average of the intensity for each 35 nm portion of the spectra. Thus, a total of 27 waveband regions were obtained. The wavebands were then combined, using each waveband combination in both the numerator and denominator, to come up with all possible ratios. Predictive models were determined using the pre- selected waveband ratios that were found to best correspond to changes in particle size as a function of pH (i.e., assuming an exponential function). Table 3 summarizes the top 10 ratios, which exhibited the highest R2 values used for the model development and based on the total R2 value from a sum of the four models which used waveband ratios (Models 6 (80 ◦ C), 6 (90 ◦ C), 7 and 10). Fig. 5 illustrates the typical light backscatter profile obtained from 200 to 1100 nm scans. Highlighted are the two sections of wavebands which were used to obtain the best models. Wavebands 6 and 7 (range of 388–458 nm) were selected as Fig. 4. (A) Model 5 temperature separated: Particle size z-average modeled as a numerator values that corresponded with models showing the highest R2 function of light backscatter maximum intensity (N = 18); (B) Model 5 (tem­ values. The number of denominators that corresponded to the highest R2 perature combined): Particle size z-average modeled as a function of light values for the models was larger, wavebands 15–20 (range of 703–878 backscatter maximum intensity (N = 18). nm). All ratios in the top 10 models fell within the previously described numerator and denominator ranges. Numerator values that produced spectroscopy techniques (Iñón et al., 2004; Wu et al., 2011) by incor­ the 10 highest R2 represent a region with relatively little change among porating a combination of chemometric techniques for analyzing the the different samples at various experimental conditions (pH and tem­ range of spectral data in order to form useful prediction models. The NIR perature of heat treatment), whereas denominator values tended to lie in region is widely used in milk analysis and therefore the spectral regions a region with more notable changes (Fig. 5). The technique of using which correspond to moisture content, milk fat, protein, lactose and numerator and denominator values presented as a ratio and/or as in­ other milk components are well-characterized (Kamishikiryo-Yamashita dividual values was further detailed in Section 2. et al., 1994; Laporte and Paquin, 1999; Robert et al., 1987; Tsenkova et al., 1999). However, it should be noted that although the light 3.2.2. Models backscatter spectral range (200–1100 nm) has been studied with respect Using the ratio (R6,18) which yielded the highest R2 for particle size, to milk particulate components, such as the fat globule and casein prediction models were developed for particle size (z-average) and light micelle (Castillo et al., 2005a,b; Fagan et al., 2008; Lamb et al., 2013), backscatter waveband ratio (R6,18) as a function of pH in temperature regions that correspond to specific characteristic and or physicochem­ “separated” and temperature “integrated” models (Table 4). In ical changes in milk have yet to be fully characterized. Table 3 R2 top 10 ratios for modeling particle size. 3.2. Waveband ratio predictive models ratio numerator denominator R2 R2avg R2 R2avg PS f(LB PS f(LB LB ratio LB ratio 3.2.1. Waveband ratio selection ratio) ratio) f(pH,T) f(pH,T) In the present study, a good modelization for particle size as a 6_18 388 808 0.979 0.952 0.991 0.993 function of light backscatter maximum intensity (model 5) has been 6_17 388 773 0.979 0.951 0.992 0.993 observed. However Lamb et al. (2013) found improved predictions 6_16 388 843 0.978 0.949 0.993 0.993 when using an alternative technique of waveband ratios for modeling 7_16 388 738 0.977 0.949 0.992 0.993 whey protein denaturation as a function of light backscatter intensity. 7_15 458 773 0.977 0.948 0.992 0.993 6_19 458 738 0.980 0.954 0.987 0.989 Lamb et al. (2013) defined waveband ratios as the average intensity at 7_17 458 703 0.978 0.951 0.989 0.991 each 25 nm portion of the light backscatter spectra, and created wave­ 6_15 458 808 0.976 0.948 0.991 0.992 band ratio combinations, which were used to form predictive models. 6_20 388 878 0.981 0.955 0.982 0.985 Implementing a technique of ratios, or a combination of parameters, has 7_18 388 703 0.976 0.951 0.986 0.989 been used successfully in prediction models found in the literature R2, determination coefficient; R2avg, average determination coefficient of 80 and (Castillo et al., 2000; M. Castillo et al., 2005a,b; Fagan et al., 2008), 90 ◦ C; LB, light backscatter intensity (bits); T, temperature of heat treatment which in certain cases was found to improve the accuracy of modeling. (◦ C); PS, particle size z-average (nm); numerator and denominator (nm). 5 H. Taterka et al. Food Research International 161 (2022) 111745 Fig. 5. Light backscatter spectra showing selected numerator and denominator regions used for modeling. Table 4 In Models using top 1 ratio R6,18 for prediction of and an exponential factor for prediction of particle size z-average (PS). Id Model Prediction equation Temperature DF err Regression coef. R2 SEP 6 In 80 C ◦ 7 β0 = 0.454 0.991 3.69 • 10-3 = βo + β1 pH Id β1 = -0.103 90 ◦ C 7 β0 = 0.488 0.996 2.57 • 10-3 β1 = -0.106 7 In Integrated 15 β0 = 0.383 0.993 0.0377 = βo + β1 pH + β2 T Id (80 and 90 ◦ C) β1 = -0.104 β2 = -1.03 • 10-3 8 PS = βo + e (α0 +α1 pH) 80 C ◦ 24 β0 = 115 0.871 3.49 α0 = -25.2 α1 = 3.52 90 ◦ C 21 β0 = 117 0.988 3.09 β1 = -41.2 β2 = 5.87 9 PS = βo + e(α0 +α1 pH+α2 T) Integrated 47 β0 = 117 0.976 3.32 (80 and 90 ◦ C) α0 = -29.7 α1 = 5.78 α2 = 0.121 N = 18 for In/Id models, N = 51 for PS models; DF err, degrees of freedom for error; coef.: coefficients; β0-2 and α0-2, prediction coefficients; R2, determination co­ In efficient; SEP, standard error of prediction for the model (dimensionless for and nm for PS). Id comparison to models that were developed to determine R6,18 (models numerator, Id is the light backscatter intensity at the pre-selected de­ 6–7; 6, 7) (Fig. 6), models for z-average (models 8–9; 8, 9) (Fig. 7) nominator, PS is the particle size z-average, T is the heat treatment incorporated an exponential factor into equations. As mentioned pre­ temperature, β0-2 and α0-2 are prediction coefficients. viously and visualized in Fig. 1B, particle size z-average tends to follow a A previous work has shown that changes in light backscatter signal more exponential curve with respect to pH when compared to light and particle size have been found to be highly correlated to pH (Taterka backscatter maximum intensity, which follows a linear trend with pH and Castillo, 2015). Thus in model development, it is reasonable that pH (Fig. 1A). Equations for predictions are as follows: has been found to be a critical factor. Castillo et al. (2006) developed models for cutting time prediction in cottage cheese manufacture also In = β0 + β1 pH (6) using pH as a variable for predictions, whereas reasonably accurate Id prediction models using pH as a variable have been developed by Tof­ In fanin, De Marchi, Lopez-Villalobos, and Cassandro (2015) in order to = β0 + β1 pH + β2 T (7) determine characteristics of milk quality and coagulation properties Id using MIR analysis. As pH is a useful tool in the present study for the PS = β0 + e(α0 +α1 pH) (8) determination of particle size and light backscatter intensity, it may be complementary to utilize a simple technique of milk pH measurement PS = β0 + e(α0 +α1 pH+α2 T) (9) for the determination of changes in the milk matrix with heat treatment. Unfortunately, simple laboratory pH meters generally contain glass where In is the light backscatter intensity at the pre-selected components and other reagents that are not approved for in-line use in 6 H. Taterka et al. Food Research International 161 (2022) 111745 the food industry (Wesstrom, 1992). A preliminary work by Arango et al. (2020) has exhibited that optical sensors are useful for the replacement of pH meters during acid-coagulation of milk in which pH predictions as a function of light scatter techniques were developed successfully with great accuracy (R2 greater than 0.99). Thus, the incorporation of an optical technology for the determination of pH-specific denaturation mechanisms may be of interest for investigation in the present study. As a result, models for particle size using the light backscatter spectrum have been developed with good correlations. Utilizing the ratio values showed to greatly improve R2 when comparing model 5, which did not use ratio values, (Table 2) (Fig. 4) to model 10, which did (Table 5) (Fig. 8). Model 10 is presented below, utilizing R6,18 in an exponential type equation to model particle size changes. (10) (α0 +α2 IIn ) PS = β0 + e d An interesting finding in prediction models for particle size z-average (Model 8–10) is that the coefficient β0 lies within the range of the initial casein micelle particle radius (~112–120 nm, Tables 4, 5) (Figs. 7, 8). The remainder of the prediction equation includes an exponential factor of other predictors (8–10), which for Model 10 corresponds to a wave­ band ratio. Since binding has been found to be the main contributor in the in­ crease in particle size (Taterka and Castillo, 2015), it may be suggested that the binding reaction follows an exponential growth curve. Thus, using the average initial casein micelle particle size plus some expo­ nential increase as a result of attachment and/or aggregation, it may be possible to predict casein micelle particle size growth and, in turn, denaturation as a function of particle size increase. In the present study, the average radius was 119.7 nm using reconstituted milk from the same Fig. 6. (A) Model 6: Light backscatter ratio R6,18 (dimensionless) modeled as a lot, comparable to average diameters reported in Martin, Williams, and function of pH (N = 18); (B) Model 7: Light backscatter ratio R6,18 (dimen­ Dunstan (2007) which saw an average of 231.0 ± 1.6 nm in recon­ sionless) modeled as a function of pH and temperature (N = 18). stituted milk stirred for 35 min. However, it should be noted that casein micelle average size varies significantly depending on batch and thus this hypothesis should be tested using different milk batches. 3.3. Comparison of single wavelength and ratio models 3.3.1. Comparison of models Two sets of models have been developed and discussed: models developed using the maximum light backscatter intensity (LB) (Tables 1, 2) and models using ratio R6,18 and an exponential factor in z-average models (Table 4, 5). Models from Table 1 are complimentary to models from Table 4 as such: Model 1 and 6, 2 and 7, 3 and 8, 4 and 9; and 5 and 10 in Table 2 and Table 5, respectively. In most cases, using a ratio value and/or adding an exponential factor to particle size prediction models acted to increase the R2 value of the model, except model 3 compared to Table 5 Predictive Model 10 for particle size z-average (PS). Model Prediction Temperature DF Regression R2 SEP equation err coef. 10 PS = β0 + 80 ◦ C 24 β0 = 115 0.911 2.89 In (α0 +α2 ) e Id α0 = -9.22 α1 = 15.6 90 ◦ C 21 β0 = 117 0.993 2.36 α0 = -12.5 α1 = 23.2 Combined 48 β0 = 119 0.979 3.04 (80 and α0 = -17.5 90 ◦ C) α1 = 37.0 Fig. 7. (A) Model 8: Particle size z-average modeled as a function of pH at 80 N = 51; DF err, degrees of freedom for error; coef.: coefficients; β0-2 and α0-2, and 90 ◦ C (N = 51); (B) Model 9: Particle size z-average modeled as a function prediction coefficients; PS, particle size z-average (nm); R2, determination co­ of pH and temperature (N = 51). efficient; SEP, standard error of prediction for the model (nm). 7 H. Taterka et al. Food Research International 161 (2022) 111745 micelle, this gives some insight into aspects of whey protein denatur­ ation in milk. 3.3.2. Analysis based on various portions of R6,18 Since two techniques were implemented for modeling, one using a single intensity at the maximum value (LB; in this section represented by Im), and the other using a ratio of intensities (R6,18), which acted to improve the models in most cases, it was of interest to test whether using only the numerator or denominator of the ratio R6,18 (In or Id, respec­ tively) may also yield good predictions individually. Note that a single wavelength prediction is always more convenient for industrial imple­ mentation of a sensor technology. Thus, in order to confirm the need for ratios in the development of prediction models, model 5 was tested using the maximum intensity (Im), the intensity used in the numerator (In) of R6,18, and the intensity of the denominator (Id) used in R6,18. As sum­ marized in Table 6, it can be seen that both Im (R2 = 0.847 and 0.992) and Id (R2 = 0.779 and 0.986) showed good R2 for both 80 and 90 ◦ C models, respectively. On the other hand, In yielded less reliable R2 values (R2 = 0.465 and 0.773). In temperature “combined” models Id was the model with the highest R2 (0.890) and lowest SEP (8.20 nm), compared to Im and In models (R2 = 0.825; 0.491, SEP = 10.3; 17.6 nm, respec­ tively). As In, in general, shows poor correlation of predicted and actual values in models, it leads us to believe that it may not be necessary for Table 6 Models using individual wavelength values to model particle size as a function of the maximum intensity value (Im = 570 nm), the intensity used in the numerator (In) and the intensity of the denominator (Id). Model Prediction Temperature DF Regression R2 SEP Fig. 8. (A) Model 10: Particle size z-average modeled as a function of light equation err coef. backscatter ratio R6,18 at 80 and 90 ◦ C (N = 51); (B) Model 10 (temperature 5 PS = β0 + 80 ◦ C 6 β0 = 1.18• 0.847 4.36 combined): Particle size z-average modeled as a function of light backscatter β1 Im + β2 Im 2 10-6 ratio R6,18 (N = 51). β1 = 0.0221 β2 = 219 90 ◦ C 6 β0 = 4.41• 0.992 3.12 model 8 where adding an exponential factor gave a slightly lower R2 for 10-6 both 80 (model 3: 0.881, model 8: 0.871) and 90 ◦ C models (model 3: β1 = 0.0899 0.992, model 8: 0.988). This may be due to the fact that there is a large β2 = 572 increase from 80 ◦ C to 90 ◦ C in particle size at pH 6.3. This difference Combined 15 β0 = 479 0.825 10.3 (80 and makes it necessary for temperature “integrated” models to include an β1 = 90 ◦ C) -0.0725 exponential factor (R2 = 0.791 -Table 1- versus 0.976 -Table 4- with β2 = exponential factor), whereas temperature “separated” models using a 3.59•10-6 simpler equation still maintain a high R2 value for the models (R2 = 11 PS = β0 + 80 ◦ C 6 β0 = -2.44• 0.465 8.14 0.881 and 0.992 -Table 1- versus 0.871 and 0.988 -Table 4-). β1 In + β2 In 2 10-5 β1 = -0.110 For the modelization of particle size z-average as a function of light β2 = 49.9 backscatter intensity, model 5 and model 10 use light backscatter in­ 90 ◦ C 6 β0 = -1.36• 0.773 16.6 tensity values (model 5: LB, model 10: R6,18). Both models show high R2 10-3 values, however it can be seen that the addition of an exponential factor, β1 = -1.99 as well as using R6,18 in model 10, acts to increase R2 in both tempera­ β2 = 847 Combined 15 β0 = 436 0.491 17.6 ture “separate” and “combined” models (R2 for 80 ◦ C; 90 ◦ C and (80 and β1 = -0.920 “combined” models were 0.847, 0.992, 0.825 for model 5: and 0.911, 90 ◦ C) 0.993, 0.979 for model 10, respectively). This effect is also observed in β2 = 6.60• the case of the temperature “integrated” models (Models 4 and 9) in 10-4 12 PS = β0 + 80 ◦ C 6 β0 = 4.18• 0.889 3.72 which the exponential addition greatly increases the R2 (0.791 and β1 Id + β2 Id 2 10-6 0.976, respectively). On the other hand, when comparing Model 3 with β1 = 0.0231 Model 8, the simple quadratic form yields a prediction with a higher R2 β2 = 148 (80 ◦ C: 0.881, 90 ◦ C: 0.992) than when the exponential form is added to 90 ◦ C 6 β0 = 1.36• 0.994 2.60 the equation (80 ◦ C: 0.871, 90 ◦ C: 0.988). In this case, models that are 10-5 β1 = 0.0879 produced at each respective temperature separately, using a linear β2 = 258 model, supports each prediction sufficiently and to a better extent than Combined 15 β0 = 216 0.890 8.20 the exponential model. Nonetheless, there is a good correlation between (80 and β1 = 4.46• actual and predicted values in all models, thus the light backscatter 90 ◦ C) 10-5 technique shows promise toward the prediction of particle size changes β2 = -0.134 as a function of milk pH and temperature treatment. Since changes in N = 18, DF err, degrees of freedom for error; β0-2, prediction coefficients; R2, particle size have been mainly found to be a result of the binding reac­ determination coefficient; SEP, standard error of prediction for the model (nm); tion of denatured whey proteins attaching to the surface of the casein Im: single intensity at the maximum value; In: intensity used in the numerator; Id: intensity of the denominator. 8 H. Taterka et al. Food Research International 161 (2022) 111745 model development. A possible reason for this will be introduced below. FP7-Reintegration-Grant, 2010. FP7-PEOPLE-2010-RG. Proposal N◦ In fact, Id yields even higher R2 values than Im; possibly showing more 268281 DenatureProbe: Development of an optical backscatter sensor for potential in the corresponding portion of the spectrum (808 nm) for the determining thermal denaturation of whey proteins during milk processing). determination of various changes in milk with heat treatment and varying milk pH as compared with the peak of maximum intensity (570 References nm). Previous work in the literature (Lamb et al., 2013) found it useful to form ratios for predictive models however in the present study the Anema, S. G. (1997). The effect of chymosin on κ-casein-coated polystyrene latex particles and bovine casein micelles. International Dairy Journal, 7(8), 553–558. regions of Im and Id used individually appear to be much more important Anema, S. G., Lowe, E. K., & Li, Y. (2004). Effect of pH on the viscosity of heated than In, at least in the case of quadratic models for particle size as a reconstituted skim milk. International Dairy Journal, 14(6), 541–548. function of light backscatter intensity. Models used for comparison Anema, S. G., & Li, Y. (2003). Association of denatured whey proteins with casein micelles in heated reconstituted skim milk and its effect on casein micelle size. (Table 6) do not include an exponential factor, as it was not possible to Journal of Dairy Research, 70, 73–83. model individual intensity values without modifying the exponential Arango, O., & Castillo, M. (2018). A method for the inline measurement of milk gel equation. Since little change in intensity occurs between sample treat­ firmness using an optical sensor. Journal of Dairy Science, 101, 3910–3917. Arango, O., Trujillo, A. J., & Castillo, M. (2020). Inline control of yoghurt fermentation ments in the regions where numerators with the highest R2 for models process using a near infrared light backscatter sensor. Journal of Food Engineering, occur (Table 3), it is suggested that the portion of the spectrum corre­ 277, Article 109885. sponding to the denominators (703–878 nm) is in fact representing the Castillo, M., Fagan, C.C., Payne, F.A., Ó Donnell, C.P., Ó Callaghan, D.J., 2007. Optical measurement of kinetic changes in curd moisture content and whey fat particle size information while numerator values act as a normalization concentration during syneresis in cheese manufacturing. J. of Dairy Sci. 90(Supple. factor for models. Additionally, exponential models tend to yield high R2 1), 272. values and include the coefficient β0, which may be representative of the Castillo, M., Payne, F. A., Hicks, C., & Lopez, M. (2000). Predicting cutting and clotting initial particle size ratio values. Considering these two observations, it is time of coagulating goat’s milk using diffuse reflectance: Effect of pH, temperature and enzyme concentration. International Dairy Journal, 10, 551–562. proposed that development of an in-line sensor using only one wave­ Castillo, M., Payne, F. A., & Lopez, M. (2005a). Preliminary evaluation of an optical length plus some predetermined normalization factor (for example, method for modeling the dilution of fat globules in whey during syneresis of cheese replacing the numerator with a constant) implemented into an expo­ curd. Applied Engineering in Agriculture, 21(2), 265–269. Castillo, M., Payne, F. A., López, M. B., Ferrandini, E., & Laencina, J. (2005b). Optical nential equation should yield a considerably accurate model, and in sensor technology for measuring whey fat concentration in cheese making. Journal of addition, would be less costly and complicated than a dual wavelength Food Engineering, 71, 354–360. optical sensor. Castillo, M., Payne, F. A., & Shea, A. P. (2005c). Development of a combined sensor technology for monitoring coagulation and syneresis operations in cheese making. Journal of Dairy Science, 88(Suppl. 1):142, (Abstr.). 4. Conclusion Castillo, M., Payne, F. A., Wang, T., & Lucey, J. (2006). Effect of temperature and inoculum concentration on prediction of both gelation time and cutting time. Cottage cheese-type gels. Int. dairy J., 16, 147–152. Two prediction models were developed to estimate particle size z- Diaz-Carrillo, E., Muñoz-Serrano, A., Alonso-Moraga, A., & Serradilla-Manrique, J. M. average only as a function of the intensities of a single wavelength/ (1993). Near infrared calibrations for goat’s milk components: Protein, total casein, waveband combinations on the light backscatter spectra (Model 5 -sin­ αs-, β-and k-caseins, fat and lactose. Journal of Near Infrared Spectroscopy, 1, 141–146. gle wavelength- and 10 -wavelength ratio-), both showing good corre­ Donato, L., Guyomarc’h, F., 2009. Formation and properties of the whey protein/ lation between actual and predicted values. However, improvements κ-casein complexes in heated skim milk – A review. Dairy Sci. Technol. 89, 3–29. were found with the incorporation of R6,18 and the addition of an Downey, G., Robert, P., Bertrand, D., & Kelly, P. M. (1990). Classification of Commercial exponential factor (Model 10) compared to the quadratic model using LB Skim Milk Powders According to Heat Treatment Using Factorial Discriminant Analysis of Near-Infrared Reflectance Spectra. Applied Spectroscopy, 44, 150–155. (i.e., Im) (Model 5). As changes in particle size diameter have been found Fagan, C. C., Castillo, M., O’Donnell, C. P., O’Callaghan, D. J., & Payne, F. A. (2008). On- to be primarily a result of the attachment of denatured whey proteins to line prediction of cheese making indices using backscatter of near infrared light. the surface of the casein micelle, this model provides useful information International Dairy Journal, 18, 120–128. Fagan, C. C., Leedy, M., Castillo, M., Payne, F. A., O’Donnell, C. P., & O’Callaghan, D. J. regarding the potential modelization of WP denaturation. Taking into (2007). Development of a light scatter sensor technology for on-line monitoring of account the investigated models, the suggested model for particle size as milk coagulation and whey separation. Journal of Food Engineering, 83, 61–67. a function of light backscatter intensity would include both an expo­ García, J., 2004. Aplicaciones de la tecnología NIRS en la industria agroalimentaria. Nuevas Tecnologías para el Control proceso y Prodroducto en la Industria nential component, an initial intercept which corresponds to an Alimentaria. Ed. Secretariado de Publicaciones en Intercambio Editorial de la approximate value of initial particle size and a light backscatter Universidad de Valladolid. Valladolid. normalized value or a waveband ratio. Further work may give more Giangiacomo, R., Braga, F., & Galliena, C. (1991). Use of near-infrared spectroscopy to detect whey powder mixed with milk powder. Near Infrared Spectroscopy, Making insight into a more improved model. Nonetheless, this experiment yields Light Work., 399–407. useful preliminary information toward the development a comprehen­ Iñón, F. A., Garrigues, S., & de la Guardia, M. (2004). Nutritional parameters of sive model for the determination of whey protein denaturation with the commercially available milk samples by FTIR and chemometric techniques. Analytica Chimica Acta, 513, 401–412. potential for implementation of an in-line optical sensor for in plant Kamishikiryo-Yamashita, H., Oritani, Y., Takamura, H., & Matoba, T. (1994). Protein processed milk. content in milk by near-infrared spectroscopy. Journal of Food Science, 59(2), 313–315. Kethireddipalli, P., Hill, A. R., & Dalgleish, D. G. (2010). Protein interactions in heat- Declaration of Competing Interest treated milk and effect on rennet coagulation. International Dairy Journal, 20(12), 838–843. The authors declare that they have no known competing financial Lamb, A., Payne, F., Xiong, Y. L., & Castillo, M. (2013). Optical backscatter method for interests or personal relationships that could have appeared to influence determining thermal denaturation of β-lactoglobulin and other whey proteins in milk. Journal of Dairy Science, 96, 1356–1365. the work reported in this paper. Laporte, M.-F., & Paquin, P. (1999). Near-Infrared Analysis of Fat, Protein, and Casein in Cow’s Milk. Journal of Agriculture and Food Chemistry, 47, 2600–2605. Acknowledgements Martin, G. J. O., Williams, R. P. W., & Dunstan, D. E. (2007). Comparison of casein micelles in raw and reconstituted skim milk. Journal of Dairy Science, 90, 4543–4551. Pedretti, N., Bertrand, D., Semenou, M., Robert, P., & Giangiacomo, R. (1993). The authors want to thank Ester Boixadera and Anabel Blasco of the Application of an experimental design to the detection of foreign substances in milk. Servei d’Estadística Aplicada at the Universitat Autònoma de Barcelona Journal of Near Infrared Spectroscopy, 1, 174. Robert, P., Bertrand, D., Devaux, M. F., & Grappin, R. (1987). Multivariate analysis for assistance in model development. applied to near-infrared spectra of milk. Analytical Chemistry, 59, 2187–2191. Sørensen, L., & Jepsen, R. (1997). Detection of cheese batches exposed to Clostridium Funding tyrobutyricum spoilage by near infrared spectroscopy. Journal of Near Infrared Spectroscopy, 5, 91. This work was supported by the European Commission (Marie Curie 9 H. Taterka et al. Food Research International 161 (2022) 111745 Taterka, H., & Castillo, M. (2015). The effect of whey protein denaturation on light Tsenkova, R., Atanassova, S., Toyoda, K., Ozaki, Y., Itoh, K., & Fearn, T. (1999). Near- backscatter and particle size of the casein micelle as a function of pH and heat- infrared spectroscopy for dairy management: Measurement of unhomogenized milk treatment temperature. International Dairy Journal, 48, 53–59. composition. Journal of Dairy Science, 82, 2344–2351. Taterka, H., & Castillo, M. (2018). Analysis of the preferential mechanisms of Vasbinder, A. J., & de Kruif, C. G. (2003). Casein–whey protein interactions in heated denaturation of whey protein variants as a function of temperature and pH for the milk: The influence of pH. International Dairy Journal, 13(8), 669–677. development of an optical sensor. International Journal of Dairy Technology, 71(1), Wesstrom, O., 1992. Inline pH Measurement for the Food / Dairy and Beverage Industry. 226–235. Wu, D., Nie, P., He, Y., & Bao, Y. (2011). Determination of Calcium Content in Powdered Toffanin, V., De Marchi, M., Lopez-Villalobos, N., & Cassandro, M. (2015). Effectiveness Milk Using Near and Mid-Infrared Spectroscopy with Variable Selection and of mid-infrared spectroscopy for prediction of the contents of calcium and Chemometrics. Food and Bioprocess Technology, 5, 1402–1410. phosphorus, and titratable acidity of milk and their relationship with milk quality and coagulation properties. International Dairy Journal, 41, 68–73. 10