Interaction between glyphosate and mitochondrial succinate dehydrogenase

Interaction between glyphosate and mitochondrial succinate dehydrogenase

Computational and Theoretical Chemistry 1043 (2014) 54–63 Contents lists available at ScienceDirect Computational and Theoretical Chemistry journal ...

2MB Sizes 0 Downloads 70 Views

Computational and Theoretical Chemistry 1043 (2014) 54–63

Contents lists available at ScienceDirect

Computational and Theoretical Chemistry journal homepage:

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 [1]. The most of the genetically modified 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. Scientific research discommends their use due to possible adverse effects on the environment and on human health [4]. 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 [14]. Cell death was estimated by inhibition of succinate dehydrogenase activity and thus of mitochondrial metabolism [16], necrosis by release of cytosolic adenylate kinase measuring membrane damage and apoptosis via ⇑ Tel.: +56 63 2293133. E-mail address: [email protected] 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 first glance, this is a logical conclusion because at physiological pH the chemical species is mainly the glyphosate dianion [17], 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 [18] 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 flavoprotein (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 first 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 figure 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 thenoyltrifluoroacetone and carboxin, binds to the ubiquinone active site; a hydrophobic pocket comprised of residues from SDHB, SDHC and SDHD subunits [19] 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 [19]. 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) [24]. Fortunately, this structure contains a co-crystallized malonate (Fig. 4), which simplifies 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 defined 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 efficient algorithms, molecular modeling can be effectively used to solve complex chemical and biological problems [18]. 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-Trifluoromethyl-Benzamide, determined at 3.24 Å resolution and retrieved from the Protein Data Bank (PDB: entry code 3ABV) [24]. The construction of the succinate binding site model was carried out with the ArgusLab molecular modeling software [25]. 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 file. 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 [26] 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 configurations. The FAD moiety was ‘‘pruned’’ to reduce the computational cost (Fig. 6c). This ‘‘FAD fragment’’ (isoalloxazine tricyclic ring derivative) along with the un-ionized configurations 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 [27], 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 [28]. Initially, the models were submitted to a PM6 semiempirical [29] geometry pre-optimization, in which the original coordinates of the PDB file 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 [30] 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 [33]. Gaussian provides a reasonable alternative to incorporate this effect, because allows that the system be studied via the self-consistent reaction field (SCRF) approach, in which the medium is represented by a polarizable continuum model (PCM) [34] 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 file. 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 [35]. 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 spacefill model. The figure 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 first to construct the in vacuo ONIOM potential energy surface (PES) [36] 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 field of the fixed 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 [37], 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 fixed 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 influence, which is reflected 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 fixed 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 field (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 five 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 significantly. 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




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 field (SCRF) calculations were carried out with the polarizable continuum model (PCM). The more stable structures of glyphosate (G) are H-bonded. Structure





7006.6526 7006.6773 7006.6675 7006.6443 7006.6447 6534.1919 RHF//PM6 6119.7429 RHF//RHF 886.8863 414.4373







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 [42] 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

In vacuo

H(16) O(5) O(15)

0.739 –0.572 –0.788

RHF/6-31+G(d) e

In vacuo




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 [43]. 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 significant, 0.12 Å and 0.066 Å, malonate distorts less the BS and turns out to be less modified. 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 figures 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 fitted 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





2.03 1.83





2.02 1.94

1.95 1.95

1.92 2.03

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 specificity. 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 [47]. 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 affinity of a ligand toward its binding site falls into a range of approximately 2–20 kcal mol1, and it is




D [DEc]a

M06 wB97xD

74.72 73.77

47.07 42.88

27.65 30.89


related to the binding constant and corresponds to the Gibbs free energy difference for the process. Consequently, enthalpic and entropic effects determine the binding affinity [48]. 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 fulfilled [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 confirmed 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 [51]. 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 fixed binding site was calculated. Five minimum energy conformers of glyphosate in the molecular field 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 first 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 confirm 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 financial support. References [1] 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. [2] 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. [3] 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. [4] 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. . [5] 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. [6] 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. [7] 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. [8] 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. [9] 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. [10] M. Hultberg, Cysteine turnover in human cell lines is influenced by glyphosate, Environ. Toxicol. Pharmacol. 24 (2007) 19–22. [11] 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. [12] 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. [13] 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. [14] 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. [15] 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. [16] F. Peixoto, Comparative effects of the roundup and glyphosate on mitochondrial oxidative phosphorylation, Chemosphere 61 (2005) 1115– 1122.

[17] 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. [18] A.R. Leach, Molecular Modelling: Principles and Applications, second ed., Prentice Hall, 2001. [19] R. Horsefield, 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. [20] K.S. Oyedotun, B.D. Lemire, The quaternary structure of the saccharomyces cerevisiae succinate dehydrogenase, J. Biol. Chem. 279 (2004) 9424–9431. [21] W.C. Kenney, The reaction of N-ethylmaleimide at the active site of succinate dehydrogenase, J. Biol. Chem. 250 (1975) 3089–3094. [22] 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. [23] O.S. Hajjawi, Succinate Dehydrogenase: Assembly, Regulation and Role in Human Disease, Eur. J. Sci. Res. 51 (2011) 133–142. [24] 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. . [25] M.A. Thompson, ArgusLab 4.0, Planaria Software, LCC, Seattle, WA. . [26] R. Dennington, T. Keith, J. Millam, K. Eppinnett, W.L. Hovell, R. Gilliland, GaussView, Version 3.09, Semichem Inc, Shawnee Mission, KS, 2003. [27] S. Joy, P.S. Nair, R. Hariharan, M.R. Pillai, Detailed comparison of the proteinligand docking efficiencies of GOLD, a commercial package and ArgusLab, a licensable freeware, Silico Biol. 6 (2006) 601–605. [28] 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. [29] J.J.P. Stewart, Optimization of parameters for semiempirical methods V: Modification of NDDO approximations and application to 70 elements, J. Mol. Model. 13 (2007) 1173–1213. [30] T. Vreven, K. Morokuma, Hybrid methods: ONIOM(QM:MM) and QM/MM, Ann. Rep. Comput. Chem. 2 (2006) 35–51. [31] D.R. Hartree, The Calculation of Atomic Structures, John Wiley & Sons Inc., New York, 1957. [32] R.G. Parr, W. Yang, Density Functional Theory of Atoms and Molecules, Oxford University Press, Oxford, 1989. [33] F. Himo, Quantum chemical modeling of enzyme active sites and reaction mechanisms, Theor. Chem. Acc. 116 (2006) 232–240. [34] 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. [35] 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 fields, J. Phys. Chem. 98 (1994) 11623–11627. [36] 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. [37] L. Koziol, N. Kumar, S.E. Wong, F.C. Lightstone, Molecular recognition of aromatic rings by flavin: electrostatics and dispersion determine ring positioning above isoalloxazine, J. Phys. Chem. A 117 (2013) 12946– 12952. [38] 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. [39] M. Walker, A.J.A. Harvey, A. Sen, C.E.H. Dessent, Performance of M06, M06-2X, and M06-HF density functionals for conformationally flexible 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. [40] J.-D. Chai, M. Head-Gordon, Systematic optimization of long-range corrected hybrid density functionals, J. Chem. Phys. 128 (2008). 084106:1-15. [41] 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. [42] R.S. Mulliken, Electronic population analysis on LCAO-MO molecular wave functions. I, J. Chem. Phys. 23 (1955) 1833–1840. [43] 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 [44] 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. [45] H. Nakamura, T. Sakamoto, A. Wada, A theoretical study of the dielectric constant of protein, Protein Eng. 2 (1988) 177–183. [46] G. Löffler, 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. [47] D.M. Lawson, C.E.M. Williams, L.A. Mitchenall, R.N. Pau, Ligand size is a major determinant of specificity in periplasmic oxyanion-binding proteins: the 1.2 Å resolution crystal structure of Azotobacter vinelandii ModA, Structure 6 (1998) 1529–1539.


[48] G. Klebe, H.-J. Bohm, Energetic and entropic factors determining binding affinity in protein-ligand complexes, J. Recep. Sig. Trans. Res. 17 (1997) 459–473. [49] J.D. Chodera, D.L. Mobley, Entropy–enthalpy compensation: role and ramifications in biomolecular ligand recognition and design, Annu. Rev. Biophys. 42 (2013) 121–142. [50] M. Korth, A quantum chemical view of enthalpy–entropy compensation, Med. Chem. Commun. 4 (2013) 1025–1033. [51] 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.