Computational and Theoretical Chemistry 1043 (2014) 54–63
Contents lists available at ScienceDirect
Computational and Theoretical Chemistry journal homepage: www.elsevier.com/locate/comptc
Interaction between glyphosate and mitochondrial succinate dehydrogenase Ricardo Ugarte ⇑ Instituto de Química, Facultad de Ciencias, Universidad Austral de Chile, Independencia 641, Valdivia, Chile
a r t i c l e
i n f o
Article history: Received 3 April 2014 Received in revised form 28 May 2014 Accepted 29 May 2014 Available online 11 June 2014 Keywords: Glyphosate Succinate dehydrogenase Potential energy surface (PES) ONIOM SCRF DFT calculations
a b s t r a c t A recent study has evaluated the toxicity of the glyphosate-based herbicide, Roundup; the contents of the report evidences that the herbicide inhibits the mitochondrial succinate dehydrogenase activity. The present investigation, tries to prove the interaction of the active compound, glyphosate, with the succinate binding site of the enzyme, by using two-layer ONIOM calculations for structure optimization. In particular, the binding energy of the glyphosate-binding site complex model, in gas phase and polarized continuum medium in order to incorporate the protein environment, is calculated and compared it with a reference malonate-binding site complex model. The binding energy difference between these two model systems, 14 kcal mol1 at B3LYP/6-31+G(d) level of theory, supports the hypothesis of the interaction of glyphosate at the succinate binding site of the mitochondrial succinate dehydrogenase enzyme, verifying thus the experimental fact. Therefore, it is possible to conclude that glyphosate is harmful for the human health. Ó 2014 Elsevier B.V. All rights reserved.
1. Introduction Roundup produced by Monsanto Company is the brand-name of an effective herbicide used worldwide. The chemical on which it is based, glyphosate (N-(phosphonomethyl) glycine), acts through competitive inhibition of the enolpyruvyl-shikimate-phosphate synthase, an enzyme essential to the synthesis of aromatic amino acids in plants . The most of the genetically modiﬁed plants has been designed to tolerate it [2,3] and as result, the expansion of the transgenic crops leads to an increment in the use of the Roundup. Scientiﬁc research discommends their use due to possible adverse effects on the environment and on human health . There is growing evidence that glyphosate (Fig. 1) and its major metabolite aminomethylphosphonic acid (AMPA), apart from the other toxic effects, to cause genotoxicity/carcinogenicity [5–15]. Also, it has been suggested that the presence of Roundup adjuvants enhances glyphosate bioavailability and/or bioaccumulation on cells and amplify the intrinsic toxicity of the compound [3,14]. The toxicity of Roundup formulations on three different human cell types has been evaluated . Cell death was estimated by inhibition of succinate dehydrogenase activity and thus of mitochondrial metabolism , necrosis by release of cytosolic adenylate kinase measuring membrane damage and apoptosis via ⇑ Tel.: +56 63 2293133. E-mail address: [email protected]
http://dx.doi.org/10.1016/j.comptc.2014.05.018 2210-271X/Ó 2014 Elsevier B.V. All rights reserved.
activation of enzymatic caspases 3/7 activity. Besides, it was concluded that the Roundup adjuvants change human cell permeability, contributing to the incorporation of glyphosate into cells; at ﬁrst glance, this is a logical conclusion because at physiological pH the chemical species is mainly the glyphosate dianion , which is not capable of crossing the cell membrane. The mitochondrial succinate dehydrogenase may be a potential target for glyphosate. Since, currently there is not structural information available for glyphosate–succinate dehydrogenase complex, a molecular modeling  study of the interaction between glyphosate and the active site of the mitochondrial enzyme would be able to be useful for demonstrating or to support the feasibility of formation of the complex. Succinate dehydrogenase (SDH) [19,20] is a membrane enzyme essential in intermediary metabolism and aerobic energy production in both prokaryotic and eukaryotic cells. SDH couples the oxidation of succinate to fumarate in the Krebs cycle to the reduction of ubiquinone to ubiquinol in the electron transport chain. The enzyme is a tetramer integrated by two hydrophilic and two hydrophobic subunits. The hydrophilic subunits, composed of a ﬂavoprotein (SDHA) and an iron–sulfur protein (SDHB), have a high degree of sequence homology across species. The SDHA and SDHB catalytic subunits are anchored to the membrane by the hydrophobic SDHC and SDHD subunits (Fig. 2). There are two distinct classes of inhibitors of SDH. The ﬁrst class, which includes succinate analogs like malate, oxaloacetate
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63
Fig. 1. Chemical structures of glyphosate zwitterion (left) and its metabolite AMPA (right).
Fig. 3. Stereoview of the crystallographically determined structure of porcine heart mitochondrial SDH in strand representation. Note that the similarity with the generic structure shown in Fig. 1.
Fig. 2. Model of the tetramer structure of SDH in a phospholipid membrane: SDHA (green), SDHB (cyan), SDHC (magenta) and SDHD (yellow). (For interpretation of the references to color in this ﬁgure legend, the reader is referred to the web version of this article.)
(tricarboxylic acid cycle intermediates) and malonate (synthetic analog), binds to the succinate binding site of the SDHA subunit [21–23]. The second class of inhibitors, which includes the ubiquinone analogs thenoyltriﬂuoroacetone and carboxin, binds to the ubiquinone active site; a hydrophobic pocket comprised of residues from SDHB, SDHC and SDHD subunits  located within the membrane domain of the enzyme. According to the previous information it is assumed that the glyphosate ion binds to the succinate binding site. This choice is based on hydrophilicity criteria and certain structural similarity of glyphosate with the succinate analogues. Certainly, we did not discard, a priori, a future study on the ubiquinone binding site. The crystal structure of human SDH is not yet available; however, we relied on the high degree of sequence homology of the SDHA hydrophilic subunit across species . The model of the succinate binding site is constructed from the X-ray structure of porcine heart mitochondrial complex II deposited in the Protein Data Bank (Fig. 3) . Fortunately, this structure contains a co-crystallized malonate (Fig. 4), which simpliﬁes the search of the succinate binding site and provides, besides, an initial structure for the construction of a reference system model: the malonate-binding site complex. The present investigation is an attempt to model the quantum mechanical behavior of glyphosate ion at the active site of the
Fig. 4. Fragment of porcine heart mitochondrial SDHA subunit. The presence of FAD indicates that malonate is linked to the succinate binding site.
SDHA subunit. The main property to determine is the interaction energy between glyphosate and succinate binding site model of SDHA (BS); this energy is deﬁned as the binding energy of the glyphosate-binding site complex (DEG-BS):
Glyphosate þ BS ! G BS
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63
DEGBS ¼ EGBS ðEGlyphosate þ EBS Þ In order to legitimate the obtained results, DEG-BS must be contrasted with the binding energy of the malonate-binding site complex (DEM-BS):
Malonate þ BS ! M BS
DEMBS ¼ EMBS ðEMalonate þ EBS Þ If the theoretical results predict a reasonable value for DEG-BS, it should be concluded that the formation of the complex is plausible, thus explaining the inhibitory capacity of glyphosate. Accordingly, from a theoretical point of view, it would be proven that glyphosate is potentially harmful. 2. Methods and models The molecular modeling is a powerful tool for the study of the molecular interactions; in particular, between a ligand and its receptor site. Quantum and molecular mechanics are the underlying theories in molecular modeling. These enable to calculate the energy of any arrangement of the atoms in the system and to determine how the energy of the system varies as the positions of the atoms change. By means of high computing power and efﬁcient algorithms, molecular modeling can be effectively used to solve complex chemical and biological problems . In order to model the interaction between glyphosate and SDH, the starting point is provided by the X-ray crystal structure of the target site. In this instance, the X-ray structure of porcine heart mitochondrial complex II bound with N-Biphenyl-3-Yl-2-Triﬂuoromethyl-Benzamide, determined at 3.24 Å resolution and retrieved from the Protein Data Bank (PDB: entry code 3ABV) . The construction of the succinate binding site model was carried out with the ArgusLab molecular modeling software . As from subunit that contains malonate, the computer program generates a model that consists of all residues that have at least one atom within 3.5 Å from any atom of the bound malonate. This generally gives a good representation of the important residues in the binding pocket for a protein target. The residues are: Alanine (A) 412, Arginine (R) 298 and 409, Glutamine (Q) 62, Glutamic acid (E) 267, Glycine (G) 63, Histidine (H) 254 and 365, Leucine (L) 264, Phenylalanine (F) 131, Threonine (T) 266 and Flavin Adenine Dinucleotide (FAD) moiety. The histidine amino acids with hydrogen on the epsilon nitrogen, Q62-G63 and T266-E267 are dipeptides. There are ten residues altogether; obviously all these residues should belong to the sucinate binding site of the SDHA subunit. The weakness of protein crystallography is the inability of the method to locate the positions of hydrogen atoms, and therefore, these are not included in the PDB coordinate ﬁle. In addition, the amino acids residues generated by the ArgusLab software lack of –OH groups (in general, the peptide bonds are broken). These issues are resolved by addition of the correspondent atoms and subsequently, deleting the hydrogen atoms in excess. In this stage, it is very important to check the molecular structures and proceed to rectify incongruences. The GaussView graphical software  was used for editing tasks of the models. The pre-optimized (vide infra) original active site (305 atoms) generated by computer is shown in Fig. 5; by default GaussView yields the amino acids in its zwitterionic conﬁgurations. The FAD moiety was ‘‘pruned’’ to reduce the computational cost (Fig. 6c). This ‘‘FAD fragment’’ (isoalloxazine tricyclic ring derivative) along with the un-ionized conﬁgurations of the amino acid residues constitute the malonate-binding site complex (M-BS), that contains 250 atoms; additionally, by virtue of the malonate elimination the binding site model (BS) is obtained.
Docking calculations, implemented by the ArgusLab software , attempt to place ligands into binding sites; in this case, glyphosate (Fig. 6b) into succinate binding site. The procedure is identical to the above, except that malonate is previously eliminated of the succinate binding site of SDHA and replaced for glyphosate ion. Thus, the glyphosate-binding site complex (G-BS) is obtained. All calculations were performed using the Gaussian 09 program package . Initially, the models were submitted to a PM6 semiempirical  geometry pre-optimization, in which the original coordinates of the PDB ﬁle are frozen; i.e., only the coordinates of the added atoms are optimized. The coordinates of these debugged models were taken as starting point for the rest of the calculations. In order to decrease the computational effort, the ONIOM method  was used to calculate the optimized geometries of the models. This hybrid method offers a solution to the calculations on very large systems because each region is treated with a different computational approach. It often turns out that expensive computational methods are only required for the ‘‘principal’’ part of the system, while the less expensive ones can be used for the ‘‘secondary’’ regions. In the present study, we carry out two-layer ONIOM calculations treating the ligand (malonate and glyphosate) at the RHF or DFT level of theory [31,32] and the binding site with PM6 semi-empirical method. In general, the following constraints were applied to the three models (BS, M-BS and G-BS): - The a-carbon position of all amino acids was frozen, in order to conserve the ‘‘framework’’ of the binding site and therefore, to prevent that it collapses during geometry optimization. However, this general condition imposed is not enough and adequate restrictions have to be introduced to each particular case. - Six dihedral angles were restricted for the purpose of maintaining the planarity of the peptide bond of both dipeptides (1 + 1) and the planar conformation of FAD fragment (4). Since the model systems are immersed in a protein environment, it is more realistic to consider the effect of the medium on the structures. This effect can be reduced to steric restraints on the binding site and polarization of the molecules; the former already has somehow been considered although for other reasons (vide supra). To account for the polarization it is assumes that the protein environment is represented by a homogenous polarizable medium with some value of dielectric constant . Gaussian provides a reasonable alternative to incorporate this effect, because allows that the system be studied via the self-consistent reaction ﬁeld (SCRF) approach, in which the medium is represented by a polarizable continuum model (PCM)  and the system is placed inside a cavity within the medium. The same method must be included in any realistic calculation involving charged species; without it, for instance, acidic groups can never become ionized and the charge distribution on the substrate will not be reasonable. 2.1. Malonate-binding site complex and binding site model The optimization of M-BS is based on the atomic coordinates obtained from the PDB ﬁle. It is supposed that the malonate conformation is the appropriate. Geometry optimization in vacuo of M-BS was performed at ONIOM(RHF/6-31+G(d):PM6) level. The obtained structure was submitted to single point calculation at RHF/ 6-31+G(d) and B3LYP/6-31+G(d) level . The optimized M-BS was re-optimized at ONIOM(B3LYP/6-31+G(d):PM6) and the respective single point calculation was carried out at B3LYP/ 6-31+G(d). The same procedure was utilized for BS, except that the optimization step was performed at PM6 level.
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63
Fig. 5. Stereoview of the binding site (stick representation). Malonate is shown in spaceﬁll model. The ﬁgure displays all residues that have at least one atom within 3.5 Å from any atom of the malonate.
Fig. 6. Chemical structures: (a) Malonate dianion. (b) Glyphosate dianion. (c) FAD fragment.
2.2. Glyphosate-binding site complex The optimization procedure is more complicated than in the previous case, because one must ﬁrst to construct the in vacuo ONIOM potential energy surface (PES)  at RHF/ 6-31+G(d):PM6 level, of the glyphosate dianion (RHF/6-31+G(d)) at the inside of the binding site model (PM6: all frozen atomic coordinates). The knowledge of the PES allows the determination of the set of minimal energy conformations. The PES, as described by the equation E = E(U1, U2), was computed by varying the consecutive dihedral angles U1:P6–C4–N9–C1 and U2:C4–N9–C1–C2 (Fig. 6b);
these angles determine the possible conformations of the molecule. Bear in mind that the dihedral angle U1 of the bonded atoms P6–C4–N9–C1 is the angle between the planes P6–C4–N9 and C4–N9– C1. The PES was calculated using a 24 24 grid generated by rotating through U1 and U2 in 15° increments from 180° to +180°. At each point glyphosate geometry optimization was carried out with U1 and U2 frozen at their respective grid values. Conformational energy map (contour diagram) was drawn in order to facilitate the search of stationary points. Glyphosate unconstrained geometry optimizations was performed to a set of conformations that lie close to the minima found on the PES. Thus, the most stable
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63
conformers of glyphosate in the molecular ﬁeld of the ﬁxed BS (all frozen atomic coordinates) were obtained. Finally, these structures were submitted to the same method utilized previously for M-BS. Once the in vacuo stable models were obtained, the medium effect on the structure and energy of these was studied. Geometries optimization and single point calculations were carried out in the presence of a structureless polarizable continuum medium characterized mainly by its dielectric constant e; different values of e were considered in this work. It is pertinent to mention that the PCM geometries optimization was a very tedious and time-consuming process; concretely, the Gaussian program was aborted frequently. 2.3. Dispersion corrected DFT functionals A recent published work , supports the notion that dispersion forces are import in non-covalent interactions between isoalloxazine ring of FAD and three prototypical aromatic ligands. In this regards, because B3LYP hybrid functional is lacking dispersion terms, a preliminary study for comparative purposes was performed. The binding energy of the complexes was calculated using M06 and wB97xD dispersion corrected DFT functionals [38–41]. M06 and wB97xD single point calculations were carried out using the following optimized geometries at the respective levels of theory: the complexes at ONIOM(B3LYP/6-31+G(d):PM6), BS at PM6 and ligands at B3LYP/6-31+G(d). 3. Results and discussion The conformational energy map (contour diagram) of glyphosate at the inside of the ﬁxed BS is shown in Fig. 7. The global minimum is located in [U1, U2] = [+90°, 60°]. A comparison with the PES of the isolated molecule (Fig. 8) provides a good semi-quantitative picture of the BS inﬂuence, which is reﬂected as: an increase in the relative energy of the conformations, tendency towards a reduction in the number of the mirror images of glyphosate, increase of the energy at the ‘‘peripheral region’’ of the MEP and as a consequence, that the more stable structures can be found in the ‘‘central region’’ and furthermore, due to the translation of the glyphosate as a whole, the PES is more irregular and less smooth in appearance. The previous results can
Fig. 7. Conformational energy map in vacuo of glyphosate at the inside of the ﬁxed binding site model. DERHF (kcal mol1) relative to the global minimum of the ONIOM-PES: [U1, U2] = [+90°, 60°].
Fig. 8. Conformational energy map in vacuo of glyphosate. DERHF (kcal mol1) relative to the global minimum of the PES: [U1, U2] = [+160°, 160°] ([160°, +160°]).
be explained simultaneously by anisotropicity of the molecular ﬁeld (BS) and, of course, the electrostatic and steric effects that accompany the interaction. The information obtained from the ONIOM-PES allowed selecting 51 G-BS conformations out of a total of 576 calculated structures. Finally, only ﬁve lowest energy conformers were chosen. It should be pointed out that the optimization procedure was not free of problems and additional restrictions had to be imposed to BS. The number of these depend on each particular case and normally, involves constraints to bonds that contain hydrogen atoms; if are not taken into account, unexpectedly the calculation crashes. As can be seen in Table 1, all complexes are stable with large negative binding energies at both RHF/6-31+G(d)//ONIOM(RHF/631+G(d):PM6) and B3LYP/6-31+G(d)//ONIOM(RHF/6-31+G(d): PM6) levels of theory. The above result is not meaningful because it is not realistic to consider isolated or ‘‘naked’’ ions in the vacuum; the energy released in the course of the formation of the complex comes mainly from the charge redistribution over the whole system. Nevertheless, in this initial stage of the study an interesting detail is revealed: the ‘‘shrinkage’’ of glyphosate in the active site (see U1 or U2 values) and the presence of an intramolecular hydrogen bond, e.g., G-BS_II (see interatomic distances), absent in the isolated molecule (G). Apparently, the steric factor, such as the strain of glyphosate (‘‘shrinkage’’) on binding site is reduced by the presence of the intramolecular hydrogen bond; i.e., in vacuo the strained conformation of glyphosate inside the binding site is stabilized to some extent by H-bond. Table 2 shows that environment had a direct consequence on binding energies of the complexes; in relation to the structural properties do not seem to contribute signiﬁcantly. In this approach the binding energy is evaluated considering the energy of the ligand in the respective medium characterized by its dielectric constant. A notorious decrease of the binding energy is observed, particularly for M-BS, and in some cases this change of sign. The stabilization energy of the isolated ligand (L), DES_L = E(L)mediumE(L)in vacuo, is higher than the stabilization energy of the ligand inside the complex (L-BS), DES_L[BS] [E(L-BS)medium–E(L-BS)in vacuo] [E(BS)medium E(BS)in vacuo], because BS already have exerted a stabilizing effect in the vacuum by the charge redistribution. For example, at RHF/6-31+G(d) level and e = 4, DES_L of malonate is 155.0 kcal mol1 and DES_L[BS] of malonate inside
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63 Table 1 In vacuo energies (Hartree), dihedral angles (degree), interatomic distances (Å) and binding energies of the complexes DE (kcal mol1) obtained by several methods. Structure
RHF//ONIOM(RHF:PM6) G-BS_I G-BS_II G-BS_III G-BS_IV G-BS_V M-BS BS G M
7006.4645 7006.5023 7006.4867 7006.4773 7006.4602 6534.0092 RHF//PM6 6119.6507 RHF//RHF 886.6689 414.1903
92.5 126.4 68.5 88.4 55.0
73.4 64.5 100.0 176.0 97.6
4.66 3.72 3.85 4.19 1.82
5.26 1.78 1.94 4.62 3.57
90.93 114.65 104.86 98.96 88.23 105.55
7046.8109 7046.8377 7046.8256 7046.8220 7046.8026 6572.9852 B3LYP//PM6 6156.3607 B3LYP//RHF 890.2906 416.4407
100.15 116.97 109.37 107.12 94.94 115.34
a The geometries were optimized at the ONIOM(RHF/6-31+G(d):PM6) level and then a single point calculation on these optimized structures was performed at the RHF/631+G(d) level of theory. b As mentioned above, except that the single point calculation at the B3LYP/6-31+G(d) level.
Table 2 Effect of the medium, characterized by dielectric constant (e), on the same properties mentioned in Table 1. The self-consistent reaction ﬁeld (SCRF) calculations were carried out with the polarizable continuum model (PCM). The more stable structures of glyphosate (G) are H-bonded. Structure
G-BS_I G-BS_II G-BS_III G-BS_IV G-BS_V M-BS
7006.6526 7006.6773 7006.6675 7006.6443 7006.6447 6534.1919 RHF//PM6 6119.7429 RHF//RHF 886.8863 414.4373
BS G M
e = 12 G-BS_I G-BS_II G-BS_III G-BS_IV G-BS_V M-BS BS G M
e = 20 G-BS_I G-BS_II G-BS_III G-BS_IV G-BS_V M-BS BS G M
88.8 127.0 68.4 89.8 57.0
71.7 64.4 103.4 178.6 100.2
4.56 3.73 3.80 4.14 1.83
5.17 1.79 1.88 4.54 3.60
7006.7096 7006.7326 7006.7316 7006.7066 7006.7019 6534.2552 RHF//PM6 6119.7778 RHF//RHF 886.9374 414.4937
89.7 128.4 69.2 87.2 59.1
71.9 62.5 104.0 177.3 96.9
4.62 3.75 3.79 4.23 1.83
5.20 1.80 1.85 4.58 3.63
7006.7224 7006.7686 7006.7426 7006.7185 7006.7207 6534.2720 RHF//PM6 6119.7871 RHF//RHF 886.9478 414.5052
89.4 128.0 70.6 -84.3 59.7
72.8 61.0 105.3 173.9 96.8
4.60 3.82 3.80 4.29 1.82
5.17 1.82 1.84 4.63 3.62
BS is 56.8 kcal mol1; i.e., the stabilization of malonate ion is the main effect of medium. Thus, the formation of the complex in the environment evolves less energy that in the vacuum. The results also show that inclusion of electron correlation, through B3LYP hybrid functional in DFT, is important in the calculation of binding energies; these become more negative. The environment affects the structure of the isolated glyphosate, enabling the formation of an intramolecular H-bond (see geometric parameters of G in Tables 1 and 2). Although the difference in the binding energy
14.68 30.18 24.03 9.48 9.73 7.34
7046.9882 7047.0019 7046.9960 7046.9860 7046.9770 6573.1574 B3LYP//PM6 6156.4435 B3LYP//RHF 890.5044 416.6827
25.29 33.89 30.18 23.91 18.26 19.58
3.51 10.92 10.29 5.40 8.35 10.23
7047.0416 7047.0537 7047.0560 7047.0424 7047.0298 6573.2160 B3LYP//PM6 6156.4757 B3LYP//RHF 890.5545 416.7380
7.15 14.75 16.19 7.66 0.25 1.44
7.84 21.15 4.83 10.29 8.91 12.74
7047.0529 7047.0857 7047.0664 7047.0545 7047.0466 6573.2326 B3LYP//PM6 6156.4840 B3LYP//RHF 890.5646 416.7493
2.70 23.28 11.17 3.70 1.26 0.44
between the in vacuo and H-bond structures of glyphosate is not large (1 kcal mol1 at RHF//RHF and 2.5 kcal mol1 at B3LYP// RHF on the average), this tends to increase with the dielectric constant. Trends in the electronic charge distribution may provide some additional insight into the intramolecular hydrogen bond in glyphosate. For G-BS_II, Table 3 displays the computed Mulliken charges  of the atoms involved in the hydrogen bond interaction, i.e. H(16), O(5) and O(15). A remarkable decrease in the absolute values of the charges is observed in going from RHF/6-31+G(d)
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63
Table 3 Mulliken charges of G-BS_II calculated at different levels of theory. The complex was optimized at ONIOM(RHF/6-31+G(d):PM6) level. Atom
H(16) O(5) O(15)
0.739 –0.572 –0.788
0.733 –0.564 –0.760
0.744 –0.566 –0.787
0.734 –0.630 –0.860
to B3LYP/6-31+G(d); in consequence, the electrostatic interaction between both oxygens is less repulsive with B3LYP/6-31+G(d). Additionally, the O5–H16 bond distance of 1.80 Å, O5–O15 interatomic distance of 2.75 Å and O5–H16–O15 angle of 170° indicates an intramolecular hydrogen bond of medium strength. Another aspect to have under consideration relates with the concept of complementarity. According to this, the binding energy of the ligand and the site will be at a maximum when the binding site has a structure complementary to that of the ligand . As can be seen in Table 2, if we focused on the geometric parameters of G-BS_II, there is a good correspondence between the isolated glyphosate and the glyphosate at the inside of the BS. That is, glyphosate does not need modifying too much its structure at the inside of the site. As a consequence of binding energies, G-BS_II was chosen to continue the study. All structures were re-optimized at the corresponding level of theory: the complexes (see Fig. 9) at ONIOM(B3LYP/6-31+G(d):PM6) and ligands at B3LYP/6-31+G(d). BS already has been previously optimized at PM6 level. The obtained structures were submitted to single point calculations at B3LYP/ 6-31+G(d) level. In Table 4 are shown the calculated energies and binding energies as a function of the dielectric constant. The dihedral angles (U1, U2) and O5–H16 bond distance, practically did not change with the dielectric constant; the obtained values for glyphosate inside the binding site were: U1 127°, U2 59°, O5–H16 = 1.70 Å. The effect of the dielectric constant on the binding energies of complexes is shown in Fig. 10. The average binding
0.655 –0.285 –0.510
B3LYP/6-31+G(d) e 4
0.650 –0.286 –0.488
0.661 –0.289 –0.507
0.663 –0.380 –0.603
energy difference between M-BS and G-BS_II is 15 kcal mol1 favoring the formation of the latter complex. The tendency of the binding energy to remain constant for e > 10 is because the effect on the change in the energy of the ligand + BS is practically the same one than in the complex. Besides, if the value of the dielectric constant is large the energy of the species is less sensitive to the variation of this parameter. Traditionally, a low dielectric constant of 4 is used, to model the protein environment [33,44]. Nevertheless, the literature reports on studies in which the dielectric constant ranged from 1 to 20 [45,46]. Owing to the hydrophilic character of SDHA, the study focuses on the results obtained for e = 8. In order to verify the distortion that can experience the binding site due to the presence of the ligand, when the complex is formed, the mean value of the deviations in the atomic positions of BS, when both complexes are formed, was calculated: BS in G-BS 0.67 Å and BS in M-BS 0.55 Å. Also, the mean value of the deviations in the atomic positions of ligands, when both complexes are formed, was calculated: G in G-BS 0.11 Å and M in M-BS 0.044 Å. As it can be appreciated, the obtained values are small, corroborating the principle of complementarity. Although the respective differences in the mean values are not signiﬁcant, 0.12 Å and 0.066 Å, malonate distorts less the BS and turns out to be less modiﬁed. Assuming that the distortion at the binding site is of similar magnitude, the difference in the binding energy of the complexes (see Fig. 10) can be attributed to several factors: number and
Fig. 9. Stereoview of G-BS_II (top) and M-BS (bottom) optimized at ONIOM(B3LYP/6-31+G(d):PM6) level (e = 8). Glyphosate and malonate are shown in ball & stick models. The ﬁgures display the residues that constitute the binding site (stick representation). Clearly is observed the glyphosate intramolecular hydrogen bond. The hydrogen bonds are indicated by dashed lines.
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63 Table 4 Calculated energies (Hartree) and binding energies (kcal mol1) at B3LYP/6-31+G(d) as a function of the dielectric constant.
4 8 12 16 20
7047.0044 7047.0421 7047.0604 7047.0719 7047.0788
6573.1593 6573.1997 6573.2167 6573.2250 6573.2329
6156.4435 6156.4673 6156.4757 6156.4805 6156.4840
890.5108 890.5483 890.5609 890.5673 890.5711
416.6865 416.7280 416.7419 416.7489 416.7531
31.44 16.63 14.94 15.12 14.87
18.38 2.80 0.56 2.75 2.69
Table 6 Calculated binding energies (kcal mol1) at M06/6-31+G(d) and wB97xD/6-31+G(d) level of theory (e = 8). For details, see Section 2.3 in Methods and models.
Fig. 10. Binding energies of complexes vs dielectric constant. The data points of Table 4 were ﬁtted by cubic spline interpolation.
Table 5 Intermolecular H-bonds. M-BS and G-BS_II optimized at ONIOM(B3LYP/631+G(d):PM6) level (e = 8). Bond distances between ligand and residues from up to 2 Å have been chosen. Residue
Malonate oxygen atoma BL
E267 G63 H254 H365 Q62 R298 R409 a b
Glyphosate oxygen atomb FR
From Fig. 9 (M-BS): BL (back left), BR (back right), FL (front left), FR (front right). The numbering of atoms as in Fig. 6b.
strength of the intermolecular H-bonds (see Table 5), size and shape of ligands, a possible interaction ligand-FAD fragment, etc. In contrast to ionic or van der Waals interactions, hydrogen bonds are directional and thus play a major role in conferring ligandbinding speciﬁcity. If the ligand is substituted by one of the different shape and size, these interactions will be compromised because the H-bond lengths and angles will change . The interaction between glyphosate and FAD fragment could be supported for the proximity of the ligand to the plane of the tricyclic ring: H(10), H(13) and N(9) atoms of glyphosate at a distance, respectively, of 2.5 Å, 2.5 Å and 3.5 Å from the plane of isoalloxazine ring (N(9) on top of the center of the rings). By contrast, this plane lies 3.1 Å from the closest atom of malonate (hydrogen). Therefore, FAD interacts more weakly with malonate. Generally speaking, the afﬁnity of a ligand toward its binding site falls into a range of approximately 2–20 kcal mol1, and it is
D [DEc] = DEM-BS DEG-BS_II.
related to the binding constant and corresponds to the Gibbs free energy difference for the process. Consequently, enthalpic and entropic effects determine the binding afﬁnity . Since the frequency calculations of BS and the complexes are extremely time-consuming, the thermal corrections and therefore the entropic effects, has not been taken into account in this work. By assuming that the corresponding entropy changes for both complexes are of about the same magnitude, then the binding energy difference could be maintained and being approximately equivalent to the Gibbs free energy difference. In other respects, if the phenomenon of enthalpy–entropy compensation is fulﬁlled [49,50], the binding energy difference could be less pronounced. Another drawback of the method is that it does not take into consideration the solvation and desolvation effects of the ligand and the binding site. In consequence, because of the nature of the study it is fundamental to have a reference system in order to compare the binding energy of the G-BS complex. The reference system is M-BS complex and it is a conﬁrmed fact that malonate interacts with the succinate binding site of the SDHA subunit. The data from Table 4 shows that the difference in binding energy (e = 8) between M-BS and G-BS is 13.83 kcal mol1. Thus, calculations suggest the possibility that glyphosate and SDHA effectively interact yielding a complex; that is, glyphosate could play a critical role in inhibiting the enzyme. Finally, in Table 6 are shown the binding energies calculated using the dispersion corrected DFT functionals. Clearly, there is a notorious tendency to overestimate the strength of these interactions; probably, the single point calculations should be performed on ONIOM(M06/6-31+G(d):PM6) and ONIOM(wB97xD/6-31+G(d): PM6) optimized complexes, and so on. For more details, a critical analysis of dispersion corrected functionals can be found in elsewhere . In any case, the binding energy difference (D[DEc]) is favorable to the formation of the G-BS complex. 4. Conclusions A general method to obtain the binding energy of a ligand to a protein binding site is formulated on the basis of ONIOM and standard calculations. Although the method is approximate, it is applicable to large systems and represents the beginning for the development of a more exact method. The potential energy surface (ONIOM-PES) of glyphosate at the inside of the ﬁxed binding site was calculated. Five minimum energy conformers of glyphosate in the molecular ﬁeld of the binding site were obtained and G-BS_II complex was chosen for the rest
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63
of calculations. I think that it could be the ﬁrst time that a PES calculation of this kind is carried out. It is important to include a dielectric medium, which has a very important effect on the energy and charge distribution of the system. For instance, the in vacuo difference in binding energy at B3LYP//RHF level of theory between M-BS and G-BS_II (see Table 1) is just of 1.6 kcal mol1. The main conclusion from the present study is that the binding energy difference between these two model systems, validates the hypothesis of the possible interaction of glyphosate to the binding site of mitochondrial succinate dehydrogenase, indicating glyphosate could play a critical role in inhibiting the enzyme. Theoretical calculations show that this difference at B3LYP/6-31+G(d) level of theory is approximately of 14 kcal mol1 for e = 8. The results conﬁrm the experimental fact and highlight the risk of glyphosate use for human health. Acknowledgments The author wish to acknowledge to the Dirección de Investigación y Desarrollo de la Universidad Austral de Chile (Grant S-201204) for ﬁnancial support. References  N. Amrhein, B. Deus, P. Gehrke, H.C. Steinrucken, The site of the inhibition of the shikimate pathway by glyphosate. II. Interference of glyphosate with chorismate formation in vivo and in vitro, Plant Physiol. 66 (1980) 830–834.  C. Bolognesi, S. Bonatti, P. Degan, E. Gallerani, M. Peluso, R. Rabboni, P. Roggieri, A. Abbondandolo, Genotoxic activity of glyphosate and its technical formulation roundup, J. Agric. Food Chem. 45 (1997) 1957–1962.  S. Richard, S. Moslemi, H. Sipahutar, N. Benachour, G.E. Seralini, Differential effects of glyphosate and roundup on human placental cells and aromatase, Environ. Health Perspect. 113 (2005) 716–720.  M. Antoniou, M. Ezz El-Din Mostafa Habib, C.V. Howard, R.C. Jennings, C. Leifert, R.O. Nodari, C. Robinson, J. Fagan, Roundup and Birth Defects: Is the public being kept in the dark? June 2011.
.  A. Paganelli, V. Gnazzo, H. Acosta, S.L. Lopez, A.E. Carrasco, Glyphosate-based herbicides produce teratogenic effects on vertebrates by impairing retinoic acid signaling, Chem. Res. Toxicol. 23 (2010) 1586–1595.  C. Elie-Caille, C. Heu, C. Guyon, L. Nicod, Morphological damages of a glyphosate-treated human keratinocyte cell line revealed by a micro- to nanoscale microscopic investigation, Cell. Biol. Toxicol. 26 (2010) 331–339.  J. George, S. Prasad, Z. Mahmood, Y. Shukla, Studies on glyphosate-induced carcinogenicity in mouse skin: a proteomic approach, J. Proteomics. 73 (2010) 951–964.  F. Mañas, L. Peralta, J. Raviolo, H. Garcıa Ovando, A. Weyers, L. Ugnia, M. Gonzalez Cid, I. Larripa, N. Gorla, Genotoxicity of AMPA, Genotoxicity of AMPA, the environmental metabolite of glyphosate, assessed by the Comet assay and cytogenetic tests, Ecotoxicol. Environ. Saf. 72 (2009) 834–837.  F. Mañas, L. Peralta, J. Raviolo, H. Garcıa Ovando, A. Weyers, L. Ugnia, M. Gonzalez Cid, I. Larripa, N. Gorla, Genotoxicity of glyphosate assessed by the comet assay and cytogenetic tests, Environm. Toxicol. Pharmacol. 28 (2009) 37–41.  M. Hultberg, Cysteine turnover in human cell lines is inﬂuenced by glyphosate, Environ. Toxicol. Pharmacol. 24 (2007) 19–22.  A. Anadón, M.R. Martínez-Larranaga, M.A. Martínez, V.J. Castellano, M. Martínez, M.T. Martin, M.J. Nozal, J.L. Bernal, Toxicokinetics of glyphosate and its metabolite aminomethyl phosphonic acid in rats, Toxicol. Lett. 190 (2009) 91–95.  M. Peluso, A. Munnia, C. Bolognesi, S. Parodi, 32P-postlabeling detection of DNA adducts in mice treated with the herbicide roundup, Environ. Mol. Mutagenesis 31 (1988) 55–59.  G.L. Poletta, A. Larriera, E. Kleinsorge, M.D. Mudry, Genotoxicity of the herbicide formulation RoundupÒ (glyphosate) in broad-snouted caiman (Caiman latirostris) evidenced by the Comet assay and the Micronucleus test, Mutation Res. 672 (2009) 95–102.  N. Benachour, G.-E. Seralini, Glyphosate formulations induce apoptosis and necrosis in human umbilical, embryonic, and placental cells, Chem. Res. Toxicol. 22 (2009) 97–105.  C. Gasnier, N. Benachour, E. Clair, C. Travert, F. Langlois, C. Laurant, C. DecroixLaporte, G.-E. Séralini, Dig1 protects against cell death provoked by glyphosate-based herbicides in human liver cell lines, J. Occupational Med. Toxicol. 5 (2010) 29.  F. Peixoto, Comparative effects of the roundup and glyphosate on mitochondrial oxidative phosphorylation, Chemosphere 61 (2005) 1115– 1122.
 D.J. Wauchope, Acid dissociation constants of arsenic acid, methylarsonic acid (MAA), Dimethylarsinic acid (Cacodylic Acid), and N-(Phosphonomethy1)glycine (Glyphosate), Agric. Food Chem. 24 (1978) 717–721.  A.R. Leach, Molecular Modelling: Principles and Applications, second ed., Prentice Hall, 2001.  R. Horseﬁeld, V. Yankovskaya, G. Sexton, W. Whittingham, K. Shiomi, S. Omura, B. Byrne, G. Cecchini, S. Iwata, Structural and computational analysis of the quinone-binding site of complex II (succinate–ubiquinone oxidoreductase), J. Biol. Chem. 281 (2006) 7309–7316.  K.S. Oyedotun, B.D. Lemire, The quaternary structure of the saccharomyces cerevisiae succinate dehydrogenase, J. Biol. Chem. 279 (2004) 9424–9431.  W.C. Kenney, The reaction of N-ethylmaleimide at the active site of succinate dehydrogenase, J. Biol. Chem. 250 (1975) 3089–3094.  F.L. Muller, Y. Liut, M.A. Abdul-Ghani, M.S. Lustgarten, A. Bhattacharya, Y.C. Jang, H. Van Remmen, High rates of superoxide production in skeletal-muscle mitochondria respiring on both complex I- and complex II-linked substrates, Biochem. J. 409 (2008) 491–499.  O.S. Hajjawi, Succinate Dehydrogenase: Assembly, Regulation and Role in Human Disease, Eur. J. Sci. Res. 51 (2011) 133–142.  S. Harada, T. Sasaki, M. Shindo, Y. Kido, D.K. Inaoka, J. Omori, A. Osanai, K. Sakamoto, J. Mao, S. Matsuoka, M. Inoue, T. Honma, A. Tanaka, K. Kita, Deposition date: 2009-12-22, Release date: 2011-02-09. To be published.
.  M.A. Thompson, ArgusLab 4.0, Planaria Software, LCC, Seattle, WA.
.  R. Dennington, T. Keith, J. Millam, K. Eppinnett, W.L. Hovell, R. Gilliland, GaussView, Version 3.09, Semichem Inc, Shawnee Mission, KS, 2003.  S. Joy, P.S. Nair, R. Hariharan, M.R. Pillai, Detailed comparison of the proteinligand docking efﬁciencies of GOLD, a commercial package and ArgusLab, a licensable freeware, Silico Biol. 6 (2006) 601–605.  M.J. Frisch, G.W. Trucks, H.B. Schlegel, G.E. Scuseria, M.A. Robb, J.R. Cheeseman, G. Scalmani, V. Barone, B. Mennucci, G.A. Petersson, H. Nakatsuji, M. Caricato, X. Li, H.P. Hratchian, A.F. Izmaylov, J. Bloino, G. Zheng, J.L. Sonnenberg, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, J.A. Montgomery Jr., J.E. Peralta, F. Ogliaro, M. Bearpark, J.J. Heyd, E. Brothers, K.N. Kudin, V.N. Staroverov, R. Kobayashi, J. Normand, K. Raghavachari, A. Rendell, J.C. Burant, S.S. Iyengar, J. Tomasi, M. Cossi, N. Rega, J.M. Millam, M. Klene, J.E. Knox, J.B. Cross, V. Bakken, C. Adamo, J. Jaramillo, R. Gomperts, R.E. Stratmann, O. Yazyev, A.J. Austin, R. Cammi, C. Pomelli, J.W. Ochterski, R.L. Martin, K. Morokuma, V.G. Zakrzewski, G.A. Voth, P. Salvador, J.J. Dannenberg, S. Dapprich, A.D. Daniels, Ö. Farkas, J.B. Foresman, J.V. Ortiz, J. Cioslowski, D.J. Fox, Gaussian 09, Revision A.01, Gaussian Inc., Wallingford, CT, 2009.  J.J.P. Stewart, Optimization of parameters for semiempirical methods V: Modiﬁcation of NDDO approximations and application to 70 elements, J. Mol. Model. 13 (2007) 1173–1213.  T. Vreven, K. Morokuma, Hybrid methods: ONIOM(QM:MM) and QM/MM, Ann. Rep. Comput. Chem. 2 (2006) 35–51.  D.R. Hartree, The Calculation of Atomic Structures, John Wiley & Sons Inc., New York, 1957.  R.G. Parr, W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, Oxford, 1989.  F. Himo, Quantum chemical modeling of enzyme active sites and reaction mechanisms, Theor. Chem. Acc. 116 (2006) 232–240.  S. Miertus, E. Scrocco, J. Tomasi, Electrostatic interaction of a solute with a continuum. A direct utilization of ab initio molecular potentials for the prevision of solvent effects, Chem. Phys. 55 (1981) 117–129.  P.J. Stephens, F.J. Devlin, C.F. Chabalowski, M.J. Frisch, Ab Initio calculation of vibrational absorption and circular dichroism spectra using density functional force ﬁelds, J. Phys. Chem. 98 (1994) 11623–11627.  R. Ugarte, C. Bustos, I. Moreno-Villoslada, Theoretical study of the isomerization of maleic acid into fumaric acid, J. Chil. Chem. Soc. 56 (2011) 656–662.  L. Koziol, N. Kumar, S.E. Wong, F.C. Lightstone, Molecular recognition of aromatic rings by ﬂavin: electrostatics and dispersion determine ring positioning above isoalloxazine, J. Phys. Chem. A 117 (2013) 12946– 12952.  Y. Zhao, D.G. Truhlar, The M06 suite of density functionals for main group thermochemistry, thermochemical kinetics, noncovalent interactions, excited states, and transition elements: two new functionals and systematic testing of four M06-class functionals and 12 other functionals, Theor. Chem. Acc. 120 (2008) 215–241.  M. Walker, A.J.A. Harvey, A. Sen, C.E.H. Dessent, Performance of M06, M06-2X, and M06-HF density functionals for conformationally ﬂexible anionic clusters: M06 functionals perform better than B3LYP for a model system with dispersion and ionic hydrogen-bonding interactions, J. Phys. Chem. A 117 (2013) 12590–12600.  J.-D. Chai, M. Head-Gordon, Systematic optimization of long-range corrected hybrid density functionals, J. Chem. Phys. 128 (2008). 084106:1-15.  J.-D. Chai, M. Head-Gordon, Long-range corrected hybrid density functionals with damped atom–atom dispersion corrections, Phys. Chem. Chem. Phys. 10 (2008) 6615–6620.  R.S. Mulliken, Electronic population analysis on LCAO-MO molecular wave functions. I, J. Chem. Phys. 23 (1955) 1833–1840.  A.R. Fersht, Catalysis, binding and enzyme-substrate complementarity, Proc. R. Soc. Lond. B 187 (1974) 397–407.
R. Ugarte / Computational and Theoretical Chemistry 1043 (2014) 54–63  M. Pavlov, M.R.A. Blomberg, P.E.M. Siegbahn, New aspects of H2 activation by nickel–iron hydrogenase, Int. J. Quantum Chem. 73 (1999) 197–207.  H. Nakamura, T. Sakamoto, A. Wada, A theoretical study of the dielectric constant of protein, Protein Eng. 2 (1988) 177–183.  G. Löfﬂer, H. Schreiber, O. Steinhauser, Calculation of the dielectric properties of a protein and its solvent: theory and a case study, J. Mol. Biol. 270 (1997) 520–534.  D.M. Lawson, C.E.M. Williams, L.A. Mitchenall, R.N. Pau, Ligand size is a major determinant of speciﬁcity in periplasmic oxyanion-binding proteins: the 1.2 Å resolution crystal structure of Azotobacter vinelandii ModA, Structure 6 (1998) 1529–1539.
 G. Klebe, H.-J. Bohm, Energetic and entropic factors determining binding afﬁnity in protein-ligand complexes, J. Recep. Sig. Trans. Res. 17 (1997) 459–473.  J.D. Chodera, D.L. Mobley, Entropy–enthalpy compensation: role and ramiﬁcations in biomolecular ligand recognition and design, Annu. Rev. Biophys. 42 (2013) 121–142.  M. Korth, A quantum chemical view of enthalpy–entropy compensation, Med. Chem. Commun. 4 (2013) 1025–1033.  D. Roy, M. Marianski, N.T. Maitra, J.J. Dannenberg, Comparison of some dispersion-corrected and traditional functionals with CCSD(T) and MP2 ab initio methods: Dispersion, induction, and basis set superposition error, J. Chem. Phys. 137 (2012). 134109:1-12.