Quantitative structure–affinity relationship study of azo dyes for cellulose fibers by multiple linear regression and artificial neural network

Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

Contents lists available at ScienceDirect

Chemometrics and Intelligent Laboratory Systems journal homepage: www.elsevier.com/locate/chemolab

Quantitative structure–affinity relationship study of azo dyes for cellulose fibers by multiple linear regression and artificial neural network Xiaojun Wang a, Yongyuan Sun a, Lei Wu a, Shaojin Gu a, Ruina Liu a, Li Liu a, Xin Liu a, Jie Xu a,b,⁎ a b

College of Materials Science & Engineering, Wuhan Textile University, 430200 Wuhan, China State Key Laboratory of New Textile Materials & Advanced Processing Technology, Wuhan Textile University, 430200 Wuhan, China

a r t i c l e

i n f o

Article history: Received 12 November 2012 Received in revised form 23 January 2014 Accepted 2 March 2014 Available online 11 March 2014 Keywords: QSPR Affinity Azo dyes Cellulose fiber Multiple linear regression Artificial neural network

a b s t r a c t A quantitative structure–property relationship study was performed to correlate descriptors representing the molecular structures to fiber affinities for azo dyes. The complete set of 51 compounds was divided into a training set of 41 compounds and a test set of 10 compounds by using DUPLEX algorithm. Multiple linear regression analysis was used to select the best subset of descriptors and to build linear models; nonlinear models were developed by means of artificial neural network. The robustness of the obtained models was assessed by different approaches, including leave-many-out cross-validation, Y-randomization test, and external validation through test set. The obtained models with four descriptors show good predictive power: for the test set, a squared correlation coefficient (r2) of 0.916 was achieved by the linear model; while the nonlinear model with r2 of 0.935 for the test set performs better than the linear model. Furthermore, the applicability domain of the models was analyzed based on the Williams plot. The donor atom number for H-bonds, group polarizability and electronegativity of the dye molecules are found to play important roles for the dye-fiber affinity. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Dyes are widely used in textile dyeing, paper printing, color photography, pharmaceutical, food, cosmetics, and other industries [1–3]. The history of dyes can be traced back to 2600 BC when there was the earliest written record of the use of dyes in China [4]. Originally, the dyes were obtained naturally until William Henry Perkin prepared the first synthetic dye pigment “Mauve” from coal-tar chemicals in 1856 [5]. Since the first azo dye spliced onto the fabric by coupling in 1875, thousands of organic chromogens were produced by the dye industry [4]. For example, in 1995, 90% of the 3000 compounds registered in the Color Index were used at the level of at least 100 t per year [6]. Currently, the world market of dyes and pigments is estimated to be around 1.3 million tons, valued at US$16 billion. With respect to both number and production volume, azo dyes are the largest group of colorants constituting 70% of all organic dyes produced in the world [7]. Cotton/cellulose fiber is a tremendous natural resource that has broad application in various productions including the textile industry. Cellulose is an isotactic β-1,4-polyacetal of cellobiose, in which two glucose molecules are linked together via a glycosidic bond. Hydrogen bond network between the polymer components is formed, involving hydroxyl groups and ether oxygen (glycoside oxygen). Cellulose does not contain ionic sites, and as such is not able to fixate cationic or ⁎ Corresponding author at: State Key Laboratory of New Textile Materials & Advanced Processing Technology, Wuhan Textile University, 430200 Wuhan, China. E-mail address: [email protected] (J. Xu).

http://dx.doi.org/10.1016/j.chemolab.2014.03.001 0169-7439/© 2014 Elsevier B.V. All rights reserved.

anionic dyes through Coulomb forces [8]. Thus, direct dyeing of cotton requires specific features of the dye molecules to result in an overall affinity for the cellulose fiber. For a successful uptake from aqueous solution and stable fixation on the textile fiber, the dye molecules with hydrogen bond donor or acceptor sites are adsorbed in intermicellar cavities, where they are thought to form large aggregates and suitable dye–fiber interaction. The interaction between the dye molecule and the cellulose fiber is affected by various factors (electrostatic fields, Van der Waals forces, formation of hydrogen bonds, hydrophobicity, etc.), and the available experimental data are rather diverse. That is why theoretical methods for calculation of quantitative structure–property/activity relationships (QSPR/QSAR) are widely used to study these interactions [8–14]. The QSPR/QSAR approach is based on the assumption that the variation of the behavior of the compounds, as expressed by any measured physicochemical properties, can be correlated with numerical changes in structural features of all compounds [15–19]. The advantage of this approach lies in the fact that it requires only the knowledge of the chemical structure and is not dependent on any experimental properties. Once a correlation is established and validated, it can help to synthesize new compounds with high affinity. Classical quantitative structure–affinity relationships have been reported in the literature for the dye affinities on the cellulose fiber [9,10]. Zhokhova et al. [10] developed a multiple linear regression (MLR) model with the squared correlation coefficient R2 of 0.971 and the cross-validated R2 (Q2) of 0.949 for the affinities of azo dyes based fragment descriptors. However, this model has not been evaluated with the external test set. In fact, validation is a crucial aspect

2

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

of any QSPR/QSAR [20]. Based on the hypothesis of specific dye–fiber interactions and as an alternative to classical QSPR studies, Comparative Molecular Field Analysis (CoMFA) [8,11] and Comparative Molecular

Surface Analysis (CoMSA) [12,13] were also used to predict dye adsorption properties. Schüürmann et al. [8] established models for the azo dye-fiber affinities using CoMFA without external validation.

Table 1 Experimental and calculated dye affinity (A, kJ mol−1) of azo dyes. No.

X

Ra

Expt. A

Cacl. A

No.

MLR

ANN

X

Ra

Expt. A

Cacl. A MLR

ANN

1

γb

22.26

21.05

22.03

27

C

9.45

9.20

9.61

2

γb

15.69

14.72

15.04

28b

B

9.2

7.92

7.74

3

γa

14.35

15.24

15.39

29

C

9.03

8.36

8.61

4

E

9.62

12.09

10.28

30

A

8.78

8.41

8.32

5

E

8.79

7.57

7.37

31

C

8.4

10.05

10.48

6

γb

13.18

14.13

12.55

32

C

8.28

8.41

8.76

7b

γb

10.92

12.33

11.30

33

B

7.15

8.15

8.32

8b

C

14.48

15.31

16.52

34

D

7.06

6.02

5.88

9

D

10.5

11.49

9.79

35

A

7.02

7.82

7.45

10b

D

7.7

7.60

7.58

36

B

6.52

7.05

6.93

11

E

5.23

6.58

6.32

37

E

6.27

5.99

5.69

12

γb

8.58

10.14

9.53

38

D

6.23

5.18

5.16

13

C

13.56

10.91

11.90

39b

D

6.02

6.87

6.72

14

E

4.48

4.39

3.87

40

E

5.81

6.19

5.81

15

E

4.6

4.39

3.87

41

D

5.18

5.23

5.22

16b

D

3.59

6.10

5.95

42

E

5.1

5.27

5.12

17

D

1.92

2.36

2.16

43

C

4.64

6.07

6.00

18

D

2.97

3.91

3.56

44b

E

4.26

5.19

5.04

19

C

9.49

9.28

9.28

45

B

4.22

4.20

4.18

20

C

7.24

7.60

7.29

46

C

4.1

6.00

5.82

21

C

6.61

6.05

5.42

47

B

4.05

4.13

4.06

22b

A

15.8

13.73

14.72

48

D

3.85

2.26

3.30

23

A

14.25

13.17

13.73

49

D

3.43

2.83

3.54

24b

A

13.08

12.36

13.21

50

E

3.22

2.34

3.28

25

A

12

10.10

10.73

51b

E

2.84

2.27

3.26

26

B

8.48

8.40

9.66

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

Their models have very high prediction accuracy (standard error of estimation s = 0.68), but the reliability is suspectable due to the poor Q2 (0.75). Polanski et al. [12] obtained very predictive models (mean relative error (MRE) = 11.93%) using the CoMSA method. Recently, FunarTimofei et al. [14] established QSPR models for the affinities of 21 azo dyes by MLR and CoMFA, but the MRE values are quite large (31.16% and 25.78% for the training and test sets, respectively). The goal of this work is to present a complete quantitative structure– affinity relationship study for a series of azo dyes by MLR and artificial neural network (ANN) methods according to the OECD Principles [21]. ANN-based modeling methods are expected to produce more accurate QSPR models compared to linear regression methods, since they have the ability to handle the possible nonlinear relationships during the training process [22]. Dye structures were modeled by molecular mechanics and molecular descriptors were derived from the optimized structures to correlate with the dye–fiber affinities. Linear and nonlinear models were developed, which could predict the affinities and help to understand the dye adsorption mechanism on fiber. The predictive ability and robustness of the obtained models were fully assessed by different approaches, including leave-many-out cross-validation, Y-randomization test, and external validation through test set. Additionally, the applicability domain (AD) of the models was analyzed based on the Williams plot. 2. Materials and method 2.1. Data set

2.2. Descriptor generation The chemical structure of each compound was sketched using the HyperChem program [23]. The conformational analysis was performed according to the reported procedure [24], using the MM + molecular mechanics method (Polak–Ribiere algorithm) and the semi-empirical AM1 method at a restricted Hartree–Fock level with no configuration interaction. A gradient norm limit of 0.01 kcal · Å−1 · mol−1 was chosen as the stopping criterion. The resulting lowest energy conformer for each compound was then submitted to the Dragon software [25] to calculate totally 1664 molecular descriptors. These descriptors include (1) 0D-constitutional (atom and group counts); (2) 1D-functional groups and atom centered fragments; (3) 2D-topological, BCUTs, walk and path counts, autocorrelations, connectivity indices, information indices, topological charge indices, and eigenvalue-based indices; and (4) 3D-Randic molecular profiles from the geometry matrix, geometrical, WHIM, and GETAWAY descriptors. To reduce redundant and useless information, constant or near constant values and descriptors found to be highly correlated pairwise (one of any two descriptors with a correlation greater than 0.99 [26]) were removed. Furthermore, the rug-like descriptors (including Notes to Table 1: a R:

.

b

Compounds in the test set.

GVWAI-80, Inflammat-80, Psychotic-80, Hypertens-80, Neoplastic-80, and Infective-80) were excluded. Finally, 797 descriptors were remained to undergo subsequent descriptor selection. 2.3. Data set splitting The data set was divided into training and test sets by the DUPLEX algorithm [27] combined with principal component analysis (PCA) due to its easy implementation and its good performance in the selection of representative training set samples. The PCA was first performed based on the 797 descriptors of the whole dataset, and the obtained principal components were put into the DUPLEX algorithm to select a training set of 41 compounds and a test set of 10 compounds. 2.4. Model development Stepwise MLR analysis combined with leave-one-out (LOO) crossvalidation was used to select the best subset of descriptors and to develop linear QSPR models based on the training set. Step-by-step descriptors were added to the equation, and a new regression was performed. If the new descriptor contributed significantly to the regression equation, it was retained; otherwise, it was excluded, hence preventing over-fitting. F-to-enter and F-to-remove were 4 and 3, respectively. In addition, a variance inflation factor (VIF) was calculated to test if multicollinearities existed among the descriptors, which is defined as VIF ¼

The experimental dye affinity (A) values of the 51 azo dyes were collected from the literature [11,14]. The compounds and their corresponding affinities are listed in Table 1.

3

1 1−R2j

ð1Þ

where R2j is the squared correlation coefficient between the jth coefficient regressed against all the other descriptors in the model. If VIF equals to one, no intercorrelation exists for each descriptor; if VIF maintains within the range 1.0–5.0, the corresponding model is acceptable; if VIF is larger than 10.0, the corresponding model is unstable and rejected [28]. Nonlinear models were then developed by submitting the selected descriptors from MLR to a three-layer, fully connected, feed-forward ANN. The number of input neurons was equal to that of the descriptors in the linear model. The number of hidden neurons was optimized by trial and error procedure on the training process. One output neuron was used to represent the experimental A. The network was trained using the quasi-Newton BFGS (Broyden–Fletcher–Goldfarb–Shanno) algorithm [29]. To avoid overtraining, one tenth data from the training set were randomly selected as separate validation set to monitor the training process; that is, during the training of the network the performance was monitored by predicting the values for the systems in the validation set. When the results for the validation set ceased to improve, the training was stopped. 2.5. Model evaluation and validation Model performance was evaluated by the following statistical parameters: the R2, the Q2, the adjusted R2, the s, the F ratio values,

4

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

and the significance level value p. The adjusted R2 is calculated using the following formula:

Selected number

0.90

ð2Þ

0.85 0.80 0.75

2

where N is the number of members of the training set and M is the number of descriptors involved in the model. The adjusted R2 is a better measure of the proportion of variance in the data explained by the correlation than R 2 (especially for correlations developed using small data sets) because R2 is somewhat sensitive to changes in N and M. The adjusted R2 corrects for the artificiality introduced when M approaches N through the use of a penalty function which scales the result. The predictive performance of the models was measured by mean absolute error (MAE) and mean relative error (MRE), defined as: 1X MAE ¼ jAð expt:Þ−Aðcalc:Þj n

0.95

Q

2

Radj

 
 N−1 2 1−R ¼ 1− N−M−1

1.00

0.70 0.65 0.60 0.55 0.50 1

2

3

4

5

6

7

Number of descriptors

ð3aÞ

Fig. 1. Q2LOO vs. the number of descriptors in the MLR equation for the training set.

MRE ¼

  1 XAð expt:Þ−Aðcalc:Þ  : n Að expt:Þ

ð3bÞ 2

The reliability and robustness of the models were assessed by the following techniques: Y-randomization test, leave-many-out (LMO) cross-validation, and external validation through test set. 2.5.1. Y-randomization test Y-randomization tests were carried out to prove the possible existence of chance correlation during the MLR model development. In this procedure, the dependent variable (Y) was randomly scrambled and used in the regression. Models were then investigated again with all members in the descriptor pool to select the best subset of descriptors. Here, fifty randomization runs were done. The resulting models obtained from the training set with the randomized A values should have significantly lower R2 and Q2 values than the proposed one because the relationship between the structure and response is broken. This is a proof of the proposed model's validity as it can be reasonably excluded that the originally proposed model was obtained by chance correlation. 2.5.2. LMO cross-validation The internal predictive ability and robustness of the developed models were evaluated by LMO cross-validation. The training set of N samples is divided into consecutive blocks of M samples, where the first M defines the first block, the following M samples are the second block, and so on [30,31]. LMO was performed for M = 2, 3, etc. This process was repeated twenty times by scrambling the samples of descriptors and dependent variable (X and Y) simultaneously for each value of M, and the average value of Q2LMO was obtained. 2.5.3. External validation Validation of the models was further performed by using the external test set composed of data not used to develop the prediction models. The external R2CV,ext for the test set is determined with the following equation: X 2

RCV;ext ¼ 1−X

2 ðyi −e yi Þ 2

ðyi −ytra Þ

2

ð5bÞ

    2 2 2 2 02 2 r −r 0 =r b0:1 or r −r 0 =r b0:1

ð5cÞ

0

0:85≤k ≤1:15 or 0:85≤k ≤1:15   X ðyi −yÞ e yi −e y r ¼ rffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi 2 X X 2 e ðyi −yÞ yi −e y

here :

2 r0

X

r 2 e yi 0 yi −e ¼ 1− X  2 e yi −e y X

02 r0

¼ 1− X

r

yi −yi 0

ð5dÞ

ð6aÞ

ð6bÞ

2 2

ðyi −yÞ

ð6cÞ

X

yi e yi k¼ X 2 yi

ð6dÞ

X

yi e yi k ¼ X 2 e yi 0

ð6eÞ

yr0 ¼ ky and yr0 ¼ k0 e y, respectively. where e yr0 and yr0 are defined as e ð4Þ

yi are the observed and the calculated response values, where yi and e respectively; and ytra is the averaged value for the response variable of the training set. According to Golbraikh and Tropsha [20], a QSPR model is successful if it satisfies the following criteria: RCV;ext N0:5

r N0:6

ð5aÞ

2.6. Applicability domain analysis The applicability domain (AD) of a QSPR model [32,33] must be defined if the model is to be used for screening new compounds. In this work, the structural AD was verified by the leverage approach. The leverage hi [33] is defined as follows:  −1 T T xi hi ¼ xi X X

ð7Þ

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

1.0

a

0.8

R

2

0.6

0.4

0.2

0.0 0.0

0.2

0.4

0.6

0.8

1.0

Correlation

b

1.0

3.1. MLR model Stepwise MLR analysis was used to select the best subset of descriptors and to develop linear models on the training set. The plot of Q2LOO (coefficient correlation of LOO cross-validation) versus the number of selected descriptors is shown in Fig. 1. The optimum subset size was reached when adding another descriptor did not improve the performance of the model significantly. Through this procedure, the four-parameter model was selected as the best model:

2

0.6

Q

The Williams plot, the plot of leverage values versus standardized residuals, was used to give a graphical detection of both the response outliers (Y outliers) and the structurally influential compounds (X outliers). In this plot, the two horizontal lines indicate the limit of normal values for Y outliers (i.e. samples with standardized residuals greater than 3.0 standard deviation units, ±3.0 s); the vertical straight lines indicate the limits of normal values for X outliers (i.e. samples with leverage values greater than the threshold value, h N h*). For a sample in the external test set whose leverage value is greater than h*, its prediction is considered unreliable, because the prediction is the result of a substantial extrapolation of the model. Conversely, when the leverage value of a compound is lower than the critical value, the probability of accordance between predicted and experimental values is as high as that for the compounds in the training set. It is noteworthy that the response outliers can be highlighted only for compounds with known responses and the possibility of a compound to be out of the structural AD of a model can be verified for every new compound, the only knowledge needed being the molecular structure information represented by the molecular descriptors selected in the model. 3. Results and discussion

0.8

0.4

0.2

A ¼ 19:908−51:675GATS3e þ 64:114GATS2p−2:115GGI3 þ 2:269nHDon

0.0 0.0

0.2

0.4

0.6

0.8

2

2

2

Fig. 2. R and Q vs. correlation coefficient between the original and permuted response data.

where xi is the descriptor row-vector the i-th compound, xTi is the transpose of xi, X is the descriptor matrix, XT is the transpose of X. The warning leverage h* is generally fixed at 3(M + 1)/N, where N is the total number of samples in the training set and M is the number of descriptors involved in the correlation. 1.0

0.8

2

0.6

0.4

0.2

0.0 0

4

8

12

16

M Fig. 3. Plot of Q2LMO vs. the M number for the LMO cross-validation.

ð8Þ

1.0

Correlation

QLMO

5

20

2

2

N ¼ 41; R ¼ 0:928; Q LOO ¼ 0:900; Radj ¼ 0:920; s ¼ 1:165; F ¼ 116:2; pb0:00001 where GATS3e is the Geary autocorrelation – lag 3/weighted by atomic Sanderson electronegativities; GATS2p is the Geary autocorrelation – lag 2/weighted by atomic polarizabilites; GGI3 is the topological charge index of order 3; and nHDon is the number of donor atoms for H-bonds (N and O), respectively. More detailed information about these descriptors can be found in Dragon software user's guide [25] and the references therein. In general, the larger the magnitude of the F ratio, the better the model predicts the property values in the training set. The large F ratio of 116.2 indicates that Eq. (8) does an excellent job in predicting the A values. Eq. (8) has an adjusted R2 value of 0.920, which indicates good agreement between the correlation and the variation in the data. The model was further validated by randomization tests. The obtained R2yrand and Q2yrand vs. the correlation coefficient between the original and permuted response data are plotted in Fig. 2. The results for all randomized models have bad quality when compared to the original model. The intercepts for R2yrand and Q2yrand were regressed to be 0.196 and − 0.043, respectively, which are within the limit values recommended in the literature, i.e., below 0.3 for R2yrand and 0.05 for Q2yrand [30]. Thus, the original model is considered being free of chance correlation. The robustness of the model was checked using LMO cross-validation. The average Q2LMO values vs. the M number are plotted in Fig. 3. It can be found that the Q 2LMO is still stable and very close to the respective Q2LOO even at M = 14, confirming the robustness of the model [30]. Some important statistical parameters (given in Table 2) were used to evaluate the involved descriptors. The t-value of a descriptor measures the statistical significance of

6

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

Table 2 Characteristics of the descriptors selected in the optimal MLR model. Descriptor Constant GATS3e GATS2p GGI3 nHDon

Descriptor type

Coefficient

Error

t-Value

t-Probability

VIF

2D autocorrelations 2D autocorrelations Topological charge indices Functional group counts

19.908 −51.675 64.114 −2.115 2.269

4.952 4.023 4.387 0.436 0.186

4.020 −12.845 14.615 −4.848 12.178

0.0000 0.0000 0.0000 0.0000 0.0000

1.126 1.117 1.314 1.407

the regression coefficient. The high absolute t-values shown in Table 2 express that the regression coefficients of the descriptors involved in the MLR model are significantly larger than the standard deviation. The t-probability of a descriptor can describe the statistical significance when combined together within an overall collective QSPR model (i.e., descriptors' interactions). Descriptors with t-probability values below 0.05 (95% confidence) are usually considered being statistically significant in a particular model, which means that their influence on the response variable is not merely by chance [34]. The smaller t-probability suggests the more significant descriptor. The t-probability values of the three descriptors are very small, indicating that all of them are highly significant descriptors. The VIF values (less than five) suggest that these descriptors are weakly correlated with each other. The sign change analysis was also performed for these descriptors and the results are listed in Table 3. According to the published criterion I [31], these descriptors are considered being free of the sign change problem. Thus, the model can be regarded as an optimal regression equation. The calculated affinity values from the MLR model for the training and test sets are given in Table 1 and Fig. 4. The errors are distributed on both sides of the zero point, thus one may conclude that there is no systematic error in the model development. The MAE for the training and test sets is 0.881 and 1.126, respectively, resulting in a MAE of 0.929 for the complete data set. The following statistical parameters were obtained for the test set, which obviously satisfy the generally accepted condition and thus demonstrate the predictive power of the present model: 2

RCV;ext ¼ 0:915N0:5

of ANN commends it as a powerful tool for pattern classification and building predictive models. A particular advantage of ANN is its inherent ability to incorporate nonlinear dependencies between the dependent and independent variables without using an explicit mathematical function. Among the neural network learning algorithms, the backpropagation (BP) method [35] is one of the most commonly used methods. The drawback of BP is that the training processes slowly, because the gradient-descent algorithm is usually used for minimizing the sum-of-squares error. In this study, the quasi-Newton BFGS algorithm [29] was used. The advantages of the BFGS algorithm are that specifying rate or momentum is not necessary and the training is much more rapid [36]. The four descriptors from the MLR model were used as inputs to the network. The number of hidden neurons is an important parameter influencing the performances of the ANN. The usual rule of thumb is that the weights and biases should be less than the samples so that the model achieved by the network is stationary [37]. Thus, a 4-2-1 network architecture was obtained after trial and error procedure, with R2 of 0.955, Q2LOO of 0.950 and s of 0.885 for the training set, respectively. The predictive results from the ANN model for the complete data set are given in Table 1 and Fig. 5. The proposed ANN model is predictive as it satisfies the conditions for the test set: 2

RCV;ext ¼ 0:929 N0:5 2

r ¼ 0:935N0:6   2 2 2 r −r 0 =r ¼ ð0:935−0:997Þ=0:935 ¼ −0:066b0:1   2 02 2 or r −r 0 =r ¼ ð0:935−0:994Þ=0:935 ¼ −0:063b0:1 0

2

r ¼ 0:916 N0:6   2 2 2 r −r 0 =r ¼ ð0:916−0:999Þ=0:916 ¼ −0:091b0:1   2 02 2 or r −r 0 =r ¼ ð0:916−1:000Þ=0:916 ¼ −0:092b0:1

0:85 ≤k ¼ 0:965≤1:15 or 0:85≤k ¼ 1:022≤1:15: The MAE for the training and test sets is 0.685 and 0.947, respectively, giving rise to a MAE of 0.736 for the complete data set.

0

0:85 ≤k ¼ 0:993≤1:15 or 0:85≤k ¼ 0:989≤1:15: 25

3.2. ANN model 20

Table 3 Statistics for the sign change problem in correlation and regression vectors. Descriptor

rca

rtb

rec

βcd

βte

GATS3e GATS2p GGI3 nHDon

−0.539 0.696 −0.297 0.107

−0.522 0.656 −0.262 0.144

−0.606 0.923 −0.382 0.008

−0.613 0.648 −0.277 0.593

−0.609 0.690 −0.248 0.646

Acalc.(kJ/mol)

The ANN has become an important and widely used nonlinear modeling technique for QSPR studies. The mathematical adaptability 15

10

5

Training Test

a

Pearson correlation coefficient between a descriptor and the dependent variable Y for the complete dataset. b Pearson correlation coefficient between a descriptor and the dependent variable Y for the training set after data split. c Pearson correlation coefficient between a descriptor and the dependent variable Y for the external validation set after data split. d Normalized regression vector for the complete dataset. e Normalized regression vector for the training set after data split.

0 0

5

10

15

20

25

Aexpt. (kJ/mol) Fig. 4. Calculated vs. experimental affinity values by the MLR model for the complete dataset.

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

25

Acalc.(kJ/mol)

20

15

10

5

Training Test 0 0

5

10

15

20

25

Aexpt. (kJ/mol) Fig. 5. Calculated vs. experimental affinity values by the ANN model for the complete dataset.

3.3. Applicability domain of the developed models No matter how robust, significant and validated a QSPR model may be, it cannot be expected to reliably predict the modeled property for

7

the complete universe of compounds. QSPR models must always be verified for their applicability with regard to chemical domain, in order to produce predicted data that can be considered reliable only for chemicals not too structurally dissimilar with the training samples. In fact, in the case of structurally dissimilar molecules, the data predicted by the model must be judged as extrapolations. The chemical AD of the model and the reliability of the predictions are verified by the leverage approach. In Williams plot, chemicals influential on the structural domain of the model, characterized by a hat threshold (h* = 0.366, vertical dotted line in Fig. 6), can be explained as compounds with peculiar features poorly represented in the training set, which could affect the variable selection for a better modeling of those chemicals. Outliers, of which the standardized residual values exceed the cut-off value (here ± 3 s, horizontal dotted line in Fig. 6), could be associated with errors in the experimental values. On analyzing the AD of the models from the Williams plot (Fig. 6), there is no outlier compound both for the training and test sets and all the compounds are located within the AD of the QSPR models, which further indicates the reliability of the model. Although Sample 1 in the training set has a relatively high leverage value, it was well fitted. Thus, the sample that has high leverage value and is fitted well by the model could reinforce the model [38]. Due to their high predictive ability, the proposed models could be used to screen existing databases. In this case, the AD will serve as a valuable tool to filter out dissimilar chemical structures. 3.4. Comparison of the results obtained by different approaches

6

a

5 4

Standardized residuals

3 2 1 0 -1 -2 -3 -4 -5 -6 0.00

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

0.45

h 6

b

5

3.5. Descriptor analysis and interpretation

4 3

Standardized residuals

The results of the present study were compared with those of the previous models and are listed in Table 4. Zhokhova et al. [10] developed the MLR model with high R2 of 0.971 and Q2 of 0.949 (Model I in Table 4) for the affinities of azo dyes based fragment descriptors. However, this model has not been evaluated with the external test set. Schüürmann et al. [8] established QSPR models for the azo dyefiber affinities using CoMFA without external validation. Their model (Model III in Table 4) is superior to others in the prediction accuracy (s = 0.68), but its reliability is suspectable due to the poor Q2 (0.75). Polanski et al. [12] obtained very predictive models (MRE = 11.93%) using the CoMSA method. Recently, Funar-Timofei et al. [14] established QSPR models for the affinities of 21 azo dyes by MLR and CoMFA, but the MRE values are quite large (31.16% and 25.78%, respectively). The number of the observed dyes in the current work is more than for the previous models, whereas the predictive ability is also very accurate. In addition, we can see that the ANN model performs better than the MLR model, which further confirms the nonlinear relationship between the structural information and the dye–fiber affinity (in agreement with that demonstrated by the poor rc values in Table 3).

Based on a previously described procedure [39], the relative contributions (Ci) of the four descriptors to the MLR and ANN models were determined and are listed in Table 5. The significance of GATS3e, GATS2p and nHDon in the proposed models is quite close. The Geary autocorrelation GATSkw [40] is calculated by applying the Geary coefficient to the molecular graph:

2 1 0 -1 -2 -3

GATSkw ¼

-4 -5 -6 0.00

0.05

0.10

0.15

0.20

0.25

0.30

0.35

0.40

h Fig. 6. Williams plots of the developed model: (a) MLR and (b) ANN.

nSK X nSK  2 1 X  δ  wi −w j 2Δ i¼1 j¼1 ij nSK X 1 2  ðw −wÞ ðnSK−1Þ i¼1 i

0.45

where k is the lag, w is any atomic property (including atomic number, masses, van der Waals volumes, Sanderson electronegativities, and polarizabilities), w is its average value on the molecule, nSK is the

8

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9

Table 4 Comparison between the different QSPR models for the dye–fiber affinity.a Model

Algorithm

N

R2

Q2

s

r2

MRE

AD

I: Zhokhova et al. [10] II: Schüürmann et al. [8] III: Schüürmann et al. [8] IV: Polanski et al. [12] V: Funar-Timofei et al. [14] VI: Funar-Timofei et al. [14] VII: The current model 1 VIII: The current model 2

MLR MLR CoMFA CoMSA MLR CoMFA MLR ANN

30 30 30 30 21 21 51 51

0.971 0.92 0.97 0.930 0.947 0.960 0.928 0.955

0.949 0.88 0.75 0.933 0.914 0.707 0.901 0.953

0.95 1.02 0.68 1.09 1.58 1.12 1.16 0.88

NR NR NR 0.923 0.959 0.979 0.915 0.929

NR 22.83% NR 11.93% 31.16% 25.78% 13.92% (14.47%/13.54%) b 11.01% (12.33%/10.09%) b

NR NR NR NR NR NR R R

a b

CoMFA, Comparative Molecular Field Analysis; CoMSA, Comparative Molecular Surface Analysis; R, reported; NR, not reported. Values in parenthesis for the same 21 samples and 30 samples as Models V and IV, respectively.

number of non-hydrogen atoms, δij is the Kronecker delta (δij = 1 if dij = k, zero otherwise, dij being the topological distance between two considered atoms). Δ is the sum of the Kronecker deltas, i.e. the number of atom pairs at distance equal to k. The computations of these descriptors involve the summation of different autocorrelation functions corresponding to the different fragment lengths and lead to different autocorrelation vectors corresponding to the lengths of the structural fragments. A weighting component in terms of physicochemical property has been embedded in this descriptor. As a result, these descriptors address the topology of the structure or parts thereof in association with a selected physicochemical property. In the descriptor's nomenclature, the character k indicates the number of consecutively connected edges considered in its computation and is called as the autocorrelation vector of lag k (corresponding to the number of edges in the unit fragment). The GATS2p descriptor characterizes the molecular polarizability and the interactions between each pair of atom in the molecules. This descriptor can be related to how the polarizability is distributed along the topological structure of the molecule. The positive sign of this descriptor points out the favorable influence of molecular polarizability on the dye–fiber affinity. This suggestion corresponds with the results reported by Funar-Timofei et al. [14] that the dye polarity is favorable for binding to cellulose. The group electronegativity of molecular substituents also has a great effect on the interaction between dye and cellulose fiber, as demonstrated by the presence of GATS3e. In addition, the participation of GATS2p and GATS3e may be viewed in terms of association of affinity information content with the two and three centered structural fragments. However, further deciphering of the information content of these descriptors is very complex as their computations involve integration of the structural fragments and due to this it is not possible to traverse backward from a higher state to a lower one [41]. The number of donor atoms for H-bonds (nHDon) is a measure of the hydrogen-bonding ability of a molecule expressed in terms of number of possible hydrogen-bond donors. Specifically, it is calculated by adding up the hydrogens bonded to any nitrogen and oxygen in the molecule. The sign of this descriptor in the MLR model is positive, indicating that the azo dye with more donor atoms for H-bonds would have larger affinity values. It is well-known that there are many hydroxyl groups on the surface of the cellulose fiber. The contribution of this descriptor to the dye–fiber affinity agrees with the contribution that one could expect for the H-bonding interactions between the dye and the cellulose fiber [14]. The other participated descriptor, GGI3, belongs to Galvez class descriptors [42]. Galvez descriptors are the Galvez topological charge indices and have their origin in the first ten eigenvalues of the Table 5 Relative contributions of the descriptors to the proposed models. Descriptor

GATS3e

GATS2p

GGI3

nHDon

MLR Ci (%) ANN Ci (%)

27.93 28.15

28.52 28.76

14.32 18.13

29.24 24.96

polynomial of corrected adjacency matrix of the compounds. All the Galvez class descriptors belong to two categories. Of this one category corresponds to the topological charge index of order n (GGIn) and the other to the mean topological charge index of order n (JGIn), where ‘n’ represents the order of eigenvalue. The topological charge index of order 3 (GGI3) contributes negatively to the affinity and suggests that a lower value of it would increase the dye–fiber affinity of azo compounds. The GGI3 descriptor evaluates charge transfer in molecule carried out from the adjacency topological matrix. Extended charge distribution leading lower value of GGI3 can improve the affinity. 4. Conclusion In this paper, linear and nonlinear QSPR models for the fiber affinity of azo dyes were successfully developed by using MLR and ANN methods. The predictive ability and robustness of the obtained models were assessed by different approaches, including leave-many-out cross-validation, Y-randomization test, and external validation through test set. The applicability domain of the models was analyzed based on the Williams plot. The MLR model provides a transparent output, which indicates that the dye affinity for cellulose increases with the increased number of donor atoms for H-bonds in the dye molecule. Molecular polarizability and electronegativity also affect the dye binding to cellulose. By comparing the performance of the MLR and ANN models, it is proved that the ANN model predicts the property more accurately than the MLR model. Thus, it can be concluded that this investigation extends the research method to predict the azo dye-fiber affinity. Conflict of interest The authors declare no conflict of interest. Acknowledgments This work was supported by the National Natural Science Foundation of China (No. 51003082 and No. 61108033) and the Young People Project of Wuhan Science and Technology Bureau (No. 201271031381). References [1] S. Venkata Mohan, N. Chandrasekhar Rao, K. Krishna Prasad, J. Karthikeyan, Waste Manag. 22 (2002) 575–582. [2] E. Forgacs, T. Cserhati, G. Oros, Environ. Int. 30 (2004) 953–971. [3] A.B. dos Santos, F.J. Cervantes, J.B. van Lier, Bioresour. Technol. 98 (2007) 2369–2385. [4] A. Bafana, S.S. Devi, T. Chakrabarti, Environ. Rev. 19 (2011) 350–371. [5] D. Wesenberg, I. Kyriakides, S.N. Agathos, Biotechnol. Adv. 22 (2003) 161–187. [6] J.R. Easton, Colour in Dye House Effluent, Society of dyers and coloursit, Oxford, 1995. [7] C.M. Carliell, S.J. Barclay, C. Shaw, A.D. Wheatley, C.A. Buckley, Environ. Technol. 19 (1998) 1133–1137. [8] G. Schüürmann, S. Funar-Timofei, J. Chem. Inf. Comput. Sci. 43 (2003) 1502–1512. [9] S. Timofei, W. Schmidt, L. Kurunczi, Z. Simon, Dye. Pigment. 47 (2000) 5–16. [10] N.I. Zhokhova, I.I. Baskin, V.A. Palyulin, A.N. Zefirov, N.S. Zefirov, Russ. J. Appl. Chem. 78 (2005) 1013–1017.

X. Wang et al. / Chemometrics and Intelligent Laboratory Systems 134 (2014) 1–9 [11] S. Funar-Timofei, G. Schüürmann, J. Chem. Inf. Comput. Sci. 42 (2002) 788–795. [12] J. Polanski, R. Gieleciak, M. Wyszomirski, J. Chem. Inf. Comput. Sci. 43 (2003) 1754–1762. [13] J. Polanski, R. Gieleciak, M. Wyszomirski, Dye. Pigment. 62 (2004) 61–76. [14] S. Funar-Timofei, W.M.F. Fabian, L. Kurunczi, M. Goodarzi, S.T. Ali, Y.V. Heyden, Dye. Pigment. 94 (2012) 278–289. [15] J. Devillers, A.T. Balaban (Eds.), Topological Indices and Related Descriptors in QSAR and QSPR, Gordon and Breach, The Netherlands, 1999. [16] M. Karelson, Molecular Descriptors in QSAR/QSPR, Wiley-Interscience, New York, 2000. [17] X.J. Yao, Y.W. Wang, X.Y. Zhang, R.S. Zhang, M.C. Liu, Z.D. Hu, B.T. Fan, Chemom. Intell. Lab. Syst. 62 (2002) 217–225. [18] J. Xu, B. Guo, B. Chen, Q. Zhang, J. Mol. Model. 12 (2005) 65–75. [19] J. Xu, L. Wang, L. Wang, X. Shen, W. Xu, J. Comput. Chem. 32 (2010) 3241–3252. [20] A. Golbraikh, A. Tropsha, J. Mol. Graph. Model. 20 (2002) 269–276. [21] The principles for establishing the status of development and validation of (quantitative) structureactivity relationships [(Q)SARs], OECD Document ENV/JM/ TG, http://www.oecd.org/chemicalsafety/risk-assessment/37849783.pdf 2004. [22] A. Afantitis, G. Melagraki, K. Makridima, A. Alexandridis, H. Sarimveis, O. Iglessi-Markopoulou, J. Mol. Struct. (Thoechem) 716 (2005) 193–198. [23] Hypercube, Inc., Gainesville, 2000. [24] L. Kurunczi, S. Funar-Timofei, A. Bora, E. Seclaman, Int. J. Quantum Chem. 107 (2007) 2057–2065.

[25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42]

9

R. Todeschini, V. Consonni, A. Mauri, M. Pavan, TALETE srl, Milan, 2006. H. Liu, P. Gramatica, Bioorg. Med. Chem. 15 (2007) 5251–5261. R.D. Snee, Technometrics 19 (1977) 415–428. G.R. Famini, C.A. Penski, L.Y. Wilson, J. Phys. Org. Chem. 5 (1992) 395–408. M.D. Wessel, P.C. Jurs, Anal. Chem. 66 (1994) 2480–2487. R. Kiralj, M.M.C. Ferreira, J. Braz. Chem. Soc. 20 (2009) 770–787. R. Kiralj, M.M.C. Ferreira, J. Chemom. 24 (2010) 681–693. A. Tropsha, P. Gramatica, V.K. Gombar, QSAR Comb. Sci. 22 (2003) 69–77. A. Atkinson, Plots, Transformations, and Regression, Clarendon Press, Oxford, UK, 1985. L.F. Ramsey, W.D. Schafer, The Statistical Sleuth, Wadsworth Publishing Company, USA, 1997. P.A. Jansson, Anal. Chem. 63 (1991) 357A–362A. L. Xu, J.W. Ball, S.L. Dixon, P.C. Jurs, Environ. Sci. Chem. 13 (1994) 851–941. Y.-H. Qi, Q.-Y. Zhang, L. Xu, J. Chem. Inf. Comput. Sci. 42 (2002) 1471–1475. P. Gramatica, QSAR Comb. Sci. 26 (2007) 694–701. F. Zheng, E. Bayram, S.P. Sumithran, J.T. Ayers, C.-G. Zhan, J.D. Schmitt, L.P. Dwoskin, P.A. Crooks, Bioorg. Med. Chem. 14 (2006) 3017–3037. R.C. Geary, Inc. Stat. 5 (1954) 115–145. P. Broto, G. Moreau, C. Vandycke, Eur. J. Med. Chem. 19 (1984) 71–78. J. Galvez, R. Garcia, M.T. Salabert, R. Soler, J. Chem. Inf. Comput. Sci. 34 (1994) 520–525.