Prediction of Benzoyl Peroxide in Flour Using Near Infrared Spectroscopy Technique

Abstract

This research proposed to design a prediction model based on radial basis function neural network and near infrared reflectance spectroscopy in detecting concentration of benzoyl peroxide in flour. Near infrared reflectance spectra acquired from 100 different concentration samples were pre-processed by the standard normal variate method, detection of leverage, and student residual. Near infrared reflectance spectroscopy models were designed to predict benzoyl peroxide in the 36 samples by means of partial least squares, back propagation neural network, and radial basis function, respectively. The results demonstrated that the radial basis function model, with prediction correlation coefficient (R), root mean squared error of prediction, and ratio of performance to standard deviate reaching 0.9937, 15.5095, and 8.8216, respectively, had optimal prediction accuracy and feasibility providing quality evaluation and dynamic monitoring service for quality inspection department and consumers.

Keywords:

• Benzoyl peroxide
• Near infrared reflectance spectroscopy
• Radial basis function

INTRODUCTION

Flour can be used to make bread, noodles, steamed buns, and other flour products, and its safety and quality are still considered as important issues all over the world. Usually, most consumers prefer to buy the flour and its by-products that looked whiter requiring a higher whiteness of flour, which forced the wheat processing enterprises to attach great importance to flour whiteness.^[1,^2] One of the components of wheat is lutein, it turns the flour yellow for a long time after it is ground. The conjugate double-bond in lutein molecules is a chromophore, it will be oxidized slowly in the air by means of absorbing light and cause the flour to fade,^[3] but a long period of storage, at least several months, is required, so it is not feasible for the sake of business.^[4] Benzoyl peroxide (BPO) which will promote flour to turn white due to its feature of strong oxidation^[5] is one of flour brighteners a lot of flour processing enterprises widely use to cater to the demand of consumers.^[6] But the excessive use of BPO will destroy the carotene, vitamin A, vitamin E, vitamin B1, and vitamin B2 in flour.^[7] If people eat the flour with excessive BPO for a long time, it will cause vitamin deficiencies and some diseases, such as angular stomatitis, neuritis, the cumulative damage of central nervous system, and even failure of the liver,^[8] and because of which many countries, including China, have issued corresponding laws and regulations to restrain the use of BPO. For example, Chinese GB2760-2007 stipulates the concentration of BPO in flour should not exceed 60 mg·Kg^–1; the maximum concentration of BPO in flour is 80 mg·Kg^–1 in the United States, which is 300 mg·Kg^–1 in Japan, and 50 mg·Kg^–1 in the United Kingdom.^[9] In March 2011, the Ministry of Public Health of China issued the latest announcement that it was forbidden to add brightener into flour From May 1, 2011. But some small flour enterprises are still using BPO as a flour brightener^[10] due to the fierce market competition and the low cost of BPO. So, it is of vital importance to detect BPO in flour rapidly by technical means, which can not only provide powerful technical support for sampling observation of quality supervision department and standardize the market order, but also ensure the safety of flour quality and the consumers’ rights and interests. Currently, China has established two national standards, GB/T18415-2001 (gas chromatography) and GB/T22325-2008 (high-performance liquid chromatography [HPLC]), to detect BPO besides other detection methods involving electrochemical analysis,^[11] Chemiluminescence,^[12] and iodometry,^[13] etc. But these methods are too complex, time consuming, and chemical reagents consuming, and they also require professionals with a variety of specialized equipments, which are problems commonly existing in the process of quantitative and qualitative analysis of other crops and their by-products. Recently, many researchers are seeking a more speedy and simple method to achieve the quantitative and qualitative analysis of the food addictive in raw material to realize the real-time monitoring quality.

Near infrared reflectance spectroscopy (NIRS) is a green and non-destructive testing technology,^[14] which has become a very effective method of speedy detection to take place of the chemical test. Ni et al.^[15] predicted the fat content, water content, acid value, and peroxide value in potatoes using the least squares support vector machine (LS-SVM) algorithm together with NIRS, and the predictive correlation coefficients of each ingredient were all above 0.949, with the maximum of the root mean squared error of prediction (RMSEP) being only 0.435, the minimum value being 0.076. Luo^[16] predicted the polyphenol and the polysaccharide in tea using the NIRS technology in combination with artificial neural network models with R, RMSEP, RSEP being 0.9804, 0.529, 0.017, and 0.9682, 0.111, 0.030, respectively. Latorre et al.^[17] designed an origin recognition model based on the NIRS and soft independent modeling of class analogy (SIMCA) analysis with the sensitivity being 93.3% and the specificity being 100%. Zhang et al.^[18] adopted the LS-SVM model based on NIRS analysis to identify the nuances scar of fresh jujube fruit, and the identification accuracy could reach 100%. Ferreira et al.^[19] used NIRS, the partial least squares (PLS) algorithm, and multivariate calibration to predict the moisture, ash, protein, and fat in soybeans, with the minimum value of the decision coefficient being 0.72, the maximum value being 0.88 and root mean square error of cross validation (RMSECV) all being below 2.09%. Zhang et al.^[20] designed a model between the spectral data and BPO using PLS algorithm, with the coefficient of calibration set and prediction set being 0.8901 and 0.8865, respectively; the RMSE being 40.85 mg·Kg^–1 and 44.69 mg·Kg^–1, respectively. Haughey et al.^[21] detected the melamine in soybean meal using NIRS in combination with PCA, with R² reaching 0.89–0.99 and the RMSE of the calibration set and prediction set being 0.081–0.276 and 0.134–0.368%, respectively. Mouazen et al.^[22] detected the glucose in a mixture of Saudi and imported honey samples adulterated with glucose syrup of five concentrations: 0, 5, 12, 19, and 33%, using visible/near infrared reflectance (VIS/NIR), with the standard deviation divided by RMSEP being 4.52, 0.85, and 2.53 g/100 g, respectively. On the basis of predecessors’ researches, this research proposed to design a prediction model for detecting the flour of different BPO concentrations using NIRS in combination with the radial basis function (RBF) algorithm. The feasibility and validity of this method is validated through a great deal of experiments.

MATERIALS AND METHODS

The Source and Preparation of the Samples

The flour in the market may have been added a variety of additives, such as gluten fortifier and brightener, in order to ensure the accuracy of experimental results, the original flour used in this experiment was hard-red-spring wheat in the North-Eastern area which was gathered and ground by the Agricultural Products Quality and Safety Research Institute of Heilongjiang Academy of Agricultural Sciences. The analytical reagent of BPO produced by Sigma in the United States was preserved in water with number 517909 and 75% concentration. The stepwise dilution method was adopted to prepare the samples instead of preparing directly to decrease experiment error, since the required concentration of BPO was lower. The procedures were as follows: First, make a 1% concentration sample (10,000 mg kg^–1); then, add flour into the 1% concentration sample to dilute it to 0.1% (1000 mg kg^–1); finally, take the 0.1% concentration sample as reference sample to make up the required concentration samples. One hundred different concentration samples were prepared, ranging from 10 to 505 mg kg^–1, with the concentration gradient being 5 mg kg^–1, the mean value being 257.5mg kg^–1, and the standard deviation being 145.06 mg kg^–1. The prepared samples were put in the same environment (room temperature 25°C, dark, and ventilation) for 48 h to make sure that the BPO reacted well to lutein.

Collection of Spectral Data

NIRS DS2500 analyzer with a scanning range of 400–2498 nm by increments 2 nm produced by FOSS company was adopted to scan the spectra at the temperature between 21–25°C. While loading the sample, the rim of cup was smoothed after a sample was packed into the sample cup to make it even with the cup’s edge in order to avoid the objective interference. To eliminate the influence of inconsistent uniformity, each sample was measured 10 times repeatedly, and the average of 10 times would be taken. One hundred groups of spectral data were obtained with each group containing 1050 spectral points.

Measurement of BPO by Physical and Chemical Methods

When the BPO was put into flour, it would react to the moisture in flour and decompose into benzoic acid under the action of enzyme. In theory, the full reaction of the BPO was required, but it was hard to control the degree of the reaction in fact. So before measuring the benzoic acid in samples, the unreacted BPO should be extracted and reduced to benzoic acid. In order to examine the accuracy of the prepared concentrations, the HPLC method was adopted to measure the actual concentrations of BPO of some samples which would be compared with the prepared concentrations. At first, 11 prepared samples were selected with their concentrations being 10, 50, 100, 150, 200, 250, 300, 350, 400, 450, 500 mg/kg, respectively. Methanol was adopted to extract BPO, and potassium iodide was used as a reducing agent to restore the BPO to benzoic acid. The experimental procedures were according to the Chinese GB/T22325-2008^[23] as follows: Measure 5 g (the accuracy is 0.1 mg) of the prepared samples and place them in colorimetric tubes of 50.0 mL; add 10.0mL of methanol; blend them in vortex mixer for 1 min; lay them still for 5 min, and add 5.0 mL of 50% potassium iodide solution; blend them in vortex mixer for 1 min; lay them still for 10 min; add water up to 50.0 mL; blend, lay still and absorb supernatant fluid through the 0.22 μm membrane filters with filtrate, place them in sample bottles for use. And then make nine copies of benzoic acid standard solutions, the concentrations of them were 0, 5.0, 10.0, 25.0, 50.0, 100.0, 200.0, 400.0, 600.0 μg/mL, respectively, with the benzoic acid as solute and the methanol as solvent. Next, measure nine copies of 5 g pure flour, place them into nine colorimetric tubes of 50 mL and put 10.00 mL of prepared benzoic acid standard solutions of different concentrations into those nine colorimetric tubes, respectively. Deal with those nine colorimetric tubes in the vortex mixer as above, and the final concentrations of the standard solutions were 0, 1.0, 2.0, 5.0, 10.0, 20.0, 40.0, 80.0, and 120.0 μg/mL, respectively. Inject 10.0 μL of benzoic acid standard solutions of different concentrations into the liquid chromatograph (Agilent 1100 series). The parameters of liquid chromatograph were set as follows: the chromatographic column was 4.6 mm × 250 mm, the test wavelength was 230 nm, the proportion between methanol and water (including 0.02 mol/L ammonium acetate) in the mobile phase was 10:90 (volume fraction), the flow velocity was 1.0 mL/min and the injection volume was 10.0 μL. Figure 1 showed the liquid chromatogram of benzoic acid standard solution of 5.0 μg/mL, it could be seen a better peak shape when the elution time was 4.805 min. Figure 2 showed the standard curve with the benzoic acid concentration being the abscissa and the peak area of benzoic acid being the ordinate. Finally, inject 10.0 μL of the sample solutions mentioned above into liquid chromatograph, the parameters of liquid chromatograph were set as above, got the liquid chromatogram of samples, look up the corresponding concentration of the benzoic acid from the curve according to the peak area of benzoic acid, and calculate the concentration of BPO in the samples by Formula 2.1.

FIGURE 1 The liquid chromatogram of benzoic acid standard solution of 5.0 μg/mL.

(2.1)

FIGURE 2 The standard curve of benzoic acid standard solution.

where D is the concentration of BPO in a sample, unit: g/kg, c is the concentration of benzoic acid detected by standard curve, unit: μg/mL, V is the volume of a sample extract, unit: mL, m is the quality of a sample, unit: g. Constant 0.992 is the conversion factor for converting between benzoic acid and BPO: 242.2/(2 × 122.1). Two significant figures would be kept in results.

The measurements of 11 samples were shown in Table 1. By analyzing the data in Table 1, the recoveries of 11 samples were over than 95% except for the 10 mg/kg concentration sample whose recovery was 90% and the average recovery of 11 samples was 98.1%, which illustrated that the prepared concentrations of samples had high accuracy. In conclusion, the concentrations of BPO prepared in this study could be used as the standard concentrations.

TABLE 1 The comparison between configured concentration and measured concentration of BPO

Configured concentration (mg/kg)	10	50	100	150	200	250	300	350	400	450	500
Measured concentration (mg/kg)	9	49	99	149	198	247	297	346	396	446	495
Recovery (%)	90.0	98.0	99.0	99.3	99.0	98.8	99.0	98.9	99.0	99.1	99.0

Calibration Methods and Evaluation of the Model

T detection of leverage and student residual

It is of vital importance to eliminate abnormal samples, since the measuring equipments, test methods, personal operation, or source of samples might lead to the abnormal samples, which could cause maladjustment of models. In general, abnormal samples were composed of two parts: concentration anomaly and spectral anomaly. If the abnormal samples were caused by improper operation or equipment abnormality, the ones should be scanned again when discovered, while they must be eliminated directly if the abnormal samples were caused by themselves. The methods commonly used for eliminating abnormal samples are residual error of predicting concentration,^[24] Mahalanobis Distance,^[25] T detection of leverage and student residual.^[26] In this study, the T detection of leverage and student residual were adopted because their ability of verifying whether concentrations are abnormal or not and the contributions to model are given by samples. All samples were analyzed according to the distribution chart of residuals distribution of the leverage values and the student residuals.

Suppose there was a spectral matrix, where n denoted the number of samples and m denoted the point number of each spectrum. The principal components matrix would be obtained from the matrix by means of PCA, and thereafter the hat matrix could be calculated as follows:

(2.2)

The diagonal element of the hat matrix was the leverage value of sample i, as shown in Formula 2.3, where i is any of one sample.

(2.3)

The leverage value of the sample calculated by this method demonstrated the effects of the sample on the model. Therefore, if the leverage value of sample i was too high, two or three times higher than the average value of the sample set, it indicated that the regression model relied on sample i too much, which was bad for the stability of the model. In other words, sample i might be an abnormal sample. The student residual of a sample demonstrated the predicting ability of the model to this sample. When the number of the principal components was f, the student residual of the sample i was calculated as follows:

(2.4)

where is the prepared concentration, is the predictive concentration, is the leverage value.

RBF

RBF can be treated as a question of the surface fitting (approximation) in a high-dimensional space, it can look for a surface best matching the training data in a multidimensional space by learning, the new data inputted into the model are analyzed by the surface. The rationale of RBF is to carry on the back propagation learning by using the recursive technology called stochastic approximation in statistics. RBF neural network is used to build model using the command “newrbe” in this study as its good robustness and self-study ability, it is composed of input layer, output layer, and hidden layer, where the Radbas and a pure linear function (Purelin) were set as the transfer function of hidden layer and output layer, respectively, and the speed was set as 1000.

Model evaluation

R, RMSEP, and ratio of performance to standard deviate (RPD) were taken as indexes to evaluate the model and calculated as follows:

(2.5)

(2.6)

(2.7)

where is the value measured by the standard method, is the value predicted by model, is the average value of preparing values, m is the sample number of the prediction set, SD is the standard deviation of the calibration set. The bigger R was and the lower RMSEP was, the higher accuracy of the model was. RPD was used to validate the stability and predictive capacity of the model. RPD > 3 demonstrated that model had a great stability and prediction capacity.^[27]

Results and discussion

The Pre-Processing of Spectral Data

The spectral data were pre-processed with different combinations of the first derivative, second derivative, and the standard normal variate (SNV), respectively, to discard irrelevant spectral elements and reduce the amounts of noise and non-significant information. Table 2 showed the predictive results from RBF model with different pre-processing methods, it could be concluded that when the spectral data in calibration set were pre-processed, the RBF models could achieve good prediction precision and all better than the model having no spectral pre-processing. This was because different pre-treatments of spectra data were performed to evaluate their ability to remove or reduce disturbing effects not related to the chemical absorption of light.^[28] Investigated from Table 2, the highest R and RPD, and the smallest RMSEP occur in the quantitative prediction by means of RBF model set up by SNV spectra, which showed SNV spectra had the priority over conventional spectra, the first and second derivative spectra. Therefore, the RBF model had the optimal outcome when the spectra data were only pre-processed by SNV.

TABLE 2 The predictive abilities of RBF model with different pre-processing means

Pre-processing means	R	RMSEP	RPD
Non-pre-processing	0.9826	23.0029	5.4504
First derivative +SNV	0.9917	17.7588	7.6803
Second derivative +SNV	0.9849	21.1271	5.7185
First derivative	0.9905	21.3118	7.0685
Second derivative	0.9835	24.7798	5.4139
SNV	0.9937	15.5095	8.8216

The Division of Sample Sets

The sample set partitioning based on joint x-y distance (SPXY) algorithm was used to divide sample sets in this study, the 100 samples were divided into two parts in a ratio of 2:1 around 65 representative samples as calibration set and other samples as prediction set, with the former being used to design the calibration model and the latter being used to inspect and evaluate the performance of the model. Table 3 showed the mean values and the standard deviations of samples in different sets, the data of calibration set were close to that of all samples, which demonstrated that the calibration set had good representativeness, and the mean value and the standard deviation of the prediction set had a big deviation with that of all samples, which illustrated that the prediction set could validate the applicability of model. Figure 3 showed the sample distribution of the calibration set and prediction set. As a general observation of Fig. 3, it is noted that the samples of the calibration set distributed in all samples evenly. Above all, SPXY algorithm was a reasonable method to divide the samples into a calibration set and prediction set.

TABLE 3 The mean value and the standard deviation

The samples’ set	Range (mg·Kg^–1)	Mean value (mg·Kg^–1)	Standard deviation (mg·Kg^–1)
The all samples	10~505	257.5	145.06
The calibration set	10~505	250.45	153.38
The prediction set	50~500	271.18	128.44

FIGURE 3 The distribution of the calibration set and prediction set. (*) the prediction set, (∙) the calibration set.

The Elimination of Abnormal Samples

The leverages and student residuals of all samples were calculated, Fig. 4 was the distribution map of them with the abscissa being the leverage value and the ordinate being the student residual. Notice from Fig. 4 that the location of the sample with number S-001 was center-right compared with other samples, and its leverage value was relatively big; the location of the samples with number S-010, S-050, S-059, S-076, and S-081 deviated from the abscissa axis compared with other samples, and the student residuals were relatively big; the leverage value and student residual of the sample with number S-095 were both relatively big; thus those seven samples above could be categorized temporarily as abnormal samples.

FIGURE 4 The distribution of T detection of leverage and student residual of all the samples.

In order to determine the influence of abnormal samples in the model, the samples would be recycled one by one. Before recycling, it was noted that there was one tentative abnormal sample in prediction set and the other six in the calibration set, since it was uncertain whether they were abnormal or not, pick them out from sets and recycle them one by one to calibration set. First, the sample S-001 was recycled from the seven tentative abnormal samples into calibration set and the other six samples were eliminated, so far, there were 60 samples in the calibration set and 34 samples in the prediction set. Then, prediction models were designed by means of PLS regression, BP neural network, and RBF, respectively, to identify the effects of this sample on models. Finally, the other six samples would be processed likewise and the results were shown in Table 4.

TABLE 4 The models performance after recycling abnormal samples one by one

	PLS			BP			RBF
Abnormal sample elimination	R	RMSEP	RPD	R	RMSEP	RPD	R	RMSEP	RPD
Non-eliminating	0.9896	20.0547	6.4042	0.9458	48.4444	2.6512	0.9906	17.4120	7.3762
Eliminating seven samples	0.9670	40.7940	3.1484	0.9321	54.2962	2.3655	0.9913	16.7741	7.6568
Recycling the S-001	0.9915	18.9090	6.7963	0.9689	44.3447	2.8963	0.9913	16.7351	7.6746
Recycling the S-010	0.9670	40.4716	3.1735	0.9111	57.90732	2.2179	0.9911	16.9843	7.5620
Recycling the S-050	0.9671	40.9219	3.1385	0.9019	60.4141	2.1259	0.9914	16.6523	7.7128
Recycling the S-059	0.9673	40.4257	3.1771	0.9230	51.8434	2.4774	0.9911	16.9645	7.5708
Recycling the S-076	0.9662	41.3527	3.1058	0.9187	58.9225	2.1797	0.9911	16.9566	7.5743
Recycling the S-081	0.9663	41.3359	3.1071	0.9246	51.7246	2.4831	0.9912	16.9267	7.5877
Recycling the S-095	0.9891	19.8177	6.4808	0.9727	43.55162	2.9490	0.9915	16.5668	7.7526

Analyzing the data in Table 4, the recycle of samples S-010, S-059, S-076, and S-081 all decreased the performance of PLS model, BP model, and RBF model with R and RPD decreasing and RMSEP increasing, so those four samples were regarded as abnormal samples and should be eliminated. On the contrary, the recycle of sample S-001 and sample S-095 had improved the performance of the three models, thus they should be recycled. As for sample S-050, though it decreased the performance of PLS model and BP model, it had contributions to the RBF model, besides, RBF model was adopted in this study, so sample S-050 was recycled at last.

The Design of Model RBF and the Analysis of the Results

Before modelling, the residual 96 samples were divided into calibration set of 60 samples and prediction set of 36 samples, the two sets were pre-processed by SNV. The RBF model was designed with the spectral data of calibration set as input and the corresponding concentrations of BPO as output, and then it was tested by prediction set. Meanwhile, as a comparison, the models based on PLS and BP neural network were designed. The predictive results of those three models were shown in Table 5, as could be seen from the table that the RBF model was the best among those three models by comparing the values of index R, RPD, and RMSEP. The fitting curve between the predicted values of RBF model and the preparing values was shown in Fig. 5, where the diagonal (A = T) meant 100% correct prediction. The more consistent the fitting curve and the diagonal and the closer the data points near the diagonal, the better the prediction performance of the model. It was demonstrated in Fig. 5 that the predictive values of samples of prediction set evenly distributed on both sides of diagonal and the fitting curve matched well with the diagonal, which illustrated good predictive effects of the RBF model. Table 6 showed the preparing values, predictive values, and the absolute errors of 36 samples of prediction set, the data demonstrated that the absolute error of each sample was relatively low except for four samples with the concentrations 50, 100, 175, and 290 mg·Kg^–1, respectively, which could be accepted in actual testing, which further verified the better predictive effect of the model.

TABLE 5 The predictive abilities of the different models by the pre-processing SNV

Predictive model	R	RMSEP	RPD
PLS	0.9890	19.8895	6.4574
BP	0.9729	42.1205	3.0492
RBF	0.9937	15.5095	8.8216

TABLE 6 The predicted results

The number	Test value (mg·Kg^–1)	Predictive value (mg·Kg^–1)	Absolute error	The number	Test value (mg·Kg^–1)	Predictive value (mg·Kg^–1)	Absolute error
1	50	10.14917	39.85083	19	290	319.4436	–29.4436
2	70	75.03393	–5.03393	20	305	295.1428	9.857182
3	75	73.43774	1.56226	21	315	329.179	–14.179
4	80	74.20337	5.796635	22	340	357.3533	–17.3533
5	95	113.9548	–18.9548	23	350	362.7771	–12.7771
6	100	133.5085	–33.5085	24	360	379.406	–19.406
7	115	105.0466	9.953373	25	365	359.1838	5.816166
8	120	114.7278	5.272221	26	370	361.677	8.323002
9	130	117.7053	12.29468	27	380	375.1965	4.803471
10	135	114.0774	20.92261	28	395	395.6013	–0.60131
11	175	151.7439	23.25611	29	400	420.4514	–20.4514
12	200	183.4998	16.50025	30	405	409.5237	–4.52368
13	210	199.9495	10.05054	31	425	416.178	8.822025
14	245	258.0813	–13.0813	32	440	430.5398	9.460209
15	250	254.9324	–4.93237	33	465	470.5422	–5.54223
16	255	253.9182	1.081791	34	470	445.1765	24.82349
17	275	284.7102	–9.7102	35	475	482.7229	–7.7229
18	285	280.8308	4.169193	36	500	487.9783	12.02173

FIGURE 5 The comparison between the predicted value and the test value.

Conclusion

The NIRS technology was adopted to detect BPO in flour taking place of the conventional chemical measurements, which had greatly improved the analysis speed of the measurement and provided theoretical basis for the speedy detection of BPO. The prediction model was designed based on RBF, abnormal samples were detected using T detection of leverage and student residual, and the predictive accuracy was improved by means of SNV. If spectra data were not pre-processed and the abnormal samples were not eliminated, R, RMSEP, and RPD were 0.9826, 23.0029, and 5.4504, respectively; then, spectra data were pro-processed by SNV, R, RMSEP, and RPD were 0.9906, 17.412, and 7.3762, respectively; on this basis, the abnormal samples were eliminated by T detection of leverage and student residual, R, RMSEP, and RPD were 0.9937, 15.5095, and 8.8216, respectively. It was noted from the evaluation indexes that the precisions were higher after pre-processing spectra data and eliminating abnormal samples. The implementation of the proposed innovative optical approach could be utilized for fast quality control purpose and greatly ensure safety of flour that was not easy reach with conventional techniques.

AcknowledgmentS

The deepest gratitude goes foremost to the Agricultural Products Quality and Safety Research Institute of Heilongjiang Academy of Agricultural Sciences. Thanks for their support on the supply of the materials (flour and BPO) and the technical assistance (the measurement of BPO by physical and chemical methods). At last I also owe my sincere gratitude to my friends and my fellow classmates who helped me overcome the difficult course of the thesis writing.