Investigation on the Binding of Polycyclic Aromatic Hydrocarbons with Soil Organic Matter: A Theoretical Approach

Polycyclic aromatic hydrocarbons (PAHs) are ubiquitous contaminants of the terrestrial environment that have been designated as Environmental Protection Agency (EPA) Priority Pollutants. In this study, molecular modeling was used to examine the physical and chemical characteristics of soil organic matter (SOM), fulvic acid (FA) and humic acid (HA), as well as their binding interactions with PAHs. The molecular structures of 18 PAHs were built by using the SYBYL 7.0 program and then fully optimized by a semiempirical (AM1) method. A molecular docking program, AutoDock 3.05, was used to calculate the binding interactions between the PAHs, and three molecular structure models including FA (Buffle's model), HA (Stevenson's model) and SOM (Schulten and Schnitzer's model). The pi-pi interactions and H-bonding interactions were found to play an important role in the intermolecular bonding of the SOM/PAHs complexes. In addition, significant correlations between two chemical properties, boiling point (bp) and octanol/water partition coefficient (Log K(ow)) and final docking energies were observed. The preliminary docking results provided knowledge of the important binding modes to FA, HA and SOM, and thereby to predict the sorption behavior of PAHs and other pollutants.


Introduction
Polycyclic Aromatic Hydrocarbons (PAHs) are organic environmental contaminants formed during the incomplete combustion of organic fuels such as coal, oil and natural gas.Some are highly toxic, mutagenic and/or carcinogenic to humans [1].The transport, fate, degradation and bioavailability of these contaminants can be strongly influenced by the presence of organic matter with various linkages between the aromatic groups [2][3][4][5][6][7].Due to many physical, chemical and biological processes that occur in soils, the understanding of their properties and behavior is essential and challenging.
Agricultural soils contain approximately 3% soil organic matters (SOM), 3% water and 94% inorganic compounds.The SOM consist of humic substances (humic acid (HA), fulvic acid (FA) and humin) and smaller amounts of carbohydrates, N-containing compounds and lipids [8].Chemically, the structure of humin is similar to that of HA but it is strongly complexed by clays and hydrous oxides and cannot be extracted by either dilute base or acid.HA and FA consist of similar amounts of aliphatic and aromatic carbon atoms.In addition, FA has components with lower molecular weight but it is richer in carboxylic groups.The components in SOM cannot be separated effectively by chemical or physical methods [9].However, from chemical and pyrolysis-mass spectrometric analyses, Schulten et al. [10] found that N-heterocylics are significant components of SOM as pyrolysis of humic substances produces fragments including pyrroles, free and substituted imidazoles, pyrazoles, pyridines, substituted pyrimidines, pyrazines, indoles, quinolines, N-derivatives of benzene, alkylamines, and alkyl and aromatic nitriles.The relation of these fragments to the parent humic structures is not known.Moreover, since the mixture of molecules is large and complex, a structural model is difficult to propose and is still being developed.
Nevertheless, the molecular modeling approach based on some analytical methods is useful for predicting a structural model of SOM now, and several models of FA, HA and SOM have been proposed [11,12].In this study, a plausible soil organic matter (SOM) structure proposed by Schulten and Schnitzer [13] was used.The structure was obtained by performing molecular mechanics and dynamics calculations.This model was constructed based on the results from temperature-programmed pyrolysis in the ion-source of the mass spectrometer combined with soft ionization in very high electric fields (Py-FIMS) and Curie-point pyrolysis-gas chromatography/mass spectrometry (Py-GC/MS).In this model, organic compounds such as carbohydrates and proteinaceous materials were adsorbed in internal voids and on external surfaces.The carbohydrates trapped in this SOM model are represented by a trisaccharide with 66 atoms.For proteinaceous materials, the hexapeptide AspGlyArgGluAlaLys consisting of amino acids typically found in soils was used.This SOM structure has proven that its structural model correlated with important properties of SOM such as mass, surface area, and volume.It can also explain some important characteristics of SOM, for example, surface activity, cation exchange capacity, binding and trapping of biological substances.In addition, two well-known proposed structural models of FA and HA were also selected to study.Fulvic acid (Buffle's model) [14] consists of naphthalene rings substituted with hydroxyl, carboxyl and short aliphatic chains containing alcohol, methyl, carboxyl and carbonyl groups.Humic acid (Stevenson's model) [15] contains phenolic −OH groups, quinone structures, and carboxyl groups substituted on the aromatic rings.Consequently, in this study, three structural models representing the components and the various size of the SOM models were used to study the effect of PAH binding.
The structural models are shown in Figure 1.In order to investigate the binding interactions of PAHs with soil, molecular docking methods were applied to the FA, HA and SOM models.Furthermore, Nègre et al. [16] have pointed out that the extent of adsorption of imidazolinone herbicides on SOM was dependent on the HA characteristics.Therefore, the effect of the carboxylic groups in the ionized form of FA and HA was investigated.
The aim of this study is to investigate binding characteristics of PAHs to three different molecular structure models, and comparing their physicochemical properties with the calculated docking energies.[18].For the ionized structures of FA and HA, their carboxylate groups were ionized and the energies then minimized.They were fully optimized at the HF/3-21G level of theory using the Gaussian03 program [19].The Gasteiger-Hückel charge in the SYBYL 7.0 program was added to the optimized structures.For SOM, the three-dimensional model of Schulten and Schnitzer was used.The model consists of humic acid, carbohydrate (saccharide), hexapeptide (proteinaceous material) and hydrating water molecules, as taken from available structures in "The Virtual Museum of Minerals and Molecules™" [20].The water molecules were removed and Gasteiger-Hückel charges were added to humic acid and saccharide.Kollman charges were added to proteinaceous material using the SYBYL 7.0 program.The molecular structures of the PAHs used in this study are shown in Figure 2. Their starting geometries were built by using the SYBYL 7.0 program and then fully optimized by semiempirical AM1 method by using the Gaussian03 program.The Gasteiger-Hückel charge in the SYBYL 7.0 program was finally added to their optimized structures.

Docking of PAHs to fulvic acid, humic acid and soil organic matter
To study the binding mode of PAHs to FA, HA and SOM, molecular docking was carried out using the AutoDock 3.05 program [21] with the Lamarckian genetic algorithm.The grid size for FA was set to 100 x 100 x 100.For HA and SOM, the grid size was set to 120 x 120 x 120.A grid spacing of 0.375 Å was used.The numbers of docking runs were set to 50 and the default values of all other parameters were used.The PAHs were also docked to the ionized structures of both FA and HA.The same docking parameters were used.

Structural comparison between neutral and ionized fulvic acid
The optimized structures of FA and ionized FA are shown in Figures 3a and 3b.Three intramolecular H-bonds were found in the optimized FA structure.The H-bond between O-11 and the hydroxyl group of O-18 was lost when this carboxylic group is ionized.The repulsive interaction causes the side chain to be displaced in the ionized FA structure.Both optimized structures are compared in Figure 4.These structures were used for further docking steps.

Structural comparison between humic acid and ionized humic acid
The optimized structures of HA and ionized HA are compared in Figure 5. Compared to the optimized HA structure, it was found that the orientation of rings A and B is different in the ionized HA structure, with these rings oriented away from the C and D rings.This change was caused by the stronger H-bond interaction formed in the ionized HA structure.Because of the O-13 in the ionized HA structure, a stronger H-bond interaction was formed between the carboxylate group O-13 and the hydrogen atoms in ring C and the O-9 hydroxyl group.The distances between O-13 and the hydrogen atoms were decreased from 2.51Å to 1.98Å and 4.07Å to 3.53Å, respectively.The other ring orientations for the two structural forms are the same.Both optimized structures were used in further docking steps.

Docking of PAHs
Final docked energies of the 18 PAHs are shown in Table 1.

Docking onto fulvic acid and ionized fulvic acid
The 18 PAHs were docked onto the optimized FA and ionized FA structures.From the docked results, the 50 docked conformations of each PAH were classified into several clusters.In the cluster with the most occurring configurations, the conformation of each PAH which corresponded to the lowest final docked energy was selected, as shown in Figure 6.The selected docked conformations onto FA are shown in Figure 6a.It was found that all docked conformation were oriented parallel to the naphthalene ring plane.Most docked conformations are superposed on the naphthalene ring of FA, except PAH8, PAH16, and PAH18.Because of the orientation of the bulky substituted group attached at ring B of FA, PAH8, PAH16, and PAH18 preferred to orient below the naphthalene ring.The calculated final docked energies of the PAHs onto FA ranged from -3.90 to -6.48 kcal/mol.The highest final docked energy (-3.90 kcal/mol) corresponds to the binding of PAH15, which contains two aromatic rings.The other PAHs with three to six aromatic rings, except PAH2, have the lower final docked energies (-4.21 to -6.48 kcal/mol).In case of ionized FA, as shown in Figure 6b, eleven docked conformations were located below the plane of the naphthalene ring; the exceptions were PAH3, PAH4, PAH9, PAH11, PAH13, PAH15, and PAH17.This finding is possibly due to the fact that because of the displacement of the chain attached to the naphthalene ring, as described in the previous section, the surface areas above and below the naphthalene ring are similar.Compared to FA, the energies of all the PAHs docked onto ionized FA have a small difference of about ±0.5 kcal/mol.The calculated final docked energies of the PAHs onto ionized FA ranged from -3.82 to -6.26 kcal/mol.The highest docked energy (-3.82 kcal/mol) onto ionized FA was again that for the binding of PAH15, which contains two aromatic rings.The increased number of rings caused the tighter binding.Most of the PAHs docked onto ionized FA have higher final docked energies than those docked onto FA, the exceptions being PAH7, PAH8, PAH16 and PAH18.π-π interactions and H-bonding interactions were found in the complex formed between the PAHs and both structures of FA. π-π interactions were formed between the naphthalene rings of both FA structures and the aromatic rings of the PAHs, and H-bonding formed between the hydrogen atoms of the PAHs and oxygen atom of carboxylic groups attached to the naphthalene ring of both FA structures.In addition, oxygen atoms (number 11, 12 and 17) also formed H-bonds to the PAHs.

Docking onto humic acid and ionized humic acid
All the PAHs were docked onto the optimized structures of HA and ionized HA.Their docked results were classified into several clusters.The cluster which has the most occurring configurations was studied.The conformation of each PAH which is the lowest final docked energy in that cluster were selected, as shown in Figure 7.The final docked energies ranged from -4.62 to -8.10 kcal/mol for HA and from -4.11 to -7.14 kcal/mol for ionized HA, respectively.The final docked energies show that all the PAHs bound onto HA more tightly than onto ionized HA.The energy differences ranged from 0.07 to 1.35 kcal/mol.The highest final docked energies of both HA structural forms corresponded to PAH15 which is the smallest compound.The increasing number of rings will decrease the final docked energies.Most docked PAHs onto both HA structural forms showed π-π interactions with rings C and D, and H-pi interactions with rings B and E.Moreover, H-bond interaction were also involved between the PAHs and the two HA structural forms.The docked orientations showed H-bond interaction to the O-6, O-7 and O-17 hydroxyl groups of both forms.Because the orientation of rings A and B twisted to face away from rings C and D in ionized HA, H-pi interactions between the PAHs and ionized HA were decreased.

Docking onto soil organic matter
All the PAHs were docked onto the SOM structure.The docked results for each PAH were classified into clusters.In the cluster with the most occurring configurations, the conformation of each PAH which has the lowest final docked energy was selected.Their selected docked conformations were shown in Figure 8.The final docked energies of the PAHs ranged from -6.04 to -9.70 kcal/mol.Nine of them were located in area-1 by forming π-π interactions with three attached HA phenyl rings and H-bonding with HA and Arg.Six of the PAHs showed π-π interactions and H-bonding to HA in area-2.In area-2, there was a wide space between two phenyl groups, therefore the bigger compounds, ex.PAH8 and PAH14, appear to be confined to this area.The others were in area-3 and formed H-bonds with HA, Asp and Gly.It was found that small compounds, PAH1, PAH2 and PAH15, preferred to bind in area-3 with high final docked energies (-7.07, -6.67 and -6.04 kcal/mol, respectively).

Relationship between some chemical properties and final docked energy of PAHs
Values of two chemical properties used in this study, boiling point (bp) and octanol/water partition coefficient (Log K ow ) describing lipophilicity, are shown in Table 2. Highly significant correlations between these properties and final docked energies were found, as shown in Table 3 and Figure 9. Ref. [22]; b.Ref. [23]; c.Ref. [24].
Therefore, statistically significant associations between boiling point of PAHs and final docked energies of them onto FA, ionized FA, HA, ionized HA, and SOM with R 2 = 0.871, 0.947, 0.916, 0.938, and 0.706 respectively, were shown.The results revealed that PAHs with higher boiling points bind tighter to FA, ionized FA, HA, ionized HA, and SOM models than those with lower boiling points.For the correlation with Log K ow , the plots also showed high correlations with R 2 = 0.872, 0.938, 0.887, 0.931, and 0.720, respectively.PAHs having high Log K ow values appear to be strongly bound to FA, ionized FA, HA, ionized HA, and SOM models.Notably, ionized forms of FA and HA showed higher correlation than those of non-ionized forms in binding PAHs, agreeing with a result from a previous study that the deprotonation of the macromolecules changes their shape, increases their polarity and thereby decreases their ability to sorb hydrophobic molecules such as PAHs [24].From the docked conformations of PAHs onto FA, HA and SOM structures, π-π interactions and H-bonding were important for the binding of PAHs.In addition, H-bonding to proteinaceous material was also found with the SOM model.The final docked energies of the PAHs showed that they bound to SOM more tightly than both forms of HA and FA.The final docked energies of SOM were lower than those of HA or FA by about 0.67 to 2.29 kcal/mol and 1.91 to 3.38 kcal/mol, respectively.With FA, 14 compounds showed lower final docked energies than with ionized FA.The energies of PAHs docked onto ionized FA differed by about ±0.5 kcal/mol.In case of the HA structure, all compounds demonstrated lower final docked energies than with the ionized HA, and ranged from 0.31 to 1.35 kcal/mol.Because of the conformational changes in optimized neutral and ionized forms of FA and HA, the results showed that PAHs tighter bound with the neutral form of FA and HA structures, as shown in the lower final docked energies in their neutral forms.These results agreed well with the experimental results of Nègre et al. [15] that the adsorption of the herbicides onto HA is dependent on the HA characteristics.The extent of adsorption is higher with the neutral form of HA than the anionic form of HA.Based on the high correlation between the physicochemical properties and docked energy, the calculation of final docked energies can be used as an useful tool for predicting the boiling point and Log K ow of PAH homologs.Thus, further research will be needed to find more properties including organic carbon partition coefficient (K oc ), as a sorption index, to predict pollutant transport.
The results indicate that the adsorption of PAHs on to FA, HA, and SOM depend on the properties of PAHs and SOM models (i.e.FA, HA, and SOM).Thus this study revealed the important interactions of PAHs onto FA, HA and SOM.The binding depends on the chemical structures of the organic compounds and also the characteristics of FA and HA.The results are also useful for suggesting the further development of the control of binding between small molecules and soil.

Figure 2 .
Figure 2. Structures of Polycyclic Aromatic Hydrocarbons (PAHs) used in this study.

Figure 8 .
Figure 8. Binding sites (a) and orientation (b) of docked PAHs onto SOM.

Figure 9 .
Figure 9. Plots between two chemical properties (bp (a) and Log K ow (b)) and final docked energies of PAHs .

Table 1 .
Final docked energies of PAHs.

Table 2 .
Boiling point (bp) and Log K ow of PAHs. a.

Table 3 .
Statistical correlations between two chemical properties (bp and Log K ow ) and final docked energy of PAHs.