Prediction of Depuration Rate Constants for Polychlorinated Biphenyl Congeners.

The integral equation formalism polarizable continuum model (IEF-PCM) for solvent effects with the default solvent (water) and solvent parameters, together with the density functional theory method at 6-31G(d) level, was used to optimize molecular structures for polychlorinated biphenyl (PCB) congeners. Four molecular descriptors were selected to develop quantitative structure-activity relationship (QSAR) models for the depuration rate constants (k d) of 63 PCB congeners in a juvenile rainbow trout (Oncorhynchus mykiss). The optimal multiple linear regression (MLR) model has the correlation coefficient R of 0.933 and the root mean square (rms) error of 0.0681 for the total set of 63 PCB congeners. The support vector regression model has R of 0.953 and rms error of 0.0576 for the total set. Both the MLR and SVM QSAR models in this paper were accurate and acceptable compared with other QSAR models for the depuration rate of PCB congeners reported in references. Thus, applying IEF-PCM and B3LYP/6-31G(d) calculations for molecular descriptor derivation of PCB congeners is successful.

Introduction

Development of predictive quantitative structure–activity relationship (QSAR) models has been the concern of chemists, pharmacists, and environmental scientists. QSARs play a significant role in evaluating the biological activity, toxicity, and physicochemical property of chemicals that may not be synthetized.[1,2] One of the main aims of QSAR studies is to develop a mathematical correlation between the activity of interest and molecular descriptors, which reflect the quantitative molecular structure information. A successful QSAR enables the exploration of the chemical reaction mechanism or biochemical process.[3,4]

Polychlorinated biphenyls (PCBs) are serious environmental pollutants, which have the characteristics of environmental persistence, bioaccumulation and long-distance global migration.[5] These characteristics are closely related to their physical and chemical properties. Animal experiments have shown that PCBs can induce lesions of skin, liver, gastrointestinal system, nervous system, reproductive system, and immune system and even lead to carcinomas. With the increasing attention to environmental and health issues, PCBs have become the research focus in the field of environmental science and toxicology. There are 209 PCB congeners that possess similar molecular structures. It is unrealistic to determine the physical and chemical properties or biological activities of each PCB because of the costs involved. The QSAR can help to predict the physical and chemical properties or activities of PCBs.

Some researchers have carried out the QSAR studies of bioconcentration factors,[6] octanol/water partition coefficients,[7] and health effects[8] for PCBs. Niu et al. developed a 5-descriptor QSAR model for the photodegradation half-life (t1/2) of 22 PCBs in the n-hexane solution under UV irradiation.[9] The ratio (4.4 = 22/5) of the samples (22) to descriptors (5) used is too small and lower than the minimum value of 5.

The research group of Liu developed 2 QSAR models for the depuration rates (the logarithm of the depuration rate constants, log kd) of 34 PCBs. One model was based on 5 descriptors and 26 PCBs in the training set.[10] The other was obtained from the training set of 24 PCBs with 8 descriptors.[11] Both the models were validated with the respective test sets. Although the sample set (34 PCBs) is relatively small compared with the descriptors used in the models, their correlation coefficients R (0.941 and 0.954, respectively) are high.[10,11]

Liu et al. further used 6 quantum chemical descriptors to build a QSAR model for the depuration rate constants (log kd) of 62 PCBs measured at 285 K in juvenile rainbow trout (Oncorhynchus mykiss). The correlation coefficient R and standard deviation are 0.853 and 0.05, respectively. The statistical results of the models are acceptable. In addition, the descriptors used have definite physical meanings.[12]

Fatemi and Chahi developed QSAR models for the half-lives (log t1/2) for 62 PCBs measured at 281 K in juvenile rainbow trout. The models were based on four descriptors and statistical methods including linear and nonlinear regression techniques. The coefficients of determination R2 from the training set including 43 PCBs are 0.855–0.921. These models are accurate and acceptable.[13]

Descriptor calculation is a key step toward developing QSAR models. All these QSAR models stated above are based on the descriptors obtained from a single molecule in vacuum and gas phase.[6−13] Due to the solvent effect, the geometrical structure, frontier orbital energy, and charge distribution of solute molecules in gas and solution may be different, which affect the dissolution and chemical reactions of solutions.[14] In this paper, the integral equation formalism polarizable continuum model (IEF-PCM) was adopted to model the solvent effects in molecular geometry optimization. After that, Dragon 6.0[15] was used to calculate the molecular descriptors. Multiple linear regression (MLR) and support vector machine (SVM) algorithm were used to develop the QSAR models of depuration rate constants for PCB congeners.[16]

Results and Discussion

Correlation analysis between the depuration rate constants (kd) of 42 PCB congeners in the training set (see Table ) and 284 descriptors calculated with Dragon 6.0[15] was carried out with the MLR analysis in IBM SPSS Statistical 19. An optimal descriptor subset was obtained, which included four molecular descriptors: LOC, SpMax3_Bh(m), SpMin8_Bh(p), and SpMin5_Bh(s). The MLR equation is as followswhere N is the number of PCB congeners, R is the correlation coefficient, rms is the root mean square error, and F is the Fischer ratio. Descriptor characteristics and definitions are shown in Tables and 3, respectively. The MLR model was used to predict the depuration rate constants kd of 21 PCB congeners in the test set. The correlation coefficient R and rms error are 0.922 and 0.0631, respectively. The total set of 63 PCB congeners has R of 0.933 and rms error of 0.0681.

Table 1

Experimental and Calculated Values of kd of 63 PCB Congeners

no.	PCB congener	kd(exp.)	kd(MLR)	kd(SVM)
Training Set
1	2,2′,3,3′,5-pentachlorobiphenyl	0.3	0.36	0.36
2	2,2′,4,4′,5-pentachlorobiphenyl	0.3	0.30	0.30
3	2,3,3′,4,4′-pentachlorobiphenyl	0.3	0.37	0.33
4	2,3′,4,4′,5-pentachlorobiphenyl	0.3	0.31	0.31
5	2,2′,3,4,4′,5′-hexachlorobiphenyl	0.3	0.31	0.31
6	2,2′,3,4′,5,5′-hexachlorobiphenyl	0.3	0.38	0.39
7	2,3,3′,4′,6-pentachlorobiphenyl	0.3	0.38	0.39
8	2,2′,4,4′,5,5′-hexachlorobiphenyl	0.3	0.35	0.38
9	2,4,4′,5-tetrachlorobiphenyl	0.3	0.35	0.31
10	2,2′,3,4,4′,5′,6-heptachlorobiphenyl	0.3	0.39	0.40
11	2,3,4′,6-tetrachlorobiphenyl	0.4	0.44	0.42
12	2,2′,3,4,5′-pentachlorobiphenyl	0.4	0.34	0.32
13	2,2′,3,3′,5,6′-hexachlorobiphenyl	0.4	0.38	0.39
14	2,2′,3,3′,4,4′,6-heptachlorobiphenyl	0.4	0.41	0.41
15	2,2′,3,3′,4,5,6′-heptachlorobiphenyl	0.4	0.39	0.40
16	2,2′,3,3′,4,5′,6-heptachlorobiphenyl	0.4	0.41	0.41
17	2,2′,3,4,4′,5,6′-heptachlorobiphenyl	0.4	0.38	0.39
18	2,2′,3,4′,5,5′,6-heptachlorobiphenyl	0.4	0.35	0.39
19	2,2′,3,3′,4,4′,5,5′-octachlorobiphenyl	0.4	0.43	0.42
20	2,2′,3,3′,4,5′,6,6′-octachlorobiphenyl	0.4	0.40	0.41
21	2,2′,3,3′,4,4′,5,6,6′-nonachlorobiphenyl	0.4	0.47	0.45
22	2,3,4′,5-tetrachlorobiphenyl	0.4	0.43	0.39
23	2,2′,3,4,4′,5-hexachlorobiphenyl	0.4	0.34	0.37
24	2,3,3′,4,4′,5-hexachlorobiphenyl	0.4	0.34	0.37
25	2,2′,3,3′,4,4′,5,6-octachlorobiphenyl	0.4	0.40	0.41
26	2,2′,3,3′,4,5,5′,6-octachlorobiphenyl	0.4	0.40	0.41
27	2,3,3′,4-tetrachlorobiphenyl	0.4	0.44	0.39
28	2,2′,3,5′-tetrachlorobiphenyl	0.4	0.43	0.40
29	2,2′,3,4′-tetrachlorobiphenyl	0.5	0.54	0.51
30	2,3,3′,6-tetrachlorobiphenyl	0.5	0.44	0.42
31	2,2′,3,3′,6,6′-hexachlorobiphenyl	0.5	0.37	0.39
32	2,2′,3,3′,4,4′,5,5′,6-nonachlorobiphenyl	0.5	0.47	0.45
33	2,2′,3,3′,4,4′,5,5′,6,6′-decachlorobiphenyl	0.5	0.55	0.49
34	2,2′,4,6′-tetrachlorobiphenyl	0.5	0.40	0.47
35	2,3′,4,4′,5,5′-hexachlorobiphenyl	0.5	0.32	0.32
36	2,4′,6-trichlorobiphenyl	0.6	0.69	0.72
37	2,2′,3,3′-tetrachlorobiphenyl	0.6	0.77	0.61
38	2,2′,3,6-tetrachlorobiphenyl	0.6	0.52	0.53
39	2,3′,4-trichlorobiphenyl	0.7	0.62	0.71
40	2,2′,6-trichlorobiphenyl	0.8	0.77	0.79
41	2,3,4′-trichlorobiphenyl	1.1	0.93	0.97
42	2,3′-dichlorobiphenyl	1.3	1.28	1.31
Test Set
43	2,3′,4′,5-tetrachlorobiphenyl	0.3	0.44	0.40
44	2,2′,4,5,5′-pentachlorobiphenyl	0.3	0.33	0.33
45	2,2′,3,3′,4,4′-hexachlorobiphenyl	0.3	0.35	0.33
46	2,2′,3,4′,6-pentachlorobiphenyl	0.3	0.38	0.39
47	2,3,3′,4,4′,6-hexachlorobiphenyl	0.3	0.37	0.39
48	2,3,3′,4′,5,6-hexachlorobiphenyl	0.3	0.31	0.37
49	2,2′,3,4′,5′-pentachlorobiphenyl	0.4	0.30	0.30
50	2,2′,3,3′,4,5,5′-heptachlorobiphenyl	0.4	0.39	0.40
51	2,2′,3,3′,4,5′,6′-heptachlorobiphenyl	0.4	0.35	0.39
52	2,3,3′,4,4′,5,5′-heptachlorobiphenyl	0.4	0.38	0.39
53	2,3,3′,4,4′,5,5′,6-octachlorobiphenyl	0.4	0.40	0.41
54	2,3,4,4′,5-pentachlorobiphenyl	0.4	0.34	0.37
55	2,2′,3,4,5,5′,6-heptachlorobiphenyl	0.4	0.35	0.39
56	2,2′,3,5-tetrachlorobiphenyl	0.4	0.47	0.45
57	2,2′,3,4,5′,6-hexachlorobiphenyl	0.4	0.38	0.39
58	2,2′,3,3′,4,6′-hexachlorobiphenyl	0.5	0.38	0.39
59	2,2′,3,3′,4,5,5′,6,6′-nonachlorobiphenyl	0.5	0.47	0.45
60	2,3′,5-trichlorobiphenyl	0.5	0.57	0.49
61	2,2′,3,6′-tetrachlorobiphenyl	0.6	0.57	0.54
62	2,2′,5-trichlorobiphenyl	0.7	0.74	0.69
63	2,2′,3-trichlorobiphenyl	1.0	0.96	0.97

Table 2

Characteristics of Molecular Descriptors in Equation

descriptor	coefficient	standard error	t-test	Sig.	VIF
constant	8.181	1.301	6.290	0.000
LOC	–6.199	0.943	–6.572	0.000	1.803
SpMax3_Bh(m)	–0.845	0.279	–3.024	0.005	4.317
SpMin8_Bh(p)	4.447	1.510	2.946	0.006	4.045
SpMin5_Bh(s)	–0.884	0.417	–2.123	0.041	2.533

Table 3

Descriptor Definitions in Equation

model	name	descriptor	block	sub-block
kd	LOC	lopping centric index	topological indices	centric indices
	SpMax3_Bh(m)	largest eigenvalue n. 3 of burden matrix weighted by mass	burden eigenvalues	smallest eigenvalues
	SpMin8_Bh(p)	smallest eigenvalue n. 8 of burden matrix weighted by polarizability	burden eigenvalues	smallest eigenvalues
	SpMin5_Bh(s)	smallest eigenvalue n. 5 of burden matrix weighted by I-state	burden eigenvalues	smallest eigenvalues

An analysis of the depuration rate constants with respect to 285 descriptors (including the average quadrupole moment (Aqm) converted to exp(0.08Aqm)) resulted in the following regressionSimilarly, eq was used for the prediction of the test set. The rms errors for the training set, test set, and total set are 0.0681, 0.0703, and 0.0688, respectively. The prediction rms errors from eq are slightly higher than those of eq . Therefore, eq was taken as the optimal MLR model in this paper.

The four molecular descriptors, LOC, SpMax3_Bh(m), SpMin8_Bh(p), and SpMin5_Bh(s), were used as input variables to develop SVM models for the depuration rate constants (kd) for 63 PCB congeners (see Table ) by applying LibSVM package[17] on the platform of MATLAB R2014a. The regularization constant C and the Gaussian function parameter γ were tuned with the particle swarm optimization (PSO) algorithm under the following experimental conditions: the cognition learning factor c1 (1.5); the social learning factor c2 (1.7); the number of particles N1 (20); the maximum number of iterations N2 (200); the searching ranges of both parameters C and γ being [1, 100]; the N of the N-fold-cross-validation being 5; and ε is the ε-insensitive loss function with a default value of 0.001.

The optimization results from the training set show that the optimal SVM model based on four molecular descriptors in the MLR model possesses parameters C of 10.874 and γ of 28.010. The correlation coefficient R and rms error from the training set are 0.957 and 0.0579, respectively. Subsequently, 21 PCB congeners in the test set were used to evaluate the SVM model. The correlation coefficient R and rms error from the test set are 0.940 and 0.0572, respectively. The R and rms error from the total set are 0.953 and 0.0576, respectively. Obviously, the results from the SVM model are more accurate than that from the MLR model (eq ). Compared with the model[12] (R = 0.853) that deals with depuration rate constants (kd) at 285 K for 62 PCB congeners, both the MLR (R = 0.937) and the SVM model (R = 0.957) in this paper show acceptable statistical results.

As can be seen from Table , the four descriptors (LOC, SpMax3_Bh(m), SpMin8_Bh(p), and SpMin5_Bh(s)) have a sig.-value (or P-value) less than 0.05, which suggest that all the descriptors are significant and related to depuration rate constants of PCBs. In addition, these descriptors possess a variance inflation factor (VIF) of less than 10, which are acceptable, since the correlation coefficients (R2) between any two descriptors are less than 0.9 according to the definition of VIF = 1/(1 – R2). Thus, the four descriptors are weakly correlated with each other.

The t-test can be used to evaluate the significance of descriptors in MLR equations. A descriptor with a larger absolute value of t-test is a more significant descriptor and makes more contribution to the dependent variable. According to the t-test in Table , the most significant descriptor appearing in eq is LOC, which is a lopping centric index. It is derived from a molecular graph with H-depletion and obtained by graph pruning, which is a process of gradually removing all terminal vertices.[15] This descriptor can reflect the degree of molecular branching: the higher the value of LOC, the more branched the graph is. For PCB congeners with similar structures, the PCB compound with more chlorine atoms tends to be inert to degradation and decreases the depuration rate constants (kd). Thus, LOC is negatively correlated with kd.

The three molecular descriptors, SpMax3_Bh(m), SpMin8_Bh(p), and SpMin5_Bh(s), belong to Burden eigenvalues. They were calculated from molecular diagrams containing hydrogen atoms by means of the Burden matrix, Bh(w). This matrix in nature is a modified adjacency matrix and calculated by the following rules: the diagonal elements are atomic properties scaled by carbon atoms; the off-diagonal elements describing the characteristics pairs of bonded atoms are equal to the square root of the conventional bond order; the elements in the matrix that reflects terminal bonds are augmented by 0.1; all other entries are set to 0.001. For each matrix in the Dragon calculation, there are first eight largest positive eigenvalues SpMaxk_Bh(w) and first eight smallest negative eigenvalues SpMink_Bh(w) (absolute values). Here, k is the eigenvalue rank and w is the atomic property such as atomic mass (m), atom polarizability (p), and intrinsic state (s). These molecular descriptors can be used to solve the problem of searching for molecular structural similarity (or diversity).[15]

Figure shows the relationship between standardized residuals (SR) and leverages and can be used to visualize the applicability domain of the SVM model.[14,18] Predictions only for the samples that fall into this domain are considered reliable. As can been seen from Figure , there are three PCBs (no. 42, 2,3′-dichlorobiphenyl; no. 37, 2,2′,3,3′-tetrachlorobiphenyl; and no. 40, 2,2′,6-trichlorobiphenyl) in the training set that possess leverage values greater than the warning leverage h* (=0.2857). But their SR values are close to zero. Thus, the three PCBs can stabilize the SVM model and make it more accurate. In addition, there is only one PCB (no. 62, 2,2′,5-trichlorobiphenyl) having h > h* and SR < 3 in the test set, which indicates that the SVM model has good generalizability.

Figure 1

Plot of standardized residuals versus leverages.

The SVM model was further evaluated with the statistical parameters (qint2, r2, k, k′, R02, R′02, and rm2).[18−23] The value of the internal correlation coefficient qint2 obtained from the training set with leave-one-out cross-validation procedure is 0.690 (>0.5); the value of the correlation coefficient r2 for the test set is 0.883 (>0.6); the slopes k and k′ of the regression lines through the origin are 0.983 and 1.003, respectively, within the range of 0.85–1.15; the values of the determination coefficients R02 (0.845) and R′02 (0.875) are close to those of r2 (0.883) with (r2 – R02)/r2 = 0.044 and (r2 – R′02)/r2 = 0.009, both less than 0.1; and the value of the correlation coefficient r02 (without an intercept for r2) is 0.878 and that of the rm2 = r2 × [1 – (r2 – r02)1/2] is 0.820 (>0.5). These statistical parameters satisfy the criteria defined for an acceptable QSPR model.[19−21] Therefore, the SVM model in this paper is successful.

Conclusions

Four molecular descriptors were successfully used to develop the QSAR models for depuration rate constants (kd) of 63 PCB congeners in juvenile rainbow trout. IEF-PCM and B3LYP/6-31G(d) calculation, with the default solvent (water) and solvent parameters, was used for descriptor derivation. The branching degree (LOC), molecular similarity indices based on atomic mass (SpMax3_Bh(m)), atom polarizability (SpMin8_Bh(p)), and intrinsic state (SpMin5_Bh(s)) can reflect the structure factors affecting the depuration rate constants (kd) of PCBs. Both the MLR and SVM models were tested to be accurate and acceptable with correlation coefficients of 0.933 and 0.953, respectively. It is feasible to develop QSAR models for depuration rate constants (kd) of PCB congeners by applying IEF-PCM and B3LYP/6-31G(d) calculations for molecular descriptor derivation. This result encourages the further application of QSAR to other properties, such as photodegradation half-life and biological activities of PCBs.

Materials and Methods

Experimental Data

The experimental values of depuration rate constants (kd) for 63 PCB congeners in juvenile rainbow trout (O. mykiss) were taken from the literature[16] and are listed in Table . Sixty-three PCB congeners were randomly split into a training set (42 PCBs) and a test set (21 PCBs). The former was used to train QSAR models, and the latter was used to validate the models. In the experiment,[16] these fishes were exposed to a high-dose treatment of PCBs and the kd values measured at 281 K ranged from 0.3 to 1.3. Because the PCB congeners possess similar molecular structures, the range of the kd values is relatively narrow. This phenomenon can be found in other properties.[12,13] Therefore, the QSAR models based on PCBs can only be used for the prediction of the properties of PCBs.

Structural Descriptor Calculations

The molecular structures of PCB congeners were sketched by ChemBioDraw Ultra 11.0 and preliminarily optimized by ChemBio3D Ultra 11.0 in ChemOffice 2008 with molecular mechanics (MM2) method with the default convergence criterion. The energy-minimized molecules were further optimized by applying IEF-PCM/B3LYP/6-31G(d) in Gaussian 09 program, with the default solvent (water) and solvent parameters. There are two reasons for carrying out the IEF-PCM and B3LYP/6-31G(d) calculations for the molecular descriptor derivation of PCB congeners. One is that this approach can be used to optimize the molecular structures of solute at a reasonable level, which has been investigated by our previous work.[14] The other is the fact that this approach can provide more structural information and molecular descriptors than other semiempirical quantum chemical calculations such as the AM1 and PM6. For example, IEF-PCM and B3LYP/6-31G(d) can be used for the calculation of the average quadrupole moment (Aqm), although these descriptors may not be present in models.

Besides the average quadrupole moment, Dragon 6.0 was used for the molecular descriptor derivation after full geometry optimization. Redundant and nonuseful structure descriptors that are equal to a constant or have a pair of correlation coefficients ≥0.90 were deleted from the 4885 descriptors calculated. In total, 284 descriptors were obtained with Dragon 6.0.[15]

SVM Principle

The SVM algorithm has many attractive characteristics and advantages. Unlike the traditional theory, the SVM adopts the principle of structural risk minimization, which makes it possible for the SVM to have a good prediction ability even with fewer training samples. The SVM was originally used in pattern recognition. With the introduction of insensitive loss function, the application of SVM extends to the prediction of properties in nonlinear regression models and the evaluation in a time-series analysis.

The main idea of the SVM used in the regression prediction is to map the input characteristic parameters into a high-dimensional feature space through a nonlinear mapping function and then carry out a linear regression analysis. The SVM algorithm is based on the following regression[2,4,24]where n is the number of samples used in the training set, φ(x) is the mapping function, x is the input vector composed of characteristic parameters reflecting molecular structure, and f(x) is the output of prediction results. The coefficients w and b can be estimated by means of the following minimizationsubject towhere C is a penalty parameter used for adjusting the training error. ε is a prescribed parameter of the ε-insensitive loss function. ξ and ξ* are the slack variables. In general, support vector regression is based on the ε-insensitive loss function to minimize the training errorThus, eq becomeswhere α and α* are the Lagrange multipliers that are used to solve the quadratic optimization problem. Equation can be obtained when the kernel function k(x, y) is introducedwhere s is equal to the number of samples whose α and α* values are above 0. In the present work, the following Gaussian radial basis function was usedThe kernel width γ affects the prediction performance of the SVM models. In this paper, the SVM parameters C and γ are tuned with the particle swarm optimization (PSO) algorithm.[25−28]