Food Chemistry 176 (2015) 271–277
Contents lists available at ScienceDirect
Food Chemistry journal homepage: www.elsevier.com/locate/foodchem
Analytical Methods
Monitoring of whey quality with NIR spectroscopy—A feasibility study Sergey Kucheryavskiy a,⇑, Carina Juel Lomborg b a b
Section of Chemical Engineering, Aalborg University, Campus Esbjerg, Denmark Arla Foods Ingredients Group P/S, Viby, Denmark
a r t i c l e
i n f o
Article history: Received 17 July 2014 Received in revised form 16 December 2014 Accepted 20 December 2014 Available online 30 December 2014 Keywords: Whey Near-infrared spectroscopy Chemometrics
a b s t r a c t The possibility of using near-infrared (NIR) spectroscopy for monitoring of liquid whey quality parameters during protein production process has been tested. The parameters included total solids, lactose, protein and fat content. The samples for the experiment were taken from real industrial processes and had a large variability for most of the parameters. Partial Least Squares (PLS) regression was used to make the prediction models based on NIR spectra taken at 30 and 40 °C. Using proper wavelength range allowed to get models for prediction of fat, protein and amount of total solids with very high precision and accuracy. The lactose was found to be the most challenging parameter. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Whey is a liquid product remaining after production of cheese or removal of fat and casein from cow milk. Being considered as a waste product for many years, nowadays it is widely used in the food industry due to its nutritional properties—whey contains about 50% of the nutrients present in milk, including milk sugar (lactose), serum proteins (whey proteins), minerals, a small amount of fat, and most of the water soluble minor nutrients from milk such as vitamins (de Wit, 2001). Production of protein is one of the applications whey is widely used for. There are several techniques to get the whey protein concentrate; one of which is a membrane filtration, where the whey is sieved through a set of filters of controlled pore size (de Wit, 2001). The concentration of proteins in liquid whey at different stages of the process is changing from 3–5% to 25–30% whereas the lactose content is decreasing. The monitoring of the different constituents during the production is first and foremost important from a quality perspective i.e., being able to control the process so that a uniform product, with respect to the protein and dry matter, is obtained from batch to batch is a significant quality parameter. Furthermore it is important financially, that the product contains the amount of protein paid for by the customer. Conventionally the different constituents may be measured applying standard chemical methods, e.g., Kjeldahl analysis (Barbano, Lynch, & Fleming, 1991) for protein, enzymatic kit ⇑ Corresponding author at: Niels Bohrs vej 8, Esbjerg DK 6700, Denmark. Tel.: +45 2787 9818. E-mail address: [email protected] (S. Kucheryavskiy). http://dx.doi.org/10.1016/j.foodchem.2014.12.086 0308-8146/Ó 2014 Elsevier Ltd. All rights reserved.
(Kleyn, 1985) or HPLC analysis for lactose, and Röse-Gottlieb (Crocker et al., 2009) for fat. However all these methods are laborious and time consuming, and thus not really a possibility for process monitoring. In an attempt to get a quicker result applicable for process monitoring and control, the degrees of brix (°Bx) based on the measurement of the refractive index of a sample from the process stream has been implemented in the production. The refractive brix was originally created to determine the amount of sucrose in water, but is now widely being used across industries for determination of dry matter in process streams. In dealing with raw materials that are not standardized and know to fluctuate over the year, the refractive brix however has it limitations for real process monitoring and understanding. As in applying the method only a measure of the total dry matter is obtained and hence not the composition with respect to lactose and protein, thus meaning a shift in the ratio between lactose and protein may not be caught by this method. As a consequence, during the last decade, analysis based on using wet chemistry tends to be replaced by modern spectroscopic methods combined with multivariate data analysis. Thus, midinfrared spectroscopy has been widely accepted as a laboratory standard for the milk nutrient analysis (Karoui & Debaerdemaeker, 2007). Numerous research has been carried out on using NIR spectroscopy for estimating quality parameters of milk and milk products (Bogomolov, Dietrich, Boldrini, & Kessler, 2012; Karoui & Debaerdemaeker, 2007; Tsenkova et al., 1999). However, to our knowledge, very few papers on estimation of liquid whey quality parameters with NIR spectroscopy exist—most of the research in this field deals with whey powders. In (Pouliot,
272
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
Paquin, Martel, Gauthier, & Pouliot, 1997) authors tried to solve the similar problem, however the data quality as well as obtained results are outside the industrial standards. There are also several commercial analyzers based in FT-NIR (see for example, (Bruker Optics, 2009)) that can analyze liquid whey quality, however they work with a narrow range of the parameters, e.g., use different models and settings for low, middle or high amount of protein content and total solids. The main purpose of the present work was to carry out a feasibility study of using NIR spectroscopy for monitoring of major whey quality parameters: total solids, protein, lactose and fat content during whey filtration processes, where the parameters may vary significantly. 2. Materials and methods 2.1. Samples and measurement of quality parameters Samples for the experiments were collected from four different productions in a real industrial process on whey filtering at Danmark Protein factory (Arla Foods Ingredients Group A/S, Videbæk, Denmark). Two 100 ml replicates were taken for each case and transported to a spectroscopic laboratory in plastic vials. Each sample was also delivered to a factory lab where reference values for total protein (Prot), fat (Fat) and lactose (EzL) content as well as amount of total solids (TS102), were evaluated using the following procedures. The total solid content was determined according to ‘‘ISO 6731:2010(E)/IDF 21:2010(E) Milk, cream and evaporated milk. Determination of the total solids content.’’ The protein content was defined as the nitrogen content times a Kjeldahl factor (6.38), and reported as g 100 g 1 sample. The Nitrogen content was determined according to the following standard: ‘‘ISO 8968-3:2004(E)/IDF 20-3:2004 Determination of nitrogen content. Block-digestion method (semi-micro rapid routine method)’’. In the destruction 0.6–0.8 g of sample was weighted and digested at 420 °C for 60 min applying a Kjeltec™ autoanalyzer 8400 (Foss, Denmark). The fat content was determined gravimetrically by applying the Röse-Gottlieb principle specified in: ‘‘ISO 1211:2010(E)/IDF 1:2010(E) Milk Determination of Fat Content—Röse-Gottlieb Gravimetric Method’’. The lactose content was determined applying a UV-spectroscopic method based on the following standard: ‘‘ISO 5765-2, IDF 79-2:2002—Dried Milk, Dried Ice-Mixes & Processed Cheese— Determination of Lactose Content, Enzymatic Methods’’. The kit applied was a Boehringer Mannheim test-kit, Cat. No. 176 303 (R-biopharm RG, Germany). There were 39 cases with different quality values in total 2 replicates for each resulting in 78 samples for the analysis, however fat content was measured only for 16 cases (32 samples). Furthermore another 30 ml sample was collected from the process and divided into subsamples. The spectra were acquired for each of the subsamples in order to validate the sampling procedure applied. 2.2. Spectra acquisition and preprocessing Preliminary experiments were carried out to find which of the two spectra acquisition techniques—reflectance and transmittance—is the best for analysis of the samples. Finally the transmittance spectra taken using 0.6 ml glass vials (3 mm inner diameter) were found as optimal. The final spectra were acquired using NIR spectrometer InfraQuant (Q-Interline, Tølløse, Denmark) with InAs detector (working
range 4000–10,000 cm 1, spectral resolution was set to 8 cm 1) coupled with a vial holder with temperature controller for keeping temperature of the samples constant. Two spectra were taken for each sample—at 30 and 40 °C. Because of using relatively small vials, the sampling procedure was tested by total fractionation of one sample and analysis of spectra taken from each fraction. The spectra were analyzed visually, statistically, and by means of principal component analysis and the analysis did not reveal any problems. The spectra had very small and random variation. 2.3. Principal component analysis Principal component analysis (PCA) (Esbensen & Geladi, 2009) was employed to explore the relationship among the spectra and to find possible correlation between the spectral variance, temperature and the quality parameters. PCA transforms data from an original variable space (in our case—spectral bands) into a new space formed by principal components. Each principal component is a linear combination of the original variables; the weights of the variables in the combination (loadings) show the influence a particular variable has to the orientation of a principal component. The components are located in the original space so to capture directions of maximum variation of the data points; therefore in most of cases each principal component corresponds to a hidden factor, which causes the variation (e.g., concentration of a particular constituent in a mixing, temperature changes, etc.). By projecting the data points to the principal component space it is possible to decrease the dimension of the data, since the principal components better represent the variation of data values. Besides that, the projections give a new look to the interrelationship of the data points in terms of the hidden factors principal components describe—scores plots. Both scores and loadings plot are very useful to reveal hidden structure (e.g., groups, trends, outliers) in the data and see which variables are particularly important for the structural part. 2.4. Regression, validation and variable selection Partial Least-Squares (PLS) regression (Wold, Sjöström, & Eriksson, 2001) was chosen to build linear regression models for prediction of lactose, total solids, fat and protein content. PLS was invented for the cases where data has a high degree of multicollinearity, therefore it is particularly useful for making regression models on spectroscopic measurements. It transforms both predictors (matrix X) and responses (matrix Y) to a latent variable space in similar way as PCA does but the directions of the latent variables in PLS are chosen so to get maximum covariance between the projected X and Y data points. Besides the main outcome—regression coefficients—PLS also provides many plots for exploring relationship among objects, original and latent variables. Thus outliers detection in the present work was carried out by investigation residual distance plot that show both how good a sample is described by the PLS latent variable space (Q2 residuals) as well as how close the projection of the sample is to the other lot (Hotelling’s T2 values). Random repeated cross-validation with four segments and four repetitions was used for validating the PLS models and chose optimal number of latent variables. Having too many cross-validation segments (e.g., like in full cross-validation) may give wrong evaluation of prediction performance while using too few segments leads to lack of reproducibility. This can be tackled by using repeated cross-validation with several iterations. The models were evaluated and compared using number of latent variables (nLV), determination coefficient (R2), root mean
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
273
Fig. 1. Acquired spectra (top—original, bottom—with noisy regions removed and after SNV transformation).
Fig. 2. Scores and loadings plots for PCA analysis of the preprocessed spectra (light gray curve on loadings plot is a mean spectrum).
squared error of cross-validation (RMSECV), and residual predictive deviation (RPD) which is a ratio between sample standard deviation and standard error of cross-validated predictions (RPD = SDy/SECV) (Williams & Sobering, 1993).
Variable selection methods being applied to a spectral data may improve regression model performance by reducing number of original variables, keeping only the once, most important for prediction of response values. Besides that, using variable selection
274
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
Table 1 Performance of PLS regression models. Parameter
T, °C
Var. sel.
Range
nLV
RMSECV
R2
RPD
TS102 TS102 TS102
30 30 30
No No CARS
Full Reduced Reduced
6 4 4
0.32 0.25 0.19
0.996 0.998 0.999
15.4 21.5 27.4
Prot Prot Prot
30 30 30
No No CARS
Full Reduced Reduced
6 5 4
0.38 0.31 0.21
0.995 0.997 0.999
14.3 17.7 26.4
EzL EzL EzL
30 30 30
No No CARS
Full Reduced Reduced
8 5 4
0.23 0.19 0.14
0.965 0.975 0.986
5.1 6.2 8.4
Fat Fat Fat
30 30 30
No No CARS
Full Reduced Reduced
6 5 4
0.04 0.04 0.03
0.993 0.993 0.996
10.8 11.2 15.1
TS102 TS102 TS102
40 40 40
No No CARS
Full Reduced Reduced
6 4 3
0.38 0.19 0.18
0.995 0.999 0.999
13.8 27.3 29.0
Prot Prot Prot
40 40 40
No No CARS
Full Reduced Reduced
6 6 3
0.36 0.19 0.18
0.996 0.999 0.999
15.0 28.4 30.0
EzL EzL EzL
40 40 40
No No CARS
Full Reduced Reduced
6 5 4
0.32 0.22 0.17
0.936 0.968 0.981
3.8 5.6 7.2
Fat Fat Fat
40 40 40
No No CARS
Full Reduced Reduced
6 3 2
0.03 0.04 0.03
0.994 0.991 0.996
12.9 10.6 15.1
often lead to a better interpretability of the spectra and their relations with responses. In this work the selection and interpretation of the most important variables has been carried out in several steps. First, regression coefficients plot have been observed in order to find out noisy and inert ranges of the spectra. After that, the results were compared with variable importance on projection (VIP) values (Wold et al., 2001). Finally, the competitive adaptive reweighted sampling (CARS) method was used to find best variables or their ranges and possible optimization of the PLS models. CARS improves regression models by selecting variables with large absolute regression coefficients (Li, Liang, Xu, & Cao, 2009). The algorithm uses Monte-Carlo sampling to build N sequential PLS models on randomly selected samples (80% of the calibration set). On each step the absolute regression coefficients are recorded and a variable importance indicator is calculated by normalizing the coefficient values. Then a ratio of variables to be kept is computed and variables with small absolute coefficients are removed. After that, a subset of variables is selected (among the ones being not removed) using adaptive reweighted sampling and RMSECV is computed and saved for the selected subset. The higher absolute regression coefficient a variable has, the more chance to be selected. Finally the algorithm ends up with N subsets of variables and corresponding RMSECV values. A subset with smallest RMSECV is selected as the optimal. Originally developed for dealing with spectra, the CARS algorithm works well for any multivariate data with high collinearity. All calculations were carried out in MATLAB R2012a (Mathworks) supplemented with PLS_Toolbox v. 7.3 (Eigenvector Research Inc., Wenatchee, USA.). 3. Results and discussion 3.1. Exploratory analysis of quality parameters The reference values for the quality parameters were relatively wide spread: from 10% to 30% for total solids, 4.9–24.5% for protein, 1.1–5.1% for lactose and 0.3–1.9% for fat. Analysis of distribu-
tion histograms have shown that for most of the parameters the distribution of the values is close to uniform, however for lactose (EzL) it looks as bimodal, which also can be seen on pairwise scatter plots. This is not unexpected as the samples applied originate from two different production processes thus covering two different product profiles. As global models covering applicable for several different products are strived for, this is to be expected in the future, i.e., the model should be able to handle the difference in lactose. Investigation of pairwise scatter plots and corresponding correlation coefficients shows high correlation between total solids, protein and fat content in the samples (r = 0.96–0.97), which is quite reasonable—total solids are increasing during the filtration process and mostly consist of fat and protein particles. This should make the parameters detectible by spectroscopy since the amount of total solids also influences the optical density of the samples. However it may cause problem with model identification—the regression models built for fat, protein and total solids may all be only sensitive to the amount of total solids in whey and give decent prediction of fat and protein just as a result of the correlation. This can be checked by investigating and comparing regression coefficients and prediction performance of the models. 3.2. Exploratory analysis of spectra The original spectra, transformed from transmittance to absorbance, are shown on the top plot of Fig. 1. Visual inspection of the spectral plots clearly showed a presence of noisy regions caused by absorption of water bands at 4000–4250 and 4800–5350 cm 1 as well as a baseline shift for the rest of the spectra, caused by light scatter. It was decided to remove the noisy regions and apply standard normal variate (SNV) preprocessing to correct the scatter effects (bottom plot of the Fig. 1). Principal component analysis was applied to explore relationship among the spectra, find possible trends and outliers. Fig. 2 shows scores plots (top) and loadings plots (bottom) for PC1–PC2 for the spectra with removed water regions and after SNV transformation (preprocessed spectra). On Fig. 2a all points are colored according to the TS102 values, Fig. 2b shows the same plot but colored according to a temperature of the samples during spectra acquisition. Obviously first principal component shows the amount of total solids whereas second component reveals the temperature effects. The loadings plot for the first component shows that the spectral range 4250–4800 cm 1 may be particularly efficient in prediction of total solids and correlated parameters. The loadings plot for the second component shows mostly one peak around 7100 cm 1 (which is mainly a small water peak combined with CAH combinations). Investigation of the scores plots with higher order components did not show significant correlation with the quality parameters. 3.3. Regression analysis Several PLS regression models were built for each response variable, separately for the spectra taken at 30 and 40 °C. For each temperature value the models were calibrated with original spectra, spectra with removed noisy parts (reduced spectra), as well as using CARS for variable selection. The choice of particular parts for reduced spectra was based on investigation of regression coefficients and VIP scores values for the full range models. The prediction performance of the final models is compared in Table 1. Fig. 3 shows values predicted by the PLS models vs. values measured by ISO chemical tests as well as regression coefficients plots for total solids, protein and fat content models, calibrated using the full range spectra taken at 30 °C. The blue points (denoted in the plot legend by ‘‘Cal’’) show predictions for calibration set and the
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
275
Fig. 3. Predicted vs. measured plots (left) and regression coefficients (right) for PLS regression models for total solids, protein and fat content calibrated on spectra taken at T = 30 °C.
red points (denoted by ‘‘CV’’)—cross-validation predictions. In all three cases the predictions were very good (RPD = 15.4, 14.3 and 10.8 correspondingly). In fact the RMSECV plot has shown a minimum at 8–9 LVs however, based on investigation of PLS-weights and regression coefficients plots, it was decided to use 6 latent variables to avoid possible overfitting. The regression coefficients plots clearly showed a noisy region at 6800–7000 cm 1 (which is another water peak) also the range 7000–10,000 cm 1 had a relatively low importance. This was also proven by calculating and visualizing VIP (variable importance in projection) values (Wold et al., 2001) for the models. Fig. 4 presents the same plots but for models calibrated using spectra reduced to 4250–4800 and 5350–6800 cm 1 range taken at the same temperature. Using the reduced range has improved the performance of the models for total solids and protein predictions (RPD = 21.5 and 17.7) and allowed to use less latent variables in each case. The performance for fat content remained almost unchanged.
CARS variable selection was finally employed in order to further improve the models (RPD = 27.4, 26.4 and 15.1 for total solids, protein and fat content accordingly), see Table 1. The predicted vs. measured and regression coefficients plots for lactose models, calibrated on full and reduced spectra are shown in Fig. 5. The results for lactose content were significantly worse than for the other parameters, which was expected, taking into account the distribution of the lactose values and lack of correlation with any of the principal components. However, by choosing a proper spectral range (4320–4800 and 7220–10,000 cm 1), and applying CARS variable selection it was possible to get model with decent performance (RPD = 8.4). It must also be noticed that the number of outliers in this case was much larger (up to 8) than for the other parameters (from 1 to 3). The models made on the spectra taken at 40 °C showed similar behavior, however performance of the models for total solids and protein content was a bit higher (RPD = 29 and 30 after reducing spectra and applying variable selection algorithm). For the fat
276
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
Fig. 4. Predicted vs. measured plots (left) and regression coefficients (right) for PLS regression models for total solids, protein and fat content calibrated on spectra taken at T = 30 °C reduced to 4250–4800 and 5350–6800 cm 1.
content the results remained the same and lactose was predicted relatively worse in this case.
3.4. Discussion All three main quality parameters—total solids, protein and fat content had cross-validated predictions with very high accuracy and precision, enough to substitute conventional analysis for atline and, possibly, on-line monitoring of the parameters during the filtration process. Apparently longer wavelength, especially the region 4250– 4800 cm 1, which correspond to NAH + CAH, CAH + CAH and CAH + CAC combinations, was the most usable for the PLS regression models. Actually it was possible (not shown in the paper) to get models with almost the same quality just by using this region. Another important region was found to be 5350–6800 cm 1, which correspond to the first CAH overtone.
The correlation between total solids, protein and fat content made it difficult to resolve the corresponding models. It means that all models may in fact predict e.g., total solids and the other parameters are predicted well because of the correlation. However, taking into account that the samples were taken from different productions based on different raw materials, it can be considered as a minor issue. It must also be noticed that the correlation between the parameters is much lower than the correlation between predicted and measured values for each of the parameters and therefore cannot be the only key factor. Investigation of regression coefficients has shown that while having similar trends in the region 5350–6800 cm 1 the coefficients vary significantly in the region 4250–4800 cm 1, in particular around 4500 and 4700 cm 1. The amount of lactose was found as the most challenging parameter to estimate, also due to the quality of the data. However using spectra with reduced wavelength range and variable selection it was possible to get decent predictions that can be used to monitor the lactose content during the production.
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
277
Fig. 5. Predicted vs. measured plots (left) and regression coefficients (right) for PLS regression models for EzL calibrated on spectra taken at T = 30 °C without noisy regions (top) and reduced to 4300–4800 and 7220–1000 cm 1 (bottom).
Temperature expectedly had an effect on the acquired spectra, but both chosen values (30 and 40 °C) gave good models. Taking into account the results for lactose it can be recommended to use 30 °C and of course it must be kept constant during spectra acquisition. 4. Conclusions The results obtained in this feasibility study have shown a high applicability of NIR spectroscopy for monitoring of total solids, protein and fat content in a liquid whey filtration process. The key to success was using a proper spectra acquisition technique (transmittance) and thin sample vials (3 mm in diameter) to make the spectroscopic signal strong enough. The experiments with reflectance and transmittance with bigger vials showed significantly worse results. This makes the method applicable to use on production scale for on-line monitoring of quality parameters during the filtration process. References Barbano, D. M., Lynch, J. M., & Fleming, J. R. (1991). Direct and indirect determination of true protein content of milk by Kjeldahl analysis: Collaborative study. Journal of the Association of Official Analytical Chemists, 74, 281–288. Bogomolov, A., Dietrich, S., Boldrini, B., & Kessler, R. W. (2012). Quantitative determination of fat and total protein in milk based on visible light scatter. Food Chemistry, 134(1), 412–418. http://dx.doi.org/10.1016/j.foodchem.2012.02.077. Bruker Optics (2009). Application note #AN-277E. FT-NIR analysis of liquid milk products. , p. 3.
Crocker, W. P., Jenkins, D. I., Provan, A. L., Macdonald, F. J., Rowland, S. J., & White, J. C. D. (2009). A comparison of the Gerber and Röse-Gottlieb methods for the determination of fat in milk. Journal of Dairy Research, 22(03), 336. http:// dx.doi.org/10.1017/S0022029900007871. De Wit, J. N. (2001). Lecturer’s handbook on whey and whey products. Eindhoven, Netherlands: European Whey Products Association, p. 91. Esbensen, K. H., & Geladi, P. (2009). Principal component analysis: Concept, geometrical interpretation, mathematical background, algorithms, history, practice. In E.-C. S. D. Brown, R. Tauler, & B. Walczak (Eds.). Comprehensive chemometrics (pp. 211–226). Oxford: Elsevier. http://dx.doi.org/10.1016/B978044452701-1.00043-0. Karoui, R., & Debaerdemaeker, J. (2007). A review of the analytical methods coupled with chemometric tools for the determination of the quality and identity of dairy products. Food Chemistry, 102(3), 621–640. http://dx.doi.org/10.1016/ j.foodchem.2006.05.042. Kleyn, D. H. (1985). Determination of lactose by an enzymatic method. Journal of Dairy Science, 68(10), 2791–2798. http://dx.doi.org/10.3168/jds.S00220302(85)81167-X. Li, H., Liang, Y., Xu, Q., & Cao, D. (2009). Key wavelengths screening using competitive adaptive reweighted sampling method for multivariate calibration. Analytica Chimica Acta, 648(1), 77–84. http://dx.doi.org/10.1016/ j.aca.2009.06.046. Pouliot, M., Paquin, P., Martel, R., Gauthier, S. F., & Pouliot, Y. (1997). Whey changes during processing determined by near infrared spectroscopy. Journal of Food Science, 62(3), 475–479. http://dx.doi.org/10.1111/j.1365-2621.1997.tb04409.x. Tsenkova, R., Atanassova, S., Toyoda, K., Ozaki, Y., Itoh, K., & Fearn, T. (1999). Nearinfrared spectroscopy for dairy management: Measurement of unhomogenized milk composition. Journal of Dairy Science, 82(11), 2344–2351. http://dx.doi.org/ 10.3168/jds.S0022-0302(99)75484-6. Williams, P., & Sobering, D. (1993). Comparison of commercial near infrared transmittance and reflectance instruments for analysis of whole grains and seeds. Journal of Near Infrared Spectroscopy, 1(1), 25. http://dx.doi.org/10.1255/ jnirs.3. Wold, S., Sjöström, M., & Eriksson, L. (2001). PLS-regression: A basic tool of chemometrics. Chemometrics and intelligent laboratory . . .. (pp. 109–130) (retrieved from http://www.sciencedirect.com/science/article/pii/ S0169743901001551).
Contents lists available at ScienceDirect
Food Chemistry journal homepage: www.elsevier.com/locate/foodchem
Analytical Methods
Monitoring of whey quality with NIR spectroscopy—A feasibility study Sergey Kucheryavskiy a,⇑, Carina Juel Lomborg b a b
Section of Chemical Engineering, Aalborg University, Campus Esbjerg, Denmark Arla Foods Ingredients Group P/S, Viby, Denmark
a r t i c l e
i n f o
Article history: Received 17 July 2014 Received in revised form 16 December 2014 Accepted 20 December 2014 Available online 30 December 2014 Keywords: Whey Near-infrared spectroscopy Chemometrics
a b s t r a c t The possibility of using near-infrared (NIR) spectroscopy for monitoring of liquid whey quality parameters during protein production process has been tested. The parameters included total solids, lactose, protein and fat content. The samples for the experiment were taken from real industrial processes and had a large variability for most of the parameters. Partial Least Squares (PLS) regression was used to make the prediction models based on NIR spectra taken at 30 and 40 °C. Using proper wavelength range allowed to get models for prediction of fat, protein and amount of total solids with very high precision and accuracy. The lactose was found to be the most challenging parameter. Ó 2014 Elsevier Ltd. All rights reserved.
1. Introduction Whey is a liquid product remaining after production of cheese or removal of fat and casein from cow milk. Being considered as a waste product for many years, nowadays it is widely used in the food industry due to its nutritional properties—whey contains about 50% of the nutrients present in milk, including milk sugar (lactose), serum proteins (whey proteins), minerals, a small amount of fat, and most of the water soluble minor nutrients from milk such as vitamins (de Wit, 2001). Production of protein is one of the applications whey is widely used for. There are several techniques to get the whey protein concentrate; one of which is a membrane filtration, where the whey is sieved through a set of filters of controlled pore size (de Wit, 2001). The concentration of proteins in liquid whey at different stages of the process is changing from 3–5% to 25–30% whereas the lactose content is decreasing. The monitoring of the different constituents during the production is first and foremost important from a quality perspective i.e., being able to control the process so that a uniform product, with respect to the protein and dry matter, is obtained from batch to batch is a significant quality parameter. Furthermore it is important financially, that the product contains the amount of protein paid for by the customer. Conventionally the different constituents may be measured applying standard chemical methods, e.g., Kjeldahl analysis (Barbano, Lynch, & Fleming, 1991) for protein, enzymatic kit ⇑ Corresponding author at: Niels Bohrs vej 8, Esbjerg DK 6700, Denmark. Tel.: +45 2787 9818. E-mail address: [email protected] (S. Kucheryavskiy). http://dx.doi.org/10.1016/j.foodchem.2014.12.086 0308-8146/Ó 2014 Elsevier Ltd. All rights reserved.
(Kleyn, 1985) or HPLC analysis for lactose, and Röse-Gottlieb (Crocker et al., 2009) for fat. However all these methods are laborious and time consuming, and thus not really a possibility for process monitoring. In an attempt to get a quicker result applicable for process monitoring and control, the degrees of brix (°Bx) based on the measurement of the refractive index of a sample from the process stream has been implemented in the production. The refractive brix was originally created to determine the amount of sucrose in water, but is now widely being used across industries for determination of dry matter in process streams. In dealing with raw materials that are not standardized and know to fluctuate over the year, the refractive brix however has it limitations for real process monitoring and understanding. As in applying the method only a measure of the total dry matter is obtained and hence not the composition with respect to lactose and protein, thus meaning a shift in the ratio between lactose and protein may not be caught by this method. As a consequence, during the last decade, analysis based on using wet chemistry tends to be replaced by modern spectroscopic methods combined with multivariate data analysis. Thus, midinfrared spectroscopy has been widely accepted as a laboratory standard for the milk nutrient analysis (Karoui & Debaerdemaeker, 2007). Numerous research has been carried out on using NIR spectroscopy for estimating quality parameters of milk and milk products (Bogomolov, Dietrich, Boldrini, & Kessler, 2012; Karoui & Debaerdemaeker, 2007; Tsenkova et al., 1999). However, to our knowledge, very few papers on estimation of liquid whey quality parameters with NIR spectroscopy exist—most of the research in this field deals with whey powders. In (Pouliot,
272
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
Paquin, Martel, Gauthier, & Pouliot, 1997) authors tried to solve the similar problem, however the data quality as well as obtained results are outside the industrial standards. There are also several commercial analyzers based in FT-NIR (see for example, (Bruker Optics, 2009)) that can analyze liquid whey quality, however they work with a narrow range of the parameters, e.g., use different models and settings for low, middle or high amount of protein content and total solids. The main purpose of the present work was to carry out a feasibility study of using NIR spectroscopy for monitoring of major whey quality parameters: total solids, protein, lactose and fat content during whey filtration processes, where the parameters may vary significantly. 2. Materials and methods 2.1. Samples and measurement of quality parameters Samples for the experiments were collected from four different productions in a real industrial process on whey filtering at Danmark Protein factory (Arla Foods Ingredients Group A/S, Videbæk, Denmark). Two 100 ml replicates were taken for each case and transported to a spectroscopic laboratory in plastic vials. Each sample was also delivered to a factory lab where reference values for total protein (Prot), fat (Fat) and lactose (EzL) content as well as amount of total solids (TS102), were evaluated using the following procedures. The total solid content was determined according to ‘‘ISO 6731:2010(E)/IDF 21:2010(E) Milk, cream and evaporated milk. Determination of the total solids content.’’ The protein content was defined as the nitrogen content times a Kjeldahl factor (6.38), and reported as g 100 g 1 sample. The Nitrogen content was determined according to the following standard: ‘‘ISO 8968-3:2004(E)/IDF 20-3:2004 Determination of nitrogen content. Block-digestion method (semi-micro rapid routine method)’’. In the destruction 0.6–0.8 g of sample was weighted and digested at 420 °C for 60 min applying a Kjeltec™ autoanalyzer 8400 (Foss, Denmark). The fat content was determined gravimetrically by applying the Röse-Gottlieb principle specified in: ‘‘ISO 1211:2010(E)/IDF 1:2010(E) Milk Determination of Fat Content—Röse-Gottlieb Gravimetric Method’’. The lactose content was determined applying a UV-spectroscopic method based on the following standard: ‘‘ISO 5765-2, IDF 79-2:2002—Dried Milk, Dried Ice-Mixes & Processed Cheese— Determination of Lactose Content, Enzymatic Methods’’. The kit applied was a Boehringer Mannheim test-kit, Cat. No. 176 303 (R-biopharm RG, Germany). There were 39 cases with different quality values in total 2 replicates for each resulting in 78 samples for the analysis, however fat content was measured only for 16 cases (32 samples). Furthermore another 30 ml sample was collected from the process and divided into subsamples. The spectra were acquired for each of the subsamples in order to validate the sampling procedure applied. 2.2. Spectra acquisition and preprocessing Preliminary experiments were carried out to find which of the two spectra acquisition techniques—reflectance and transmittance—is the best for analysis of the samples. Finally the transmittance spectra taken using 0.6 ml glass vials (3 mm inner diameter) were found as optimal. The final spectra were acquired using NIR spectrometer InfraQuant (Q-Interline, Tølløse, Denmark) with InAs detector (working
range 4000–10,000 cm 1, spectral resolution was set to 8 cm 1) coupled with a vial holder with temperature controller for keeping temperature of the samples constant. Two spectra were taken for each sample—at 30 and 40 °C. Because of using relatively small vials, the sampling procedure was tested by total fractionation of one sample and analysis of spectra taken from each fraction. The spectra were analyzed visually, statistically, and by means of principal component analysis and the analysis did not reveal any problems. The spectra had very small and random variation. 2.3. Principal component analysis Principal component analysis (PCA) (Esbensen & Geladi, 2009) was employed to explore the relationship among the spectra and to find possible correlation between the spectral variance, temperature and the quality parameters. PCA transforms data from an original variable space (in our case—spectral bands) into a new space formed by principal components. Each principal component is a linear combination of the original variables; the weights of the variables in the combination (loadings) show the influence a particular variable has to the orientation of a principal component. The components are located in the original space so to capture directions of maximum variation of the data points; therefore in most of cases each principal component corresponds to a hidden factor, which causes the variation (e.g., concentration of a particular constituent in a mixing, temperature changes, etc.). By projecting the data points to the principal component space it is possible to decrease the dimension of the data, since the principal components better represent the variation of data values. Besides that, the projections give a new look to the interrelationship of the data points in terms of the hidden factors principal components describe—scores plots. Both scores and loadings plot are very useful to reveal hidden structure (e.g., groups, trends, outliers) in the data and see which variables are particularly important for the structural part. 2.4. Regression, validation and variable selection Partial Least-Squares (PLS) regression (Wold, Sjöström, & Eriksson, 2001) was chosen to build linear regression models for prediction of lactose, total solids, fat and protein content. PLS was invented for the cases where data has a high degree of multicollinearity, therefore it is particularly useful for making regression models on spectroscopic measurements. It transforms both predictors (matrix X) and responses (matrix Y) to a latent variable space in similar way as PCA does but the directions of the latent variables in PLS are chosen so to get maximum covariance between the projected X and Y data points. Besides the main outcome—regression coefficients—PLS also provides many plots for exploring relationship among objects, original and latent variables. Thus outliers detection in the present work was carried out by investigation residual distance plot that show both how good a sample is described by the PLS latent variable space (Q2 residuals) as well as how close the projection of the sample is to the other lot (Hotelling’s T2 values). Random repeated cross-validation with four segments and four repetitions was used for validating the PLS models and chose optimal number of latent variables. Having too many cross-validation segments (e.g., like in full cross-validation) may give wrong evaluation of prediction performance while using too few segments leads to lack of reproducibility. This can be tackled by using repeated cross-validation with several iterations. The models were evaluated and compared using number of latent variables (nLV), determination coefficient (R2), root mean
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
273
Fig. 1. Acquired spectra (top—original, bottom—with noisy regions removed and after SNV transformation).
Fig. 2. Scores and loadings plots for PCA analysis of the preprocessed spectra (light gray curve on loadings plot is a mean spectrum).
squared error of cross-validation (RMSECV), and residual predictive deviation (RPD) which is a ratio between sample standard deviation and standard error of cross-validated predictions (RPD = SDy/SECV) (Williams & Sobering, 1993).
Variable selection methods being applied to a spectral data may improve regression model performance by reducing number of original variables, keeping only the once, most important for prediction of response values. Besides that, using variable selection
274
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
Table 1 Performance of PLS regression models. Parameter
T, °C
Var. sel.
Range
nLV
RMSECV
R2
RPD
TS102 TS102 TS102
30 30 30
No No CARS
Full Reduced Reduced
6 4 4
0.32 0.25 0.19
0.996 0.998 0.999
15.4 21.5 27.4
Prot Prot Prot
30 30 30
No No CARS
Full Reduced Reduced
6 5 4
0.38 0.31 0.21
0.995 0.997 0.999
14.3 17.7 26.4
EzL EzL EzL
30 30 30
No No CARS
Full Reduced Reduced
8 5 4
0.23 0.19 0.14
0.965 0.975 0.986
5.1 6.2 8.4
Fat Fat Fat
30 30 30
No No CARS
Full Reduced Reduced
6 5 4
0.04 0.04 0.03
0.993 0.993 0.996
10.8 11.2 15.1
TS102 TS102 TS102
40 40 40
No No CARS
Full Reduced Reduced
6 4 3
0.38 0.19 0.18
0.995 0.999 0.999
13.8 27.3 29.0
Prot Prot Prot
40 40 40
No No CARS
Full Reduced Reduced
6 6 3
0.36 0.19 0.18
0.996 0.999 0.999
15.0 28.4 30.0
EzL EzL EzL
40 40 40
No No CARS
Full Reduced Reduced
6 5 4
0.32 0.22 0.17
0.936 0.968 0.981
3.8 5.6 7.2
Fat Fat Fat
40 40 40
No No CARS
Full Reduced Reduced
6 3 2
0.03 0.04 0.03
0.994 0.991 0.996
12.9 10.6 15.1
often lead to a better interpretability of the spectra and their relations with responses. In this work the selection and interpretation of the most important variables has been carried out in several steps. First, regression coefficients plot have been observed in order to find out noisy and inert ranges of the spectra. After that, the results were compared with variable importance on projection (VIP) values (Wold et al., 2001). Finally, the competitive adaptive reweighted sampling (CARS) method was used to find best variables or their ranges and possible optimization of the PLS models. CARS improves regression models by selecting variables with large absolute regression coefficients (Li, Liang, Xu, & Cao, 2009). The algorithm uses Monte-Carlo sampling to build N sequential PLS models on randomly selected samples (80% of the calibration set). On each step the absolute regression coefficients are recorded and a variable importance indicator is calculated by normalizing the coefficient values. Then a ratio of variables to be kept is computed and variables with small absolute coefficients are removed. After that, a subset of variables is selected (among the ones being not removed) using adaptive reweighted sampling and RMSECV is computed and saved for the selected subset. The higher absolute regression coefficient a variable has, the more chance to be selected. Finally the algorithm ends up with N subsets of variables and corresponding RMSECV values. A subset with smallest RMSECV is selected as the optimal. Originally developed for dealing with spectra, the CARS algorithm works well for any multivariate data with high collinearity. All calculations were carried out in MATLAB R2012a (Mathworks) supplemented with PLS_Toolbox v. 7.3 (Eigenvector Research Inc., Wenatchee, USA.). 3. Results and discussion 3.1. Exploratory analysis of quality parameters The reference values for the quality parameters were relatively wide spread: from 10% to 30% for total solids, 4.9–24.5% for protein, 1.1–5.1% for lactose and 0.3–1.9% for fat. Analysis of distribu-
tion histograms have shown that for most of the parameters the distribution of the values is close to uniform, however for lactose (EzL) it looks as bimodal, which also can be seen on pairwise scatter plots. This is not unexpected as the samples applied originate from two different production processes thus covering two different product profiles. As global models covering applicable for several different products are strived for, this is to be expected in the future, i.e., the model should be able to handle the difference in lactose. Investigation of pairwise scatter plots and corresponding correlation coefficients shows high correlation between total solids, protein and fat content in the samples (r = 0.96–0.97), which is quite reasonable—total solids are increasing during the filtration process and mostly consist of fat and protein particles. This should make the parameters detectible by spectroscopy since the amount of total solids also influences the optical density of the samples. However it may cause problem with model identification—the regression models built for fat, protein and total solids may all be only sensitive to the amount of total solids in whey and give decent prediction of fat and protein just as a result of the correlation. This can be checked by investigating and comparing regression coefficients and prediction performance of the models. 3.2. Exploratory analysis of spectra The original spectra, transformed from transmittance to absorbance, are shown on the top plot of Fig. 1. Visual inspection of the spectral plots clearly showed a presence of noisy regions caused by absorption of water bands at 4000–4250 and 4800–5350 cm 1 as well as a baseline shift for the rest of the spectra, caused by light scatter. It was decided to remove the noisy regions and apply standard normal variate (SNV) preprocessing to correct the scatter effects (bottom plot of the Fig. 1). Principal component analysis was applied to explore relationship among the spectra, find possible trends and outliers. Fig. 2 shows scores plots (top) and loadings plots (bottom) for PC1–PC2 for the spectra with removed water regions and after SNV transformation (preprocessed spectra). On Fig. 2a all points are colored according to the TS102 values, Fig. 2b shows the same plot but colored according to a temperature of the samples during spectra acquisition. Obviously first principal component shows the amount of total solids whereas second component reveals the temperature effects. The loadings plot for the first component shows that the spectral range 4250–4800 cm 1 may be particularly efficient in prediction of total solids and correlated parameters. The loadings plot for the second component shows mostly one peak around 7100 cm 1 (which is mainly a small water peak combined with CAH combinations). Investigation of the scores plots with higher order components did not show significant correlation with the quality parameters. 3.3. Regression analysis Several PLS regression models were built for each response variable, separately for the spectra taken at 30 and 40 °C. For each temperature value the models were calibrated with original spectra, spectra with removed noisy parts (reduced spectra), as well as using CARS for variable selection. The choice of particular parts for reduced spectra was based on investigation of regression coefficients and VIP scores values for the full range models. The prediction performance of the final models is compared in Table 1. Fig. 3 shows values predicted by the PLS models vs. values measured by ISO chemical tests as well as regression coefficients plots for total solids, protein and fat content models, calibrated using the full range spectra taken at 30 °C. The blue points (denoted in the plot legend by ‘‘Cal’’) show predictions for calibration set and the
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
275
Fig. 3. Predicted vs. measured plots (left) and regression coefficients (right) for PLS regression models for total solids, protein and fat content calibrated on spectra taken at T = 30 °C.
red points (denoted by ‘‘CV’’)—cross-validation predictions. In all three cases the predictions were very good (RPD = 15.4, 14.3 and 10.8 correspondingly). In fact the RMSECV plot has shown a minimum at 8–9 LVs however, based on investigation of PLS-weights and regression coefficients plots, it was decided to use 6 latent variables to avoid possible overfitting. The regression coefficients plots clearly showed a noisy region at 6800–7000 cm 1 (which is another water peak) also the range 7000–10,000 cm 1 had a relatively low importance. This was also proven by calculating and visualizing VIP (variable importance in projection) values (Wold et al., 2001) for the models. Fig. 4 presents the same plots but for models calibrated using spectra reduced to 4250–4800 and 5350–6800 cm 1 range taken at the same temperature. Using the reduced range has improved the performance of the models for total solids and protein predictions (RPD = 21.5 and 17.7) and allowed to use less latent variables in each case. The performance for fat content remained almost unchanged.
CARS variable selection was finally employed in order to further improve the models (RPD = 27.4, 26.4 and 15.1 for total solids, protein and fat content accordingly), see Table 1. The predicted vs. measured and regression coefficients plots for lactose models, calibrated on full and reduced spectra are shown in Fig. 5. The results for lactose content were significantly worse than for the other parameters, which was expected, taking into account the distribution of the lactose values and lack of correlation with any of the principal components. However, by choosing a proper spectral range (4320–4800 and 7220–10,000 cm 1), and applying CARS variable selection it was possible to get model with decent performance (RPD = 8.4). It must also be noticed that the number of outliers in this case was much larger (up to 8) than for the other parameters (from 1 to 3). The models made on the spectra taken at 40 °C showed similar behavior, however performance of the models for total solids and protein content was a bit higher (RPD = 29 and 30 after reducing spectra and applying variable selection algorithm). For the fat
276
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
Fig. 4. Predicted vs. measured plots (left) and regression coefficients (right) for PLS regression models for total solids, protein and fat content calibrated on spectra taken at T = 30 °C reduced to 4250–4800 and 5350–6800 cm 1.
content the results remained the same and lactose was predicted relatively worse in this case.
3.4. Discussion All three main quality parameters—total solids, protein and fat content had cross-validated predictions with very high accuracy and precision, enough to substitute conventional analysis for atline and, possibly, on-line monitoring of the parameters during the filtration process. Apparently longer wavelength, especially the region 4250– 4800 cm 1, which correspond to NAH + CAH, CAH + CAH and CAH + CAC combinations, was the most usable for the PLS regression models. Actually it was possible (not shown in the paper) to get models with almost the same quality just by using this region. Another important region was found to be 5350–6800 cm 1, which correspond to the first CAH overtone.
The correlation between total solids, protein and fat content made it difficult to resolve the corresponding models. It means that all models may in fact predict e.g., total solids and the other parameters are predicted well because of the correlation. However, taking into account that the samples were taken from different productions based on different raw materials, it can be considered as a minor issue. It must also be noticed that the correlation between the parameters is much lower than the correlation between predicted and measured values for each of the parameters and therefore cannot be the only key factor. Investigation of regression coefficients has shown that while having similar trends in the region 5350–6800 cm 1 the coefficients vary significantly in the region 4250–4800 cm 1, in particular around 4500 and 4700 cm 1. The amount of lactose was found as the most challenging parameter to estimate, also due to the quality of the data. However using spectra with reduced wavelength range and variable selection it was possible to get decent predictions that can be used to monitor the lactose content during the production.
S. Kucheryavskiy, C.J. Lomborg / Food Chemistry 176 (2015) 271–277
277
Fig. 5. Predicted vs. measured plots (left) and regression coefficients (right) for PLS regression models for EzL calibrated on spectra taken at T = 30 °C without noisy regions (top) and reduced to 4300–4800 and 7220–1000 cm 1 (bottom).
Temperature expectedly had an effect on the acquired spectra, but both chosen values (30 and 40 °C) gave good models. Taking into account the results for lactose it can be recommended to use 30 °C and of course it must be kept constant during spectra acquisition. 4. Conclusions The results obtained in this feasibility study have shown a high applicability of NIR spectroscopy for monitoring of total solids, protein and fat content in a liquid whey filtration process. The key to success was using a proper spectra acquisition technique (transmittance) and thin sample vials (3 mm in diameter) to make the spectroscopic signal strong enough. The experiments with reflectance and transmittance with bigger vials showed significantly worse results. This makes the method applicable to use on production scale for on-line monitoring of quality parameters during the filtration process. References Barbano, D. M., Lynch, J. M., & Fleming, J. R. (1991). Direct and indirect determination of true protein content of milk by Kjeldahl analysis: Collaborative study. Journal of the Association of Official Analytical Chemists, 74, 281–288. Bogomolov, A., Dietrich, S., Boldrini, B., & Kessler, R. W. (2012). Quantitative determination of fat and total protein in milk based on visible light scatter. Food Chemistry, 134(1), 412–418. http://dx.doi.org/10.1016/j.foodchem.2012.02.077. Bruker Optics (2009). Application note #AN-277E. FT-NIR analysis of liquid milk products. , p. 3.
Crocker, W. P., Jenkins, D. I., Provan, A. L., Macdonald, F. J., Rowland, S. J., & White, J. C. D. (2009). A comparison of the Gerber and Röse-Gottlieb methods for the determination of fat in milk. Journal of Dairy Research, 22(03), 336. http:// dx.doi.org/10.1017/S0022029900007871. De Wit, J. N. (2001). Lecturer’s handbook on whey and whey products. Eindhoven, Netherlands: European Whey Products Association, p. 91. Esbensen, K. H., & Geladi, P. (2009). Principal component analysis: Concept, geometrical interpretation, mathematical background, algorithms, history, practice. In E.-C. S. D. Brown, R. Tauler, & B. Walczak (Eds.). Comprehensive chemometrics (pp. 211–226). Oxford: Elsevier. http://dx.doi.org/10.1016/B978044452701-1.00043-0. Karoui, R., & Debaerdemaeker, J. (2007). A review of the analytical methods coupled with chemometric tools for the determination of the quality and identity of dairy products. Food Chemistry, 102(3), 621–640. http://dx.doi.org/10.1016/ j.foodchem.2006.05.042. Kleyn, D. H. (1985). Determination of lactose by an enzymatic method. Journal of Dairy Science, 68(10), 2791–2798. http://dx.doi.org/10.3168/jds.S00220302(85)81167-X. Li, H., Liang, Y., Xu, Q., & Cao, D. (2009). Key wavelengths screening using competitive adaptive reweighted sampling method for multivariate calibration. Analytica Chimica Acta, 648(1), 77–84. http://dx.doi.org/10.1016/ j.aca.2009.06.046. Pouliot, M., Paquin, P., Martel, R., Gauthier, S. F., & Pouliot, Y. (1997). Whey changes during processing determined by near infrared spectroscopy. Journal of Food Science, 62(3), 475–479. http://dx.doi.org/10.1111/j.1365-2621.1997.tb04409.x. Tsenkova, R., Atanassova, S., Toyoda, K., Ozaki, Y., Itoh, K., & Fearn, T. (1999). Nearinfrared spectroscopy for dairy management: Measurement of unhomogenized milk composition. Journal of Dairy Science, 82(11), 2344–2351. http://dx.doi.org/ 10.3168/jds.S0022-0302(99)75484-6. Williams, P., & Sobering, D. (1993). Comparison of commercial near infrared transmittance and reflectance instruments for analysis of whole grains and seeds. Journal of Near Infrared Spectroscopy, 1(1), 25. http://dx.doi.org/10.1255/ jnirs.3. Wold, S., Sjöström, M., & Eriksson, L. (2001). PLS-regression: A basic tool of chemometrics. Chemometrics and intelligent laboratory . . .. (pp. 109–130) (retrieved from http://www.sciencedirect.com/science/article/pii/ S0169743901001551).
