Quantitative Structure Activity Relationship Studies and Molecular Dynamics Simulations of 2-(Aryloxyacetyl)cyclohexane-1,3-Diones Derivatives as 4-Hydroxyphenylpyruvate Dioxygenase Inhibitors

Introduction

4-Hydroxyphenylpyruvate dioxygenase (HPPD), a Fe(II)-dependent non-heme oxygenase, belongs to the α-ketoacid family and plays different roles in organism and plant cells (Rocaboy-Faquet et al., 2014; Huang et al., 2016). It catalyzes the conversion of 4-hydroxyphenylpyruvate (HPPA) into homogentisate (HGA), which is the first committed metabolism of the tyrosine catabolism pathway in humans (Raspail et al., 2011; Moran, 2014; Silva et al., 2015). In plants, HPPD is an essential element in the biosynthesis of plastoquinone and tocopherol; both of them are significant cofactors in the photosynthesis. Inhibition of HPPD will lead to a deficiency of the isoprenoid redox cofactors, followed by the presence of bleaching in plants, eventually bringing necrosis and death (Zou et al., 2007; Wang et al., 2015a).

HPPD has been the subject as an important target for development of new herbicides and multiple series of compounds have been designed and synthesized (Wang et al., 2014, 2015b; Ndikuryayo et al., 2019). When applied pre- or post-emergence, HPPD inhibitors provide control to the important broad leaf weeds in maize and a certain amount of annual weeds (López-Ramos and Perruccio, 2010). HPPD inhibitor herbicides manifest many advantages, for instance, good activity, broad-spectrum weed control, low mammalian toxicity, low residual rate, desired selectivity, and environment friendly (Beaudegnies et al., 2009; Cho et al., 2013; Schultz et al., 2015). However, the first case of HPPD inhibitor herbicide resistance was confirmed in Iowa and Illinois simultaneously in 2010 (Hausman et al., 2011; Kohlhase et al., 2018). Monoculture production systems and multiple uses of herbicides with similar mechanism of action contributed to the generation of weeds resistance to the existing HPPD herbicides (Duke, 2012; Larran et al., 2017; Ye et al., 2018). In response to the evolution of herbicide resistance in weeds, discovering novel inhibitors with high efficiency is urgent. Triketone compounds represent one of the HPPD herbicides, and its substructure is typically based on the 2-benzoyl or 2-heteroaroyl cyclohexane-1,3-dione (Matringe et al., 2005; Roy and Paul, 2010). The activity of triketone HPPD inhibitor was better than any other categories, and they can directly exert effects in the weeds, causing plants to die (Ndikuryayo et al., 2017; Lin et al., 2019).

In this research, a series of 2-(aryloxyacetyl)cyclohexane-1,3-diones derivatives were selected to establish three-dimensional quantitative structure activity relationship (3D-QSAR), applying comparative molecular field analysis (CoMFA) and comparative molecular similarity indices analysis (CoMSIA). Subsequently, molecular docking and molecular dynamics (MD) study was applied to analyze the robustness of the ligands inside the receptor cavity and to learn more about the binding interactions. The analysis strategy is shown in Figure 1. The obtained information will contribute to the rational design of novel HPPD inhibitors with powerful activity in the future.

FIGURE 1

Figure 1. The computational workflow was applied.

Materials and Methods

Data Collection and Preparation

A total of 43 2-(aryloxyacetyl)cyclohexane-1,3-diones derivatives as effective inhibitors were collected to build 3D-QSAR models based on the published literature (Wang D.W. et al., 2016). The activity range of the inhibitors was 0.029–5.571 μM. The structures of these compounds were built and optimized by SYBYL 6.9 to generate 3D structures with appropriate conformation (SYBYL, 2006; Arvind et al., 2014). Simulations were carried out by employing Tripos force field with energy termination of 0.005 kcal/mol, and a maximum of 1,000 iterations. Gasteiger-Hückel charges were used to calculate the partial atomic charges (Zhang et al., 2010).

Molecular Docking

Molecular docking study was applied to obtain corresponding active conformations and analyze receptor–ligand interaction. During the docking operation, 43 HPPD inhibitors were docked into the active pocket of Arabidopsis thaliana HPPD (AtHPPD) using the Accelrys Discovery Studio v3.5 (Catalyst, 2005). The x-ray crystal structure of the AtHPPD (PDB code: 1TFZ) was obtained from the RCSB Protein Data Bank (Yang et al., 2004). All the redundant water molecules and co-crystallized ligand were deleted from the complex before docking study; hydrogen was added to the protein (Yang et al., 2013; Wang D.W. et al., 2016). CHARMm force field was added to the receptor and ligands, and the binding site was defined from the known ligand pose (Fu et al., 2019a,b). Docking operation was performed by CDOCKER protocol with the default docking setting, in which 10 conformations were saved about each ligand based on docking score values (Wu et al., 2003). The postures of the ligands were checked manually, comparing with the co-crystallized ligand (DAS869) in the 1TFZ and other reported inhibitor-enzyme complex crystals (Lin et al., 2019). The chemical structure of the DAS869 is shown in Figure 2. Molecules removed the unreasonable conformations and were used to build the QSAR models. The ligands with the best CDOCKER_ENERGY were employed for the analysis of the binding mode.

FIGURE 2

Figure 2. The chemical structure of DAS869, mesotrione, and compounds 01 and 02.

Alignment of Compounds

Alignment step was extremely important in the process of the development of 3D-QSAR models. The whole data set was divided into training set and test set to develop and validate the model. Random selection is a popular utilized method to build the QSAR models, and the diversity of chemical structures and activities was also considered. Nine compounds were selected as test set, and their structures were abundant, while their pK_i values were uniformly distributed in terms of the value range of the whole set.

To develop an ideal model, two different alignment measures were employed. The first alignment rule was a common framework approach, which appointed the molecule 12 with the best activity as the framework template (Wang J.H. et al., 2016). In this strategy, a multi-search method was applied to search aligning postures with the lowest energy, followed by using the “align database” tool in SYBYL 6.9; all the other compounds were superimposed to the template with the form of common scaffold. Differing from the previous protocol, the second strategy was a receptor-based approach, which states that all molecular conformations were obtained from docking simulation rather than the previous one on the basis of atoms. The best active molecule 12 with docking conformation was chosen as the template molecule.

3D-QSAR Model Generation and Validation

A standard development of CoMFA or CoMSIA model was performed by the partial least squares (PLS) regression analysis to select interrelated components and set up the optimal 3D-QSAR model (Dong et al., 2017). A sp3 carbon atom, as the steric probe, was used with a charge of +1.0 in the process of steric and electrostatic field in CoMFA generation. For CoMSIA analysis, five descriptor fields, namely steric, electrostatic, hydrophobic, H-bond donor, and H-bond acceptor field, were selected to simulate models (Kothandan et al., 2011). pK_i values, which were negative logarithm converted from AtHPPD inhibition K_i values, were carried out as the dependent variable for model development. To establish a model with excellent prediction ability, the leave-one-out (LOO) strategy was used to carry out the cross-validation analysis. The optimum number of components (ONC) was calculated and the cross-validated coefficient (Q²) was obtained to evaluate the model. The model was followed by the non-cross-validation analysis and the coefficient of determination (R²), the standard error of estimate (SEE), and the F value were calculated based on the ONC originated from LOO (Arvind et al., 2014). The predicted correlation coefficient for the test set ($R_{pred}^2$) was used to examine the predictive power of the model. In addition, the reliability and effectiveness of the model were measured through comparing the experimental pK_i values with the predicted pK_i values of the data set.

Molecular Dynamics

A portion of docking complex adopted for MD was designed to explore the major driving force for ligand–receptor interactions and analyze the related amino acid residues. Two representative inhibitors, compound 01 and 02, and commercial herbicide, mesotrione, were selected to perform simulation in the best pose of docking. The chemical structure of these three compounds is shown in Figure 2. Compound 01 represented the backbone of this class of compounds, and the change in the activity of compound 02 was attributed to the introduction of methyl groups at the 5-position of R₁. Mesotrione, a widely used herbicide, was used as a control compound in this study.

The ligand, receptor, and complex information of the two docking structures were introduced in Amber 16 (Case et al., 2016). Initially, the ligands were formatted in Antechamber program and AM1-BCC protocol was employed to calculate the partial atomic charges of molecules. The metal ion, Fe(II), in the protein needed to be specially treated, which was critical to build non-bonding model with simple form and excellent transferability implemented. The metal center parameter builder (MCPB) module of Amber was used to modify Fe(II)-amino acids interaction including His205, His287, and Glu373 (Peters et al., 2010; Li, 2014). The side chain connecting Fe(II) was treated by the restrained electrostatic potential (RESP) tool of Gaussian03. Meanwhile, the atomic partial charges and the geometry optimization were calculated (Frisch et al., 2004). Angle, bond, torsion, improper torsion, van der Waals, and other information parameters were performed through the MCPB. The charge neutralized and solvated progress were generated in the “LEaP” module. In order to produce the appropriate topologic and coordinate files required for the MD simulations, the generalized Amber force field gaff and ff14SB force field were used for ligand and receptor, respectively (Hornak et al., 2006). A rectangular box of TIP3P water was added to the system with a boundary of 10 Å from the edge of the box to the complex atom, and sodium ions that assisted to maintain the electrical properties reflected the neutral state (Gadd et al., 2017). The optimization process was divided into three parts with different constraints. Each section included the steepest descent method of 2,500 steps and the conjugate gradient method of 2,500 steps as well. Heating of the system was a gentle rise in temperature from 0 to 298 K in the canonical (NVT) ensemble with 20 kcal mol⁻¹ Å⁻². A density balance achieved in 500 ps with fixed protein backbone atoms to allow relaxation of the solvent and overall equilibrium lasting 1 ns was performed to ensure the equilibrium state of the MD simulation conditions. Subsequently, the whole simulation was over the course of 10 ns with a 2-fs step.

The procedure of combining free energy calculation was applied to the molecular mechanics method based on all atoms and Poisson-Boltzmann solvation area (MM-PBSA) measure (Hao et al., 2011). The average over the extracted snapshots from the MD stable trajectories was used to compute the free energy. Based on the following equations, the correlative binding free energies were obtained:

$\begin{array}{llll}\Delta G_{bind}& =& G_{cpx}-G_{rec}-G_{lig}& (1)\end{array}$

$\begin{array}{llll}\Delta G_{bind}& =& \Delta E_{MM}+\Delta G_{sol}-T\Delta S& (2)\end{array}$

$\begin{array}{llll}\Delta E_{MM}& =& \Delta E_{int}+E_{ele}+E_{vdw}& (3)\end{array}$

$\begin{array}{llll}\Delta G_{sol}& =& \Delta G_{PB}+\Delta G_{SA}& (4)\end{array}$

where ΔE_MM is determined by the internal energy (ΔE_int) contributed from bonds, angles, and torsions, the van der Waals energy (ΔE_vdw), and electrostatic force (ΔE_ele). ΔGsol denotes the solvation free energy, which consists of the polar solvation contribution (ΔG_PB) and non-polar solvation contribution (ΔG_SA). As the contribution of entropy is insignificant for a series of similar systems, TΔS items are excluded in our study (Fu et al., 2017, 2018).

To obtain the detailed interactions between the HPPD and inhibitors, the binding free energy was decomposed onto each individual residue using the MMPBSA.py module. In the decomposition process, the van der Waals contribution (ΔE_vdw), the electrostatic contribution (ΔE_ele), and the free energy of solvation (ΔG_sol) in the binding process of enzyme and ligands were calculated and the contribution of entropy was omitted.

$\begin{array}{ll}\Delta G_{\text{inhibitor}\_\text{residue}}=\Delta E_{\text{vdw}}+\Delta E_{\text{ele}}+\Delta G_{\text{sol}}& (5)\end{array}$

Results and Discussion

3D-QSAR Models

The framework of molecules, each molecular structure, and the activity values are shown in Table 1. Six statistical parameters including the Q², ONC, R², SEE, F, and $R_{pred}^2$ value are obtained to assess the creditability of each 3D-QSAR model. As far as Q² and R² are concerned, they are considered as two vital standards to evaluate the quality and predictive capability of the QSAR models. In addition, a low SEE value and good F and $R_{pred}^2$ values are also crucial for a reliable model.

TABLE 1

Table 1. The structure of 2-(aryloxyacetyl)cyclohexane-1,3-diones derivatives and corresponding experimental and predicted activity.

The parameters of the obtained models are listed in Table 2. The best CoMFA model based on common framework was established with a best cross-validated correlation coefficient value (Q² = 0.872) and a high conventional correlation coefficient (R² = 0.999). The optimum number of components (ONC) was 10 and the contributions of steric and electrostatic fields were 52.3 and 47.7%, respectively. The standard error of estimate (SEE) was 0.024, the F value was 1776.949, and the predicted correlation coefficient ($R_{pred}^2$) was 0.863, which proved that the model possessed great predictable capability. The CoMFA model based on molecular docking was built with Q² = 0.693 and R² = 0.998, and at this time, the ONC value was 10; the SEE value of 0.034, the F value of 898.323, and the $R_{pred}^2$ value of 0.828 were also obtained. The contribution rate of the steric field was 83.4% and that of the electrostatic field was 16.6%, which suggested that the steric field had more influence on the inhibition than the electrostatic field. The statistical values of the CoMFA model from molecular docking were found to be inferior to those from the common framework, especially on the cross-validated correlation coefficient value (0.693 and 0.872, respectively).

TABLE 2

Table 2. Results of CoMFA and CoMSIA models.

The CoMSIA model, based on common framework, gave a good Q² value of 0.864 and an ideal R² value of 0.990 with 10 components. All the parameters of the model, containing the SEE value of 0.069, the F value of 215.356, and the $R_{pred}^2$ value of 0.850, are shown in Table 2. The model was generated through a combined use of five fields, steric field, electrostatic field, hydrophobic field, hydrogen bond donor, and hydrogen bond acceptor. The contributions were 9.4, 24.2, 28.5, 26.5, and 11.3%, respectively. Electrostatic, hydrophobic, and hydrogen bond donor field had a greater impact on the CoMSIA results, and by modifying these characteristics, the activity would be promoted. The CoMSIA model provided more comprehensive information than CoMFA. The CoMSIA model based on molecular docking led to a satisfactory Q² value of 0.823 using 10 components and an R² value of 0.995 with SEE = 0.050, F = 425.569, and $R_{pred}^2$ = 0.801. The contributions of steric field, electrostatic field, hydrophobic field, hydrogen bond donor, and hydrogen bond acceptor were 9.7, 16.4, 31.4, 30.6, and 11.9%, respectively, which was roughly similar as the proportion of the CoMSIA model from the common framework. The parameters indicated that the CoMSIA models generated by two strategies had both satisfying conventional statistical correlation and good predictive ability of bioactivity.

The plots of the experimental vs. the predicted activity values for all of the compounds are shown in Figure 3. The linear relation between the experimental and predicted pK_i was excellent for either the CoMFA or CoMSIA model from the common framework or molecular docking, indicating closeness of the experimental and predicted biological activity values, and the strong predictive power of the model could be verified. The alignments of the molecules are shown in Figure 4. The molecules used for the common framework alignment were derived from multiple search, which resulted in the conformation of the molecule exhibiting a low energy fold (Figure 4A). The alignment from molecular docking was not the same. During the interaction of the protein and the inhibitors, the conformations of the molecules were in a stretched state (Figure 4B). The conformation of the oxygen atoms used for chelation was highly similar, resulting in a high degree of overlap in this portion. The remaining molecular groups exhibit different postures under the influence of proteins due to their different properties.

FIGURE 3

Figure 3. The alignment of the molecules using (A) the common framework and (B) the docking simulation. Molecules are displayed in white for common C, red for O, blue for N, yellow for S, and green for F and Cl atoms, respectively. For a clear observation, hydrogen atoms are hidden.

FIGURE 4

Figure 4. The plot of experimental and predicted activity based on the common framework (A,B) and molecular docking (C,D).

CoMFA Results

In order to analyze the general feature of the steric and electrostatic field contribution, the structure–activity relationship calculation results of the CoMFA were demonstrated using the contour maps. The steric field result from common framework is shown in Figure 5A; the green color represented that the bulky group was favorable to the bioactivity of the HPPD inhibitors. On the contrary, the less bulky substituent, which was a benefit to the bioactivity, was marked in yellow. Comparing compound 01 with compound 06, it was found that the activity was increased with the change in pK_i values from 5.906 to 6.005 when the hydrogen atom at 4-position of the benzene ring was replaced by methyl. A small yellow area was in the near 5-position of the benzene ring, which suggested that bulky substituents at this site exerted an adverse impact on inhibition. For example, compounds containing a hydrogen atom always displayed better activity than the derivatives (comp. 02, 04, 10, 18, 25, 26, 29, 32, and 33) bearing one or two methyl groups as side chain. The contour map of the electrostatic descriptor based on the common framework is presented in Figure 5B, where the blue region indicated that the electropositive group was favorable to enhance the efficiency of the compounds; in contrast, the red region represented that the electronegative substituents would be conducive to the activity of the compounds. This map meant that the substituents at 2- and 3-position of R₂ would have an electropositive effect, and it was better to have an optimum electronegative action at 4- and 6-position of R₁.

FIGURE 5

Figure 5. CoMFAStDev*Coeff contour maps based on the common framework. (A) Steric contour map. Green and yellow contours show regions where steric bulk has favorable and unfavorable effects on the inhibition ability, respectively. (B) Electrostatic contour map. Blue contours indicate regions where electro-positive groups increase the activity, while white contours indicate regions that were electro-negative.

The CoMFAStDev^*Coeff contour maps, based on molecular docking, are shown in Figure 6 and provided some additional guidance. The steric effects of the substituents need to be increased at the 4- and 5-position of R₂, while the introduction of bulky groups should be avoided. The supplement offered by the electrostatic field was that the 4- and 6-position of R₂ were more suitable for negative groups. For example, compound 42 (R₂ = 2,4,6-tri-Cl, pK_i = 6.735) showed better activity than compounds 01 (R₂ = H, pK_i = 5.906) and 09 (R₂ = Cl, pK_i = 5.942).

FIGURE 6

Figure 6. (A) Steric and (B) electrostatic contours of the CoMFA from molecular docking.

CoMSIA Results

To visualize the generated results, contour maps of CoMSIA based on the common framework is presented in Figure 7. The steric field and electrostatic field of the CoMSIA model based on the common framework (Figures 7A,B) provided the spatial and electrical impact of the substituents on the inhibitor, which were basically similar to the information obtained by the CoMFA contours. Figure 7C depicted the hydrophobic field of CoMSIA, in which white and yellow regions reflected the preference of hydrophilic substitutions and hydrophobic groups. Two white regions at the 3- and 5-position of R₂ symbolized that the addition of the hydrophilic group would enhance the activity; however, introducing a hydrophobic group in the 4-position of R₂ wrapped in yellow would also increase the inhibition, which was supported by compound 41 (R₂ = 3,5-diF-4-CN, pK_i = 7.086) being more active than compound 24 (R₂ = 3,5-diCl, pK_i = 6.020). The hydrogen bond donor is displayed in Figure 7D. In this plot, the cyan displayed positions where a H-bond donor group would be favorable for higher activity. In contrast, purple indicated positions where the H-bond donor of the target molecules is unfavored. There was a cyan contour near the 5-position of the six-membered ring and a small purple contour a little further from the 4-position. The content of Figure 7E showed the effect of the H-bond acceptor on the activity of the molecule, where the magenta and red contours stand for the promotion and suppression of inhibition effect, respectively. The characteristic contours were not at the substituent site, so we could infer that the influence of the H-bond acceptor was minimal to the activity of these series of compounds.

FIGURE 7

Figure 7. CoMSIAStDev*Coeff contour maps based on the common framework. (A) Steric contour map. (B) Electrostatic contour map. (C) Hydrophobic contour map. Yellow and white regions suggest the preference of hydrophobic groups and hydrophilic substitutions, respectively. (D) H-bond donor contour map. Cyan illustrates regions in which the introduction of a H-bond donor group is favored. Purple illustrates regions where the introduction of a H-bond donor group is disfavored. (E) H-bond acceptor contour map. Purple areas are the regions where H-bond acceptor is conducive to the activity; red areas are unfavorable.

The CoMSIA results of molecular docking are shown in Figure 8, and the following discussion focused on the parts that were not obtained previously. The hydrophobic field of CoMSIA gave us a new perspective of a yellow contour with small size surrounding the 5-position of R₁. It suggested that a hydrophobic substituent at this position would increase the inhibitory efficiency. The favorable results of hydrogen bond donors were formed around the hydroxyl group on the six-membered ring, while there were also favorable regions for hydrogen bond acceptors covering the ketone carbonyl of the triketone structure. These results were in line with the actual active data and could prove the accuracy and credibility of our CoMSIA model based on docking.

FIGURE 8

Figure 8. (A) Steric, (B) electrostatic, (C) hydrophobic, (D) H-bond donor, and (E) H-bond acceptor contours of the CoMSIA from molecular docking.

Molecular Docking Analysis

The structure of compounds 01 (pK_i = 5.906) and 02 (pK_i = 5.254) only differed from two methyl groups, but their activities were slightly different. Both were not as active as mesotrione (pK_i = 7.886), which aroused our interest. The overall orientation of these three molecules within the active site pocket of AtHPPD is shown in Figure 9, and it was found that all molecules were fit well into the active cavity. In the process of complexing enzymes and inhibitors, the binding mode of the compound being studied was similar to that of the co-crystallized ligand (DAS869). The three amino acids (His205, His287, and Glu373) involved in chelation with the metal ion remain the same as the co-crystal complex (Yang et al., 2004). The two coordinating water molecules were displaced by different inhibitors. The distances from the 1,3-diketone moiety of the DAS869 inhibitor to the Fe(II) were measured to be 2.3 and 2.4 Å. The chelation distance of compounds 01 and 02 and mesotrione was refined to a range of 2.3–2.4 Å. It is worth noting that Phe360 and Phe403 formed π-π stacking interaction with the benzoyl moiety of DAS869, and similar effects occurred in the benzene of compounds 01 and 02 and mesotrione.

FIGURE 9

Figure 9. Binding model of (A) compound 01, (B) compound 02, and (C) mesotrione.

The conformations of the same part in compounds 01 and 02 were similar. Due to the presence of methyl, the activity of compound 02 was significantly weakened because the two methyls occupied a large pocket space. This inference was consistent with the QSAR results that on the 5-position of R₁, the smaller group was beneficial to increase the inhibitor activity. The docking result of mesotrione showed that the conformation of the six-membered ring was different from that of compounds 01 and 02, and it fitted more closely to the active pocket. The activity of the compounds in this study was lower than that of mesotrione, probably because the oxygen atom in the framework structure affected the conformation of the molecule. The presence of an oxygen atom reversed the six-membered ring of compounds 01 and 02, which, although not affecting its coordination with the iron atom, reduced the activity of the inhibitor. At the same time, the carbon chain was elongated, causing the benzene ring to move back, and π-π interaction was weaker than that of mesotrione. To further explore the factors influencing activity, MD simulations were applied to these three compounds.

MD Analysis

In order to verify whether the systems reached equilibrium during the dynamics simulation, the root-mean-square deviation (RMSD) was calculated, which reflected the dynamic change of the entire structure in the simulation process. RMSD values included the backbone Cα atoms of the protein, active pocket with residues of 5 Å around ligand, and the heavy atoms of ligand (Figure 10). All systems were dynamically changing throughout the kinetics. The RMSD values of the backbone of the compound 02 and mesotrione systems were small, showing higher stability throughout the simulation, and the skeleton structure of compound 01 was more unstable. It is worth noting that mesotrione in the protein complex was not as stable as ligands 01 and 02 during the simulation. Mesotrione entered an unstable phase at 5 ns and eventually restabilized at 9–10 ns. All the RMSDs were steady in the last 1-ns simulation process maintained within the 0.5-Å range. The equilibrium stage of the MD simulation was taken for the binding free energy and free energy decomposition analysis of each compound.

FIGURE 10

Figure 10. RMSD of (A) compound 01, (B) compound 02, and (C) mesotrione.

The calculated results are given in Table 3 including the van der Waals interaction energy (ΔE_vdw), the electrostatic energy (ΔE_ele), the polar solvation free energy (ΔG_PB), the non-polar solvation free energy (ΔG_SA), the interaction energy (ΔE_MM), the solvation contribution (ΔG_sol), and the overall binding free energy (ΔG_bind). It could be seen that the total binding free energies of compounds 01 and 02 and mesotrione were −19.81, −3.65, and −28.34 kcal mol⁻¹, respectively. The calculated binding energy was in good agreement with the experimental activity order. As shown in Table 3, the electrostatic terms occupied the principal driving forces for the three complexes, which made a supreme contribution to the binding free energy. The ΔE_vdw, ΔE_ele, and ΔG_SA calculated by the MM-PBSA approach were the favorable contributions to ΔG_bind; in contrast, ΔG_PB had a certain passive effect. By comparing with systems 01 and 02, the addition of methyl groups at the 5-position of R₁ led to significant distinction in each term and thus its herbicidal activity is poor. It was found that ΔE_ele of compound 02 (−76.36 kcal mol⁻¹) was lower than that of compound 01 (−92.50 kcal mol⁻¹). The unfavorable contribution, ΔG_PB, of compound 02, which was 107.68 kcal mol⁻¹, was stronger than that of compound 01, which was 102.45 kcal mol⁻¹. Interestingly, the change in ΔE_vdw tended to increase the activity of compound 02, and ΔG_SA slightly increased, which had no impact on the overall trend. The ΔE_vdw contribution of mesotrione was significantly greater than 01 and 02, while the inhibition of ΔG_PB was also small; the contribution of ΔE_ele was similar to 02. The binding free energies (ΔG_exp = −RTln K_i) for the compounds were also calculated using the K_i values. The ΔG_exp of compounds 01 and 02 and mesotrione were −8.05, −7.16, and −10.74 kcal mol⁻¹, respectively. We noted that MM-PBSA calculations systematically overvalued the binding free energies between ligand and protein for compound 01 and mesotrione systems. However, the value of ΔG_cal was qualitatively consistent with ΔG_exp, confirming the reliability of MD simulation.

TABLE 3

Table 3. Binding free energy (kcal mol⁻¹) of compounds 01 and 02 and mesotrione^a.

The amino residue contributions of HPPD binding with the ligand at the active site cavity are given in Figure 11. It was generated to understand the binding mechanism of protein–ligand. As listed in the plot, the residue groups including Val207, Leu244, His287, Ala289, Phe371, Glu373, Lys400, and Phe403 participated in the binding with molecule 01. Interestingly, the His205 and Glu373 produced a positive number combination of free energy, which was an unfavorable element even if they were involved in chelation with Fe(II). The candidates promoting contributions to the binding free energy in compound 02 were Val207, Leu244, Pro259, Asn261, His287, Phe360, and Phe403. The residue His205 played the same role as it did in system 01. Residues Val207, Leu244, Leu347, Phe398, Gly399, Phe403, and Phe407 made the greatest contribution to the binding energy for mesotrione, and Glu373 had a negative effect. The contribution of residue Phe403 to the binding of the three compounds obviously increased, which indicated the importance of π-π interaction between Phe403 and candidate ligands. Gly Leu llePhe Pro and Val, which belonged to the nonpolar amino acid family that contributed to nonpolar interactions, act as positive drivers of receptor–ligand binding. ArgAsnGln Glu His and Lys were polar amino acids, which was consistent with the conclusions of binding free energy and electrostatic played a major role in the interaction of molecules with protein.

FIGURE 11

Figure 11. Inhibitor–residue interaction spectrum of complexes.

Conclusion

In the current work, 3D-QSAR models including CoMFA and CoMSIA with ideal cross-validated correlation coefficient values and best correlation coefficient values were established to analyze the 2-(aryloxyacetyl)cyclohexane-1,3-diones derivatives as valid HPPD inhibitors. The structural features conducive to enhance the activity are summarized in Figure 12. The methyl at the 5-position posed an adverse effect on the inhibitor by forming a steric hindrance as well as an effect of oxygen atom in the backbone on the molecular conformation, which was demonstrated by the result of molecular docking. The MD simulation and MM-PBSA energy calculation revealed that the electrostatic energy was the major driving force for ligand binding. It also illuminated the amino acid residues involved in inhibitor–HPPD interaction, in which the Phe403 was prominent in the systems. This study not only is helpful in clarifying the binding mechanism of the HPPD inhibitor but also provides useful information to the discovery of novel HPPD inhibitors.

FIGURE 12

Figure 12. The summary of structure–activity relationship.

Author Contributions

YF and FY conceived and designed the workflow. Y-XL built 3D-QSAR models. K-HY performed molecular docking. M-QL and J-ZL developed the molecular dynamics. YF and Y-XL contributed to the discussion and analysis of the results. YF completed the manuscript. All authors read and approved the manuscript.

Funding

This work was supported by the National Nature Science Foundation of China (no. 31772208), Natural Science Foundation of Heilongjiang Province (ZD2017002), and the Research Science Foundation in Technology Innovation of Harbin (2015RAXXJ032).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

References

Arvind, K., Solomon, K. A., and Rajan, S. S. (2014). QSAR studies on diclofenac analogues as potent cyclooxygenase inhibitors using CoMFA and CoMSIA. Med. Chem. Res. 23, 1789–1796. doi: 10.1007/s00044-013-0771-5

CrossRef Full Text | Google Scholar

Beaudegnies, R., Edmunds, A. J. F., Fraser, T. E. M., Hall, R. G., Hawkes, T. R., Mitchell, G., et al. (2009). Herbicidal 4-hydroxyphenylpyruvate dioxygenase inhibitors-a review of the triketone chemistry story from a Syngenta perspective. Bioorg. Med. Chem. 17, 4134–4152. doi: 10.1016/j.bmc.2009.03.015

CrossRef Full Text | Google Scholar

Case, D. A., Betz, R. M., Cerutti, D. S., Cheatham, T. E. I. I. I., Darden, T. A., Duke, R. E., et al. (2016). AMBER2016. San Francisco, CA: University of California.

Catalyst (2005). Version 4.10. San Diego, CA: Biovia Inc.

Cho, J. E., Kim, J. T., Kim, E., Ko, Y. K., and Kang, N. S. (2013). The structure-based three-dimensional pharmacophore models for Arabidopsis thaliana HPPD inhibitors as herbicide. Bull. Korean Chem. Soc. 34, 2909–2914. doi: 10.5012/bkcs.2013.34.10.2909

CrossRef Full Text | Google Scholar

Dong, H. H., Liu, J., Liu, X. R., Yu, Y. Y., and Cao, S. W. (2017). Combining molecular docking and QSAR studies for modeling the anti-tyrosinase activity of aromatic heterocycle thiosemicarbazone analogues. J. Mol. Struct. 1151, 353–365. doi: 10.1016/j.molstruc.2017.08.034

CrossRef Full Text | Google Scholar

Duke, S. O. (2012). Why have no new herbicide modes of action appeared in recent years? Pest Manag. Sci. 68, 505–512. doi: 10.1002/ps.2333

PubMed Abstract | CrossRef Full Text | Google Scholar

Frisch, M. J., Schlegel, H. B., Scuseria, G. E., Robb, M. A., Cheeseman, J. R., Montgomery, J. A., et al. (2004). Gaussian 03. Wallingford, CT: Gaussian Inc.

Fu, Y., Sun, Y. N., Yi, K. H., Li, M. Q., Cao, H. F., Li, J. Z., et al. (2017). 3D pharmacophore-based virtual screening and docking approaches toward the discovery of novel HPPD inhibitors. Molecules 22:959. doi: 10.3390/molecules22060959

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Sun, Y. N., Yi, K. H., Li, M. Q., Cao, H. F., Li, J. Z., et al. (2018). Combination of virtual screening protocol by in silico toward the discovery of novel 4-hydroxyphenyl pyruvate dioxygenase inhibitors. Front. Chem. 6:14. doi: 10.3389/fchem.2018.00014

CrossRef Full Text | Google Scholar

Fu, Y., Wang, K., Wang, P., Kang, J. X., Gao, S., Zhao, L. X., et al. (2019a). Design, synthesis, and herbicidal activity evaluation of novel aryl-naphthyl methanone derivatives. Front. Chem. 7:2. doi: 10.3389/fchem.2019.00002

PubMed Abstract | CrossRef Full Text | Google Scholar

Fu, Y., Zhang, S. Q., Liu, Y. X., Wang, J. Y., Gao, S., Zhao, L. X., et al. (2019b). Design, synthesis, SAR and molecular docking of novel green niacin-triketone HPPD inhibitor. Ind. Crop. Prod. 137, 566–575. doi: 10.1016/j.indcrop.2019.05.070

CrossRef Full Text | Google Scholar

Gadd, M. S., Testa, A., Lucas, X., Chan, K. H., Chen, W. Z., Lamont, D. J., et al. (2017). Structural basis of PROTAC cooperative recognition for selective protein degradation. Nat. Chem. Biol. 13, 514–525. doi: 10.1038/NCHEMBIO.2329

PubMed Abstract | CrossRef Full Text | Google Scholar

Hao, G. F., Tan, Y., Yu, N. X., and Yang, G. F. (2011). Structure-activity relationships of diphenyl-ether as protoporphyrinogen oxidase inhibitors: insights from computational simulations. J. Comput.-Aided Mol. Des. 25, 213–222. doi: 10.1007/s10822-011-9412-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Hausman, N. E., Singh, S., Tranel, P. J., Riechers, D. E., Kaundun, S. S., Polge, N. D., et al. (2011). Resistance to HPPD-inhibiting herbicides in a population of water hemp (Amaranthus tuberculatus) from Illinois, United States. Pest Manag. Sci. 67, 258–261. doi: 10.1002/ps.2100

CrossRef Full Text | Google Scholar

Hornak, V., Abel, R., Okur, A., Strockbine, B., Roitberg, A., and Simmerling, C. (2006). Comparison of multiple Amber force fields and development of improved protein backbone parameters. Proteins 65, 712–725. doi: 10.1002/prot.21123

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, C. W., Liu, H. C., Shen, C. P., Chen, Y. T., Lee, S. J., Lloyd, M. D., et al. (2016). The different catalytic roles of the metal-binding ligands in human 4-hydroxyphenylpyruvate dioxygenase. Biochem. J. 473, 1179–1189. doi: 10.1042/BCJ20160146

PubMed Abstract | CrossRef Full Text | Google Scholar

Kohlhase, D. R., Edwards, J. W., and Owen, M. D. K. (2018). Inheritance of 4-hydroxyphenylpyruvate dioxygenase inhibitor herbicide resistance in Amaranthus tuberculatus. Plant Sci. 274, 360–368. doi: 10.1016/j.plantsci.2018.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Kothandan, G., Gadhe, C. G., Madhavan, T., and Cho, S. J. (2011). Binding site analysis of CCR2 through in silico methodologies: docking, CoMFA, and CoMSIA. Chem. Biol. Drug Des. 78, 161–174. doi: 10.1111/j.1747-0285.2011.01095.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Larran, A. S., Palmieri, V. E., Perotti, V. E., Lieber, L., Tuesca, D., and Permingeat, H. R. (2017). Target-site resistance to ALS-inhibiting herbicides in Amaranthus palmeri from Argentina. Pest Manag. Sci. 73, 2578–2584. doi: 10.1002/ps.4662

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, P. F. (2014). Taking into account the ion-induced dipole interaction in the nonbonded model of ions. J. Chem. Theory Comput. 10, 289–297. doi: 10.1021/ct400751u

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, H. Y., Yang, J. F., Wang, D. W., Hao, G. F., Dong, J. Q., Wang, Y. X., et al. (2019). Molecular insights into the mechanism of 4-hydroxyphenylpyruvate dioxygenase inhibition: enzyme kinetics, X-ray crystallography and computational simulations. FEBS J. 286, 975–990. doi: 10.1111/febs.14747

PubMed Abstract | CrossRef Full Text | Google Scholar

López-Ramos, M., and Perruccio, F. (2010). HPPD: ligand- and target-based virtual screening on a herbicide target. J. Chem. Inf. Model. 50, 801–814. doi: 10.1021/ci900498n

PubMed Abstract | CrossRef Full Text | Google Scholar

Matringe, M., Sailland, A., Pelissier, B., Rolland, A., and Zink, O. (2005). p-Hydroxyphenylpyruvate dioxygenase inhibitor-resistant plants. Pest Manag. Sci. 61, 269–276. doi: 10.1002/ps.997

PubMed Abstract | CrossRef Full Text | Google Scholar

Moran, G. R. (2014). 4-Hydroxyphenylpyruvate dioxygenase and hydroxymandelate synthase: Exemplars of the α-keto acid dependent oxygenases. Arch. Biochem. Biophys. 544, 58–68. doi: 10.1016/j.abb.2013.10.022

PubMed Abstract | CrossRef Full Text | Google Scholar

Ndikuryayo, F., Kang, W. M., Wu, F. X., Yang, W. C., and Yang, G. F. (2019). Hydrophobicity-oriented drug design (HODD) of new human 4-hydroxyphenylpyruvate dioxygenase inhibitors. Eur. J. Med. Chem. 166, 22–31. doi: 10.1016/j.ejmech.2019.01.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Ndikuryayo, F., Moosavi, B., Yang, W. C., and Yang, G. F. (2017). 4-Hydroxyphenylpyruvate dioxygenase inhibitors: from chemical biology to agrochemicals. J. Agric. Food Chem. 65, 8523–8537. doi: 10.1021/acs.jafc.7b03851

PubMed Abstract | CrossRef Full Text | Google Scholar

Peters, M. B., Yang, Y., Wang, B., Fusti-Molnar, L., Weaver, M. N., and Merz, K. M. (2010). Structural survey of zinc-containing proteins and development of the Zinc AMBER Force Field (ZAFF). J. Chem. Theory Comput. 6, 2935–2947. doi: 10.1021/ct1002626

PubMed Abstract | CrossRef Full Text | Google Scholar

Raspail, C., Graindorge, M., Moreau, Y., Crouzy, S., Lefèbvre, B., Robin, A. Y., et al. (2011). 4-Hydroxyphenylpyruvate dioxygenase catalysis Identification of catalytic residues and production of a hydroxylated intermediate shared with a structurally unrelated enzyme. J. Biol. Chem. 286, 26061–26070. doi: 10.1074/jbc.M111.227595

PubMed Abstract | CrossRef Full Text | Google Scholar

Rocaboy-Faquet, E., Noguer, T., Romdhane, S., Bertrand, C., Dayan, F. E., and Barthelmebs, L. (2014). Novel bacterial bioassay for a high-throughput screening of 4-hydroxyphenylpyruvate dioxygenase inhibitors. Appl. Microbiol. Biotechnol. 98, 7243–7252. doi: 10.1007/s00253-014-5793-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Roy, K., and Paul, S. (2010). Docking and 3D QSAR studies of protoporphyrinogen oxidase inhibitor 3H-pyrazolo[3,4-d] [1,2,3] triazin-4-one derivatives. J. Mol. Model. 16, 137–153. doi: 10.1007/s00894-009-0528-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Schultz, J. L., Weber, M., Allen, J., and Bradley, K. W. (2015). Evaluation of weed management programs and response of FG72 Soybean to HPPD-inhibiting herbicides. Weed Technol. 29, 653–664. doi: 10.1614/WT-D-14-00067.1

CrossRef Full Text | Google Scholar

Silva, T. C., Pires, M. D., de Castro, A. A., de Cunha, E. F. F., Caetano, M. S., and Ramalho, T. C. (2015). Molecular insight into the inhibition mechanism of plant and rat 4-hydroxyphenylpyruvate dioxygenase by molecular docking and DFT calculations. Med. Chem. Res. 24, 3958–3971. doi: 10.1007/s00044-015-1436-3

CrossRef Full Text | Google Scholar

SYBYL (2006). Version 6.9. St. Louis, MO: Tripos Inc.

Wang, D. W., Lin, H. Y., Cao, R. J., Chen, T., Wu, F. X., Hao, G. F., et al. (2015a). Synthesis and herbicidal activity of triketone-quinoline hybrids as novel 4-hydroxyphenylpyruvate dioxygenase inhibitors. J. Agric. Food Chem. 63, 5587–5596. doi: 10.1021/acs.jafc.5b01530

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D. W., Lin, H. Y., Cao, R. J., Ming, Z. Z., Chen, T., Hao, G. F., et al. (2015b). Design, synthesis and herbicidal activity of novel quinazoline-2,4-diones as 4-hydroxyphenylpyruvate dioxygenase inhibitors. Pest Manag. Sci. 71, 1122–1132. doi: 10.1002/ps.3894

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D. W., Lin, H. Y., Cao, R. J., Yang, S. G., Chen, Q., Hao, G. F., et al. (2014). Synthesis and herbicidal evaluation of triketone-containing quinazoline-2,4-diones. J. Agric. Food Chem. 62, 11786–11796. doi: 10.1021/jf5048089

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D. W., Lin, H. Y., He, B., Wu, F. X., Chen, T., Chen, Q., et al. (2016). An efficient one-pot synthesis of 2-(aryloxyacetyl)cyclohexane-1,3-diones as herbicidal 4-hydroxyphenylpyruvate dioxygenase inhibitors. J. Agric. Food Chem. 64, 8986–8993. doi: 10.1021/acs.jafc.6b04110

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J. H., Yang, Y. F., Li, Y., and Wang, Y. H. (2016). A computational study exploring the interaction mechanism of benzimidazole derivatives as potent cattle bovine viral diarrhea virus inhibitors. J. Agric. Food Chem. 64, 5941–5950. doi: 10.1021/acs.jafc.6b01067

PubMed Abstract | CrossRef Full Text | Google Scholar

Wu, G., Robertson, D. H., Brooks, C. L. III, and Vieth, M. (2003). Detailed analysis of grid-based molecular docking: a case study of CDOCKER—a CHARMm-based MD docking algorithm. J. Comput. Chem. 24, 1549–1562. doi: 10.1002/jcc.10306

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, C., Pflugrath, J. W., Camper, D. L., Foster, M. L., Pernich, D. J., and Walsh, T. A. (2004). Structural basis for herbicidal inhibitor selectivity revealed by comparison of crystal structures of plant and Mammalian 4-hydroxyphenylpyruvate dioxygenases. Biochemistry 43, 10414–10423. doi: 10.1021/bi049323o

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, S. G., Hao, G. F., Dayan, F. E., Tranel, P. J., and Yang, G. F. (2013). Insight into the structural requirements of protoporphyrinogen oxidase inhibitors: molecular docking and CoMFA of diphenyl ether, isoxazole phenyl, and pyrazole phenyl ether. Chin. J. Chem. 31, 1153–1158. doi: 10.1002/cjoc.201300449

CrossRef Full Text | Google Scholar

Ye, F., Ma, P., Zhang, Y. Y., Li, P., Yang, F., and Fu, Y. (2018). Herbicidal activity and molecular docking study of novel ACCase inhibitors. Front. Plant Sci. 9:1850. doi: 10.3389/fpls.2018.01850

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L., Tan, Y., Wang, N. X., Wu, Q. Y., Xi, Z., and Yang, G. F. (2010). Design, syntheses and 3D-QSAR studies of novel N-phenyl pyrrolidin-2-ones and N-phenyl-1H-pyrrol-2-ones as protoporphyrinogen oxidase inhibitors. Bioorg. Med. Chem. 18, 7948–7956. doi: 10.1016/j.bmc.2010.09.036

PubMed Abstract | CrossRef Full Text | Google Scholar

Zou, J. W., Luo, C. C., Zhang, H. X., Liu, H. C., Jiang, Y. J., and Yu, Q. S. (2007). Three-dimensional QSAR of HPPD inhibitors, PSA inhibitors, and anxiolytic agents: Effect of tautomerism on the CoMFA models. J. Mol. Graph. 26, 494–504. doi: 10.1016/j.jmgm.2007.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar