Biological enrichment prediction of polychlorinated biphenyls and novel molecular design based on 3D-QSAR/HQSAR associated with molecule docking

Polychlorinated biphenyls (PCBs for short) are the general names of compounds based on the substitution of a hydrogen atom on a benzene ring with a chlorine atom, including 209 congeners characterized by the number and position of the chlorine atoms on the biphenyl core. PCBs are one of the persistent organic pollutants (POPs); PCBs are high in stability, toxicity, environmental persistence, bioaccumulation, long-distance migration ability, and other characteristics. PCBs have the advantages of low solubility, high dielectric constant, low vapor pressure, good heat resistance, and excellent insulating properties. PCBs widely served as flame retardant, paint, dielectric fluids in capacitors, insulating oil in transformers, and as a plasticizer agent and rubber sealants [1].

In the 1970s, PCBs were prohibited from being produced and used worldwide [2]. However, because of the high stability of PCBs and PCBs themselves being resistant to their natural degradation process, the discontinuation does not mean that the hidden dangers are eliminated, and the destructiveness of PCBs can persist for several years. Today, 40 years after the PCBs have been banned globally, PCB residues are still universally detected in organisms and various environmental media in the world, affecting various organisms including humans through the cumulative effect of the food chain. The hazards to environment posed by PCBs are longer term and more complex and have become global environmental issues that pose a serious threat to human health and the ecological environment. Therefore, studies on PCBs are of great significance in the environmental pollution and biological health aspects and have received substantially more global attention.

3D-QSAR, which could reveal the structure–activity relationship of compounds by analyzing the changes of 3D spatial fields around the superposed molecules, overcomes the limitations of the conventional 2D model in characterizing the relationship between property and substance and have a clearer physical meaning and more abundant information of the molecular field energy, thus 3D-QSAR is more widely used in environmental science and toxicology and other fields than 2D-QSAR [3]. Some studies have also used some kind of bacteria to establish the method of toxicity test of some chemicals, and obtained good research results [4,5]. Su et al. [6] determined the combined toxicity of phenolic compounds and heavy metal lead to the photobacterium and then conducted a QSAR study. In addition to the use of 3D-QSAR method, this dissertation also combined hologram quantitative structure–activity relat

ionship (HQSAR) method together to explore the mechanisms that affect the activity of PCB compounds.

In recent years, various methods have been used to estimate the biological enrichment coefficients of PCBs in various mediums, and to assess the biological enrichment degree of PCBs in sediments and various biological mediums. The concentration of PCBs detected in main species in the food chain of the Antarctic and Arctic has reached moderate levels of contamination [7]. Engel et al. [8] used logistic regression analysis to find that some PCBs in serum were associated with an increased risk of lymphoma. Experts in environmental field at home and abroad have also conducted a great deal of research on the migration and conversion behavior, ecological effects, toxicity as well as pollution elimination methods of PCBs in water bodies, sediments, soils, and organisms [9–14], including microbial degradation and bioremediation methods of PCBs [15]. At present, the researches on PCBs in terms of toxicology, exposure characteristics, adsorption, and degradation have been relatively comprehensive, but there is a lack of research on the modification of the POPs characteristics at the molecular structure level.

Although PCBs have been banned, but there are still residues on the global scale, and the modification of PCBs in the present paper was based on the residual PCBs. PCBs have the advantages of various commercial functional properties, and modification of PCB molecules makes sense if the POP properties of newly designed PCB molecules are reduced or even removed from POPs without changing their functional properties.

In the present paper, HQSAR and 3D-QSAR models were established by using with 13 PCB octanol-water partition coefficients (K_ow) from Hawker and Connell [16]. The substitution sites and substituent groups that affect the octanol-water partition coefficients of PCBs were determined by combining these two models. And carry out a modification of target molecular PCB-207 to lower its K_ow (biological enrichment). Only using the 3D contour plot of 3D-QSAR to determine the substitution sites and substituent groups, the fewer sites and groups can be obtained, and almost impossible to do the double-site substitution. Less scheme on molecular modification, it is not enough to provide reference for future research. The octanol-water partition coefficient (K_ow), which can reflect the distribution ability of organic compounds between octanol and aqueous phases, is one of the effective index of the distribution of organic compound in environmental media (water, soil, and sediment) and is also one of the important property parameters used to study the environmental behavior of organic compounds. Additionally, K_ow can simulate the distribution of organic between biological phase and water phase, and it is closely related to the toxicity [17], biological enrichment [18], and solubility [19] of compounds.

Through the evaluation of toxicity, persistence, and long-range mobility, the low biological enrichment newly designed PCB-207 molecules, all the POP characteristic parameters are reduced while the actual functional properties are not changed, fulfilling their industrial and commercial functional requirements and greatly reducing the environmental impact. Furthermore, the correlation analysis between the number of Cl atoms in PCBs and the average value of logK_ow of isomer in chlorobiphenyl to decachlorobiphenyl, respectively, showed that with the increase in the number of Cl atoms, the logK_ow values also increase and the biological enrichment capability becomes stronger, and the results of molecular docking between PCBs and degrading enzyme BphA also confirmed this conclusion.

The target pollutant PCB-207 selected in the present paper has large observed logK_ow value of 7.52. After modification, the logK_ow values of newly designed molecules decreased and the lowest reached 5.489. When logK_ow <5, the compounds were judged to be divorced from POPs, if the modification is carried out by using PCBs with less chlorine atoms as target molecular, which the logK_ow values of the compounds are relatively small, the logK_ow values can also be reduced by the same extent, then that of modified compounds must be less than 5, which can make the newly designed compounds out of the scope of POPs. In the present paper, we mainly come up with a modified method for the residual PCBs range from chlorobiphenyl to decachlorobiphenyl, which provides a theoretical basis in methods for the conversion of residual PCBs into low environmental impact compounds. In addition, the molecular docking technique was used to further study the biological enrichment of the modified molecules on the liver tissue where the liver enzymes exist in, and it was found that not all of bioaccumulation of newly designed molecules of low biological enrichment of PCB-207 was reduced on liver tissue. But it can be screened out that the several kinds of low biological enrichment newly designed PCB-207 molecules from all newly designed PCB-207 molecules have lower biological enrichment in liver tissue where most liver enzymes exist in, and provide an important theoretical method for the study of bioaccumulation of PCBs in the future. The present paper is expected to provide a theoretical foundation for further study of biological enrichment of PCBs.

The comparative molecular field analysis (CoMFA), comparative molecular similarity indices analysis (CoMSIA) models, HQSAR model that was used for K_ow value prediction and low biological enrichment molecular design were established with Sybyl-X 2.0 software to analyze the 3D/HQSAR models and perform the molecular docking.

The lowest energy conformation of the molecule was used as the advantage stable conformation; the geometry of these compounds was subsequently optimized using Tripos force field with the Gasteiger-Huckel charges. Repeated minimizations were performed by the Powell method with a maximum iteration of 10000 to reach an energy convergence gradient value of 0.005 kJ/mol, and others are default values.

To establish the 3D-QSAR models, the whole dataset (containing 13 compounds) was divided 3:1 into a training set (containing ten compounds) for 3D-QSAR model generation and a test set (containing three compounds) for model validation for its accuracy and stability. Compounds in the training set were selected based on their ability to appropriately represent the structural diversity of the whole dataset and cover the range of logK_ow values [20,21]. In this paper, we have chosen a set of training set and test set that both internal and external verifications were qualified to establish 3D-QSAR model. In this process, the 2,2′,3,3′,4,4′,5,5′-octachlorobiphenyl (PCB-194) molecule with the largest logK_ow was used as the template and the benzene ring skeleton was used as a model skeleton to align the other compounds using the Align Database command in Sybyl. The 3D-QSAR of CoMFA and CoMSIA was directly yielded by partial least squares (PLS) analyses in which the ten logK_ow values of PCBs in training set served as the dependent variable, and the 3D structure of ten PCBs served as the independent variable. In this process, the leave-one-out (LOO) cross-validation procedure was performed to determine the optimum number of components (n) and the highest cross-validation correlation coefficient (q²) for the correlation models. Simultaneously, noncross-validated analysis was performed. The quality of the models was measured by cross-validation correlation coefficient (q²), noncross-validation correlation coefficient (r²), standard error of estimate (SEE), and Fisher test (F) values. Finally, the relationship between the logK_ow of the PCB compounds and each field was expressed as a 3D contour plot.

HQSAR model does not require the selection of active conformations and the alignment of PCB molecules, HQSAR and PLS analysis was conducted using the Sybyl-X 2.0 package by Tripos Company. HQSAR calculations were carried out using three distinct parameters the fragment distinction, the fragment size, and the hologram length. Use SYBYL-HQSAR module to generate molecular hologram. The HQSAR calculation provides 12 default prime number (53, 59, 61, 71, 83, 97, 151, 199, 257, 307, 353, and 401) molecular hologram lengths by default. The HQSAR model can be optimized by changing the fragment distinction and the molecular fragment size. The fragment discrimination represents the topology parameter mapped in the molecular holography, which included atoms (A), bonds (B), connections (C), hydrogen atoms (H), chirality (Ch), and donor and acceptor (DA), in which atoms (A) can distinguish the different types of atoms; the bonds (B) can identify the difference amongst the chemical bonds formed by the atoms; connections (C) can reveal the hybridization state of atoms inside the fragment; Ch can obtain stereochemical information of atomic Ch and chemical bonds in fragments; hydrogen bond DA can ascertain hydrogen bond donor or acceptor of fragments. In general, the default setting A/B/C for the fragment distinction has already contained the basic information required to distinguish different fragments. On this basis, adding other fragment distinction parameters, establishing different combinations of fragment distinction parameters, selecting different fragment sizes. And the HQSAR models were obtained in specific preset parameters. The fragment size is the number of atoms contained in the fragment: (1–3) means the smaller atomic fragments, which can characterize the atomic type and functional groups, approximately (4–7) means medium atomic fragments that could distinguish the chain length of the alicyclic hydrocarbon or the characters and the structural features such as the substitution sites and substituent groups of the aromatic hydrocarbon; (8–10) means larger atomic fragments.

HQSAR models (method) were conducted according to the different parameter settings, including the hologram length (HL) values (53, 59, 61, 71, 83, 97, 151, 199, 257, 307, 353, and 401), and the fragment size (4–7) by default. Based on the default setting A/B/C of the fragment distinction, eight kinds of fragment distinction combinations are constituted after adding other fragment distinction parameters, and the combination modes are A/B/C, A/B/C/H, A/B/C/Ch, A/B/C/DA, A/B/C/H/Ch, A/B/C/H/DA, A/B/C/Ch/DA, and A/B/C/H/Ch/DA. Molecular hologram of 13 training set molecules and activity were used as the independent and dependent variables, respectively. A series of HQSAR models were obtained by taking linear regression analysis through PLS. Then, the prediction ability, stability, and fitting capacity of these models were validated by a high cross-validated correlation coefficient q², noncross-validated correlation coefficient r² values, and a low-SEE along with a low cross-validated SEE (SE_cv). The HQSAR analysis results can be graphically displayed in the form of contribution maps, in which the color coding of each atom reflects its contribution to the activity of the whole molecule. Hence, we can obtain the favorable information of molecular modification [22].

The substitution sites identified by HQSAR models that have a significant effect on the biological enrichment activity of the target molecular, combined with the distribution characteristics of the force field (hydrophobic field and hydrophobic group) that has a significant effect on the bioaccumulation of the target molecular determined by CoMSIA model of 3D-QSAR. It can precisely locate the low biological enrichment activity modified substitution sites and the substituent groups of target molecular, and then modify the target molecules.

In the present paper, the molecular docking was performed by using SYBYL-X 2.0 software from Tripos Company in the United States. The molecular docking of PCB-207/newly designed PCB-207 that have low biological enrichement with target protein of liver enzymes were both the docking that between small molecular and protein. Therefore, use Surflex-dock module by SYBYL software for the semiflexible docking [23,24]. Open the Surflex-dock module to correct and modify the protein that is about to carry on molecular docking, and extract the ligand from it to determine the binding sites. Set up docking pockets for the prepared protein, put the ligand molecules that are about to carry on molecular docking into the docking pocket, prepare for the molecular docking to take place. The evaluation standards of molecular docking results include Crash, Polar, and Total Score three expression functions. Crash represents the inadaptation of docking for the ligand into the receptor, the closer to 0 the better, that is the smaller the absolute value the better; Polar stands for the polarity function score, which the higher the scores the better when the binding molecule sites on the binding surface while the lower the scores the better when the binding site is inside the molecule; Total Score denotes the comprehensive score of above parameters, the higher the score the better. According to Total Score, it could be to determine the appropriateness of docking between molecules and proteins.

Quantum chemical descriptor of PCB-207 molecules and low biological enrichment newly designed PCB-207 molecules were calculated using Gaussian09, including energy gap (insulativity parameter) and C–Cl bond dissociation enthalpy (flame retardancy parameter). The geometry of all compounds was optimized using B3LYP/6-31g (d, p) level of density functional theory [25].

The calculation orbit (highest occupied molecular orbital (HOMO) and lowest unoccupied molecular orbital (LUMO)), which is significant to the study of molecular physics and chemical properties, is an important parameter of quantum chemistry. The energy gap, which is the energy difference of the HOMO and LUMO, is the required lowest energy of the electron excitation process, and is also an important parameter of conductivity and luminescence [25]. In addition, the larger the energy gap, the weaker the conductivity. The insulation information of newly designed molecules can be obtained by calculating the energy gap value.

By comparing the C–Cl bond dissociation enthalpy of PCB-207 molecules and low biological enrichment newly designed PCB-207 molecules, we obtained the flame retardancy information about low biological enrichment newly designed PCB-207 molecules. The lower the C–Cl bond dissociation enthalpy, indicating that the newly designed PCB-207 molecules after modification are more likely to release Cl-free radicals and HCl for inflaming retarding.

The CoMFA in the 3D-QSAR model mainly reflects the non-bonding interaction between the drug molecules and the receptor, and is a widely used indirect drug design method. The steric fields’ energy and electrostatic fields’ energy are calculated using the Lennard-Jones and Coul potential energy functions. However, the lattice points of molecular surfaces are ignored due to the rapid increase in van der Waals repulsion. Compared with the CoMFA method, the CoMSIA method uses Gaussian similarity functions, which avoids drastic changes of the potential energy at the lattice points on the molecular surface, improves the sensitivity to the molecular alignment method and the spatial orientation, and increases hydrophobic fields in addition to hydrogen bond-donor and hydrogen bond-acceptor fields. The results of the CoMSIA method show a relatively small influence from compound matching rules and can more intuitively explain the quantitative structure–activity relationship of a compound. The use of CoMSIA can overcome inherent defects of CoMFA but does not necessarily gain better results [26]. Therefore, in this study, two methods (CoMFA and CoMSIA) were used to verify and supplement each other to gain reliable predicting models.

The statistical parameters of the CoMFA model and the CoMSIA model were shown in Table 1. According to Table 1, the CoMFA model and the CoMSIA model both had an optimum n of 2, the cross-validated q² of 0.784, 0.883 (>0.5), the SEE of 0.215, 0.177, the noncross-validated r² of 0.954, 0.969 (>0.8), the F-value of 72.595, 108.026, and the scrambling stability test parameters Q² of 0.715, 0.787, cSDEP of 0.525, 0.452, and dQ²/dr²yy of 1.625, 1.309, respectively, illustrating that these two models had both suitable fitting predictive abilities [27,28].

As is shown in Table 1, the corresponding contribution percentages of the steric (S) and electrostatic (E) fields to the CoMFA model were 33.90 and 66.10%, respectively, indicating that the steric and electrostatic distributions of the groups may influence the K_ow values of the PCB homologs. Verifying electrostatic interaction was the major contribution to the logK_ow of PCBs. The CoMSIA defines explicit hydrophobic and hydrogen-bond DA descriptors in addition to the steric and electrostatic fields in CoMFA. The contributions of the steric (S), electrostatic (E), hydrophobic (H), hydrogen bond-donor (D), and hydrogen bond-acceptor (A) fields to the CoMSIA model were 0.90, 84.40, 14.70, 0.0, and 0.00%, respectively.

Take above discussion into consideration, the steric distributions, electrostatic distributions, and hydrophobic properties of the groups may influence the K_ow values of the PCBs homologs, the electrostatic distribution had the most influence, and the hydrogen bond-donor and hydrogen bond-acceptor fields had no effect.

External validation was also conducted to further assess the reliabilities and the predictive ability of the built models. The test set with three compounds was used for this validation. The r²_pred of 0.896 (>0.6), 0.937 (>0.6), the SEP of 0.301, 0.234, the Q²_ext of 0.996 (>0.5), 0.891 (>0.5) were achieved, verifying the good external predictive ability of the two models [28]. The correlations between the experimental logK_ow values and logK_ow values predicted by the CoMFA and CoMSIA models are depicted in Figure 1A,B, respectively. In further analysis of the scatter plot of experimental versus logK_ow values predicted by the CoMFA and CoMSIA models revealed a fine linear dependence (r² of 0.9579 and 0.9641, respectively) among experimental values and predicted values, all data were concentrated around the trend line. And the slopes of the linear equations for experimental and predicted values were 0.96105 and 0.96343, respectively, which can be observed from the scatter plot of experimental versus predicted values of logK_ow, which can be observed from the scatter plot of experimental versus predicted values of logK_ow, which indicates the two models have good internal prediction ability and can be used for the prediction of PCBs K_ow values.

CoMFA and CoMSIA models were used to predict 209 types of PCBs, and the experimental and predicted logK_ow values and standard deviations for the PCBs are given in Table 2.

From Table 2, in the CoMFA and CoMSIA model prediction, the relative errors of the test set (3 compounds) were 0.12, 2.20, 6.41% and 2.86, 1.72, 7.31%, respectively, which are acceptable for small base values.

Based on the default setting A/B/C of the fragment distinction, eight kinds of fragment distinction combinations (A/B/C, A/B/C/H, A/B/C/Ch, A/B/C/DA, A/B/C/H/Ch, A/B/C/H/DA, A/B/C/Ch/DA, and A/B/C/H/Ch/DA.) are constituted after adding other fragment distinction parameters. The statistical results of the PLS analyses for the training set using different fragment distinctions in combinations with default fragment size (4–7) and 12 default hologram lengths (53, 59, 61, 71, 83, 97, 151, 199, 257, 307, 353, and 401) were summarized in Table 3.

It can be seen from Table 3 that each parameter result of fragment distinction combinations among A/B/C and A/B/C/Ch, A/B/C/DA, A/B/C/Ch/DA was exactly the same, and each parameter result was identical for fragment distinction combinations A/B/C/H and A/B/C/H/Ch, A/B/C/H/DA, A/B/C/H/Ch/DA, indicating that the fragment distinction parameter Ch cannot distinguish the fragment, and it will not make any differences in Ch of molecules with the introduction of fragment distinction parameter DA. That is to say, PCBs molecules do not contain chiral atoms, by the same token, the fragment distinction parameter DA cannot distinguish fragments, and the introduction of fragment distinction parameter Ch did not make the statistical results of model differences, that is, there were no hydrogen bond donor atoms and hydrogen bond acceptor atoms in PCB molecules. This result is the same as that of the CoMFA model, confirming each other.

It was found that the r² of above eight fragment distinction combinations were greater than 0.8, and q² were greater than 0.5. Often, a high value of this statistical characteristic (r² > 0.8, q² > 0.5) is considered as a proof of the high predictive ability of the model. Taking into account the predictive SEE and SE_cv values, A/B/C/H was proved as the most effective fragment distinction combinations to build the HQSAR model. According to HQSAR contribution maps, we can get the favorable information for molecular modification.

As shown in Figure 2B, the order of color about the contribution to activity is: green > green blue > yellow > white > orange > red orange > red. The color from green to red represents the positive contribution (PC) to negative contribution to the activity. As can be seen from Figure 2, the green-labeled Cl₄, Cl₅, Cl₆ molecular Cl-substitutions have a PC to the bioaccumulation activity.

Taking the most active compound PCB-207 as an example, the properties of PCB molecules were analyzed by using 3D contour map of the model. The steric field is denoted by yellow and green colored contours, where the yellow regions represent the small volume groups near these regions favorable to logK_ow, and the green regions indicate that the sterically bulkier substituents close to these regions may increase logK_ow. In the electrostatic field, blue-colored contours represent regions where the positive charge increases the logK_ow values, whereas red-colored regions display areas where the negative charge enhances the logK_ow values. Moreover, the yellow and white contours depict hydrophobic and hydrophilic favored regions, respectively. To visualize the field effects on the target compounds in 3D space, the contour diagrams of the final CoMFA and CoMSIA models are shown in Figure 3.

The 3D contour maps of CoMFA model were seen from Figure 3. In steric fields (Figure 3A), large yellow contours located at the 3, 4, 5-positions of the A ring and 3′, 4′-positions of B ring and indicated that the sterically smaller substituent is favored at these positions. A medium-sized region of yellow contour covered the 6′-position of the B ring, indicating that these positions were preferred for smaller substituents to increase logK_ow, that is, increase its solubility in organic phase, reduce its solubility in water, as illustrated by the fact that the logK_ow value of PCB-207 was stronger than that of compound PCB-201, that is to say, introducing bulky Cl-substituent that at site 5 decreased the solubility of PCBs in the organic phase. Cl₂, Cl₄, Cl₅, Cl₆, Cl₃′, Cl₆′ atoms were located in the blue contours of the electrostatic contour map, and Cl₃, Cl₃′, Cl₄′ atoms were located in the red contours of the electrostatic contour map (Figure 3B), indicating that imbedding electron withdrawing groups at 2, 4, 5, 6, 3′, 6′-positions and imbedding electronegative groups at 3, 3′, 4′-positions can increase the logK_ow of PCBs. The higher logK_ow of PCB-194 compared with PCB-207 was an example of such a case.

As shown in Figure 4, Cl₂, Cl₂′ atoms and carbon skeleton near the Cl₅, Cl₆, Cl₆′ atoms were located in yellow contours of the steric contour map, Cl₃′, Cl₄, atoms were located in green contours of the steric contour map, and the small green colored contours were mapped near the 3- and 5-positions, indicating that small volume groups at 2-, 2′-positions and bulkier groups at 3-, 5-, 3′-, 4′-positions may increase logK_ow. In electrostatic fields, two small blue colored contours were mapped the carbon skeleton near the Cl₄ and Cl₄′ atoms, the 3-, 4-, 5-, 3′-, and 4′-positions are encompassed by medium-sized region of red contours, representing the electronegative groups at 3-, 4-, 5-, 3′-, 4′-Substitution were favorable to increase logK_owof PCBs. According to Figure 4C, two large white contour located at the 3-, 4-, 5-positions and 3′-, 4′-, 5′-positions, respectively, indicated that the sterically bulkier substituent is favored at above positions to increase the logK_ow of PCBs. From the contribution rates of the descriptor fields, the two models mutually verify and prove that electronic effects primarily influence the logK_ow.

Although the models exhibited a satisfactory fitting ability, stability, and predictive ability, only the CoMFA model contained two descriptor fields (steric and electrostatic fields), presenting certain limitations in terms of the analysis of the effect of the descriptor field and the modification of the information on the compounds. The CoMSIA model contained five descriptor fields (steric fields, electrostatic fields, hydrophobic fields, and hydrogen bond donor/acceptor fields), which provided a more comprehensive understanding of the effect on the physical and chemical properties of PCBs. For example, an analysis of the hydrophobic field could reveal that hydrophobic groups on the 3-, 4-, 5-, 3′-, 4′-, and 5′-positions of the benzene ring would increase the logK_ow.

Furthermore, by comparing 3D contour maps of the CoMFA and CoMSIA models, we found that massive regions of the CoMFA contour maps were larger, and the distribution of molecular fields was wide, making it difficult to pinpoint the influence area. Therefore, contour maps of the CoMFA model were selected to modify compounds.

The effect of each field on the substitution characteristics of each site of the target molecular can be obtained by 3D-QSAR model, and the HQSAR model can determine the substitution sites that contribute most to the bioaccumulation of the compounds. Combined with the 3D contour plot of 3D-QSAR CoMSIA model and activity contribution diagram of HQSAR model, the target molecular can be modified to the low biological enrichment newly designed PCBs-207 molecules by determining the substitution sites and substituent groups.

As can be known from target molecular activity contribution diagram of HQSAR, Cl₄, Cl₅, Cl₆ molecular Cl-substitutions have an uppermost PC to the bioaccumulation activity. A synthetic analysis of the effect of steric, electrostatic, and hydrophobic fields in the CoMSIA model showed that Cl₅ molecular Cl-Substitution was located in green contours of the steric contour map; Cl₄, Cl₅ atoms were located in the red contours of the electrostatic contour map; Cl₄, Cl₅ atoms were located in the white contours of the hydrophobic contour map. Know then, Cl₄ molecular Cl-Substitution was affected by the electrostatic fields and hydrophobic fields, Cl₅ atom was affected by three molecular fields of the steric fields, electrostatic fields and hydrophobic fields. According to the PLS analysis in Table 1, the contribution rate of the steric fields under the CoMSIA model is only 0.90%, so only the influence of the electrostatic fields and the hydrophobic fields on the molecular modification is considered. Therefore, it was conducive to reducing the K_ow value of monosubstituted compounds introducing electron withdrawing or hydrophobic groups at Cl₄, Cl₅ two positions, respectively, and reducing the K_ow value of the bis-substituted compounds introducing electron withdrawing or hydrophobic groups at Cl₄, Cl₅ positions at same time.

Based on the above analysis, ten kinds of groups with less electronegativity than Cl atom and hydrophobic were introduced to modify the target molecular (PCB-207) at Cl₄ and Cl₅ Cl-substitutions. A total of 63 monosubstituted and bis-substituted substitution schemes were established by combining substitution sites and substituent groups. According to the Predicted logK_ow values of newly designed PCB-207 molecules through CoMFA and CoMSIA models, the K_ow values of newly designed PCB-207 compounds decreased by embedding the groups with less electronegativity than Cl atom (-OH, -CH₂OH, -NH₂, -NO), and groups hydrophobic besides less electronegativity (-CH₂CH₃, -Br, -CH₃, -NO₂, -Phenyl, -OCH₃), indicating that the introduction of groups is the reason for the descending K_ow values, the other words descending biological enrichment of the modified compounds.

From the above analysis, it is known that the bioaccumulation of newly designed PCB-207 molecules is reduced. However, in addition to the bioaccumulation of the four POPs characteristics, the toxicity, persistence, and long-range mobility of the low biological enrichment newly designed PCB-207 molecules needs to be evaluated and analyzed. A total of 32 low biological enrichment newly designed PCB-207 molecules with a K_ow value reduced by more than 10% were selected from 63 newly designed PCB-207 molecules mentioned above and evaluated for toxicity, persistence and long-range mobility to further screen out the optimally modified compounds, and evaluate its functional properties that is to say the flame retardancy, insulation evaluation.

Li et al. [29] constructed the 3D-QSAR model to predict the toxicity parameter pEC₅₀ of PCBs, Xu et al. [30] constructed the 3D-QSAR model to predict the persistence parameter t_1/2 of PCBs, and Chen et al. [31] constructed the 3D-QSAR model to predict the long-distance mobility parameter K_OA of PCBs, the present paper will use above model to predict the toxicity, long-lasting, and long-distance mobility of newly designed PCB-207 molecules.

Based on the analysis of 32 new molecules, we knew that the logK_ow values of 32 newly designed PCB-207 molecules decreased significantly compared with PCB-194 and PCB-201, indicating that the introduction of groups is the reason for the descending K_ow values of the modified compounds and not the decrease in Cl. Thirty-two kinds of low biological enrichment newly designed PCB-207 molecules with a logK_ow value reduced by more than 10% are mostly obtained after being modified and substituted at two substitution sites at same time, thus the logK_ow reduction in bis-substituted newly designed molecules is greater than that of monosubstituted newly designed compounds. 12-low biological enrichment newly designed PCB-207 molecules with reduced toxicity, persistence, and long-distance mobility were screened from 32 compounds, and these 12 newly designed PCB-207 molecules were calculated using Gaussian09 for their functional properties. (63 substitution schemes of new designed PCB-207 molecules are shown in the supplementary material.)

Energy gap (insulativity parameter) and C–Cl bond dissociation enthalpy (flame retardancy parameter) of PCB-207 molecules and low biological enrichment newly designed PCB-207 molecules were calculated using Gaussian09 showed in Table 4. The insulation information of newly designed molecules can be obtained by calculating the energy gap value, in the same way, the flame retardancy information about low biological enrichment newly designed PCB-207 molecules can be obtained by calculating the C–Cl bond dissociation enthalpy.

As shown in Table 4, the dissociation enthalpies of C–Cl bonds are lower than those of the target molecular PCB-207 (70.706 kcal/mol) under the condition that four POPs characteristic parameters of 12 low biological enrichment newly designed PCB-207 molecules are reduced. Studies have shown that PCBs will release HCl during combustion process, capturing and reacting with highly active H and OH radicals produced during the combustion reaction of polymer materials sequentially, eventually leading to the slowdown or termination of combustion [32,33]. The dissociation of C–Cl bond is closely related to the flame retardant efficiency of PCBs, indicating that the low biological enrichment newly designed PCB-207 molecules after being modified are more likely to release Cl-free radicals and HCl for inflaming retarding. The energy gap values increase in four low biological enrichment new designed PCB-207 molecules (5-methyl-PCB-207, 4-bromine-5-hydroxymethyl-PCB-207, 4-ethyl-5-hydroxymethyl-PCB-207, and 4-methyl-5-hydroxymethyl-PCB-207) was 0.40–1.35%, indicating that the insulation increased; and the energy gap values decreased of the rest eight low biological enrichment newly designed PCB-207 molecules was from −0.37 to 0.20%, the amplitude of increase or decrease both are small, thus explaining that the effect on the degree of energy gap of low biological enrichment modified PCB-207 compounds improves the flame retardancy while the insulation remained essentially unchanged.

There are ten Cl isomers of PCBs, one to ten chlorinated biphenyls. The values of 209 kinds of PCBs molecules in Table 2 are used for predicting the mean value of logK_ow of PCBs from one to ten chlorine number. The values of fitting correlation coefficients R² for the value of logK_ow and the number of chlorine atom under CoMFA and CoMSIA model are 0.70911 and 0.81723 (n=10, P=0.01, r > r₀), respectively. The values of R² for the mean value of logK_ow and the number of chlorine atom under each category are 0.9907 and 0.9992 (n=209, P=0.01, r > r₀), respectively. These results show that the values of logK_ow gradually increased with the increasing number of Cl atoms, and stronger K_ow values of PCBs imply a stronger biological enrichment ability of PCBs. Thus, the biological enrichment increased with the increasing number of Cl atoms.

The names of the five most common PCB mixtures commodity are Aroclor1016 (41.5% Cl), Aroclor1242 (42% Cl), Aroclor1248 (48% Cl), Aroclor1254 (54% Cl), and Aroclor1260 (60% Cl), respectively. Except for Aroclor1016, the first two digits of other commodity name are 12, which represent the number of carbon atoms in PCB molecules. The last two digits represent the isomer ratios of Cl in the commercial mixtures of PCBs, indicating the name of commodity according to the different degree of chloride.

Table 5 [31] presents the ratio of the isomers contained in each PCBs commodity mixtures, which show that the Cl isomers from trichlorodiphenyl to heptachlorobiphenyl make up most of the commercial products of PCBs, and chlorobiphenyl and dichlorobiphenyl also account for a certain percentage of the total products, while the Cl range from octachlorobiphenyl to decachlorobiphenyl is rarely used in commercial mixture of PCBs.

The following PCB congeners were detected in Lampetra fluviatilis in Europe by Merivirta et al. [34]: PCB-28, 52, 101, 114, 118, 128, 138, 141, 149, 151, 153, 170, 180. The range of these congeners is from trichlorodiphenyl to heptachlorobiphenyl. The study of PCB accumulation and dietary exposure to fish in Hyderabad, India, by Ahmed et al. [35], found that PCB congeners in all fish species are distributed in the dichlorobiphenyl to decachlorobiphenyl. Russo et al. [36] traced PCBs enriched in water samples found that PCB homologs recovered from water samples contained (PCB-1, 15, 31, 44, 138, 180, 195), in which the range is chlorobiphenyl, dichlorobiphenyl, trichlorodiphenyl, tetrachlorobiphenyl, hexachlorobiphenyl, heptachlorobiphenyl, and octachlorobiphenyl. This is consistent with the finding in Table 5 that the range of almost all PCBs that can be detected in all aquatic organisms cover trichlorodiphenyl to heptachlorobiphenyl. The results showed that the more the certain PCBs homolog used, the more the detection.

Take 209 kinds of PCBs bioaccumulation in fish as the target molecule docking with its degrading enzyme, getting the corresponding scoring function associated with its bioaccumulation. The higher the scoring function, the better the PCBs and degrading enzymes bind, and the easier degradable PCB molecules are, the lower the bioaccumulation is. Therefore, the relationship between the number of Cl atoms and the bioaccumulation is verified by this way. BphA is the only enzyme direct contact with PCBs in biodegradation. Thus, BphA was used to dock with 209 PCBs using protein structure of 2GBX [37].

According to the scoring functions of all the Cl isomers from chlorobiphenyl to decachlorobiphenyl binding with BphA enzyme, the average Total Score under each Cl atom classification was 3.087, 2.564, 2.187, 1.560, 1.177, 0.817, 0.633, 0.432, 0.427, and 0.130, respectively, the average Total Score under each Cl atom classification of PCBs binding with BphA enzyme and the number of Cl atoms fit the curve y= −0.32417x + 3.08427, R² = 0.93251, thus r = 0.9656. When the significance level P=0.01, r₀ = 0.7645, there was a significant linear relationship between the total score of PCBs binding to BphA enzyme and the number of Cl atoms in each classification (n=10, P=0.01, r > r₀). The higher the total score, the better the binding of PCB molecules with the BphA enzyme, and the lower the biological enrichment. As the number of Cl atoms increases, the total score shows a decreasing trend, which indicating that highly chlorinated PCBs are not easily degraded and more easily enriched in the fish body, that is, the number of Cl atoms increases, their biological enrichment follows enhanced. In addition, a linear correlation between the total score obtained by molecular docking between 209 PCBs and BphA enzyme and the number of corresponding Cl atoms was obtained, the correlation coefficient r² was 0.6757.

It can be seen from Table 6 that when the significance level P=0.01, the r of the four linear fitting relationships is greater than the correlation coefficient test critical value r₀, indicating that the logK_ow values of PCB molecules and the total scores of docking with BphA enzymes (represent biological enrichment) and the number of Cl atoms have a strong correlation.

The logK_ow values used for modeling comes from experimental determinations and belongs to the property values of PCB molecules, but the effect of bioaccumulation on specific receptor proteins is unknown. In order to further investigate the effect of low biological enrichment newly designed PCBs-207 molecules on the bioaccumulation of liver enzymes in fish, 12 newly designed PCB-207 molecules, in which all the POPs characteristic parameters are reduced while the actual functional properties are not changed were respectively docked with aryl hydrocarbon receptor (AHR) [38], superoxide dismutase (SOD) [39], cytochrome P450-dependent monooxygenases (CYP1A) [40], ethoxyresorufin-o-deethylase (EROD) [41], glutathione s-transferase (GST) [42], catalase (CAT) [43], and the aryl hydrocarbon hydroxylase (AHH) [44] to verify whether the 12 newly designed PCB-207 molecules have low bioaccumulation for various enzymes. The structures of seven receptor proteins are all from the Protein Data Bank (http://www.rcsb.org/pdb), SOD is an antioxidant enzyme located in erythrocytes of the liver, which the enzymatic activity reduced due to PCBs entering [45]. EROD enzyme was shown to directly contact PCBs and is regulated by CYP1A enzymes. The higher the concentration of pollutants entering the fish, the more the enzyme quantity, and the higher the enzyme activity [46]. When the pollutants enter the body of the tested zebrafish, the GST in liver activity decreases [47].

As shown in Figure 5, taking the modified molecule 5-amino-PCB-207 as an example, compared with the target molecule PCB-207, the total score of 5-amino-PCB-207 molecular docking with SOD, GST, and AHH had declined. According to the color of the legend, the decreasing order of total score values for 5-amino-PCB-207 docking with enzymes is GST > AHH > SOD; the total score of 5-amino-PCB-207 molecular docking with CYP1A, EROD, AHR, and CAT had increased. According to the color of the legend, the increasing order of total score values for 5-amino-PCB-207 docking with enzymes is CAT > CYP1A > EROD > AHR.

The decrease in total score indicated that the binding of 5-amino-PCB-207 with SOD, GST, and AHH had got worse, which could reduce PCBs’ enrichment in liver. Although the increase in total score indicated that the degree of biological enrichment for 5-amino-PCB-207 in CYP1A, EROD, AHR, and CAT did not decrease, that is, not all of the 12 low biological enrichment new designed PCB-207 molecules reduce the bioaccumulation at liver enzymes, the present paper only provides a method to study the bioaccumulation effects of new molecules on specific receptors by using molecular docking.

Table 7 showed the total scores of 12 low biological enrichment newly designed PCB-207 molecules after docking with seven enzymes. It can be seen from the total scores in Table 7, 5-methyl-PCB-207, 5-amino-PCB-207 and 4-amino-5-ethyl-PCB-207, the binding extent of each newly designed molecular mentioned above to the three enzymes was decreasing compared with that of the target molecular PCB-207. For example, the total scores of 5-methyl-PCB-207 docking with SOD, GST and CAT enzymes were lower than the total scores of the target molecular PCB-207 docking with the above three enzymes; the total scores of 5-amino-PCB-207 docking with SOD, GST, and AHH enzymes were lower than the total scores of the target molecular PCB-207 docking with the above three enzymes; the total scores of 4-amino-5-ethyl-PCB-207 docking with EROD, GST, and AHH enzymes were lower than the total scores of the target molecular PCB-207 docking with the above three enzymes.

1. In the present study, HQSAR model was constructed to determine that the Cl-substitutions Cl₄, Cl₅, and Cl₆ have an uppermost PC to the bioaccumulation activity. Based on 3D contour plot of CoMSIA model, indicated that introducing electron withdrawing or hydrophobic groups into Cl₄ and Cl₅ substitution sites of PCB molecules could decrease the bioaccumulation of target molecular. Coupled with HQSAR and CoMSIA model, reduced the K_ow value of monosubstituted compounds by introducing electron withdrawing or hydrophobic groups at Cl₄, Cl₅ two positions, respectively, and reducing the K_ow value of the bis-substituted compounds by introducing electron withdrawing or hydrophobic groups at Cl₄, Cl₅ positions at same time.

2. Based on HQSAR/CoMSIA model, 12 kinds of novel PCB-207 molecules were designed, which have lower biological enrichment, toxicity, persistence, and long-distance migration ability while the actual functional properties (insulation and flame retardancy) are not changed. Based on the molecular docking, the bioaccumulation mechanism of PCBs shows that the biological enrichment effect of highly chlorinated PCBs is significantly higher than that of PCBs with less Cl substituents. Through the molecular docking, it can also reflect the biological enrichment effect of low biological enrichment new designed PCB-207 molecules on the specific enzymes in the liver of fish and further screen out the environmentally-friendly features of the more obvious low biological enrichment new designed PCB-207 molecules.

Establishment of models, collection, analysis, and interpretation of the data: J.Y. Drafting the article or revising it critically for important intellectual content: Y.L. and W.G. All authors approved the final version for publication and agree to be accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. All persons designated as authors qualify for authorship, and all those who qualify for authorship are listed.

The authors declare that there are no competing interests associated with the manuscript.

This work was supported by the Fundamental Research Funds for the Central Universities in 2013 [grant number JB2013146]: the Study on Toxicity Control of PBDEs Based on Quantum Chemistry and QSAR Model; and the Fundamental Research Funds for the Central Universities in 2017 [grant number 2017XS058]: the Study on Enhancement of Vibrational Spectra by Derivation of Phthalic Acid Esters using Pharmacophore Model.

AHH

aryl hydrocarbon hydroxylase

AHR

aryl hydrocarbon receptor

CAT

catalase

Ch

chirality

CoMFA

comparative molecular field analysis

CoMSIA

comparative molecular similarity indices analysis

CYP1A

cytochrome P450-dependent monooxygenase

EROD

ethoxyresorufin-o-deethylase

GST

glutathione s-transferase

HOMO

highest occupied molecular orbital

HQSAR

hologram quantitative structure–activity relationship

PCB

polychlorinated biphenyl

PLS

partial least square

POP

persistent organic pollutant

SE_cv

cross-validated standard error of estimate

SEE

standard error of estimate

SOD

superoxide dismutase