Combination of Virtual Screening Protocol by in Silico toward the Discovery of Novel 4-Hydroxyphenylpyruvate Dioxygenase Inhibitors

4-Hydroxyphenylpyruvate dioxygenase (EC 1.13.11.27, HPPD) is a potent new bleaching herbicide target. Therefore, in silico structure-based virtual screening was performed in order to speed up the identification of promising HPPD inhibitors. In this study, an integrated virtual screening protocol by combining 3D-pharmacophore model, molecular docking and molecular dynamics (MD) simulation was established to find novel HPPD inhibitors from four commercial databases. 3D-pharmacophore Hypo1 model was applied to efficiently narrow potential hits. The hit compounds were subsequently submitted to molecular docking studies, showing four compounds as potent inhibitor with the mechanism of the Fe(II) coordination and interaction with Phe360, Phe403, and Phe398. MD result demonstrated that nonpolar term of compound 3881 made great contributions to binding affinities. It showed an IC50 being 2.49 μM against AtHPPD in vitro. The results provided useful information for developing novel HPPD inhibitors, leading to further understanding of the interaction mechanism of HPPD inhibitors.


INTRODUCTION
The success probability for novel herbicides discovery in agricultural field is cutting down and weeds are becoming widely resistant to most common used herbicide in recent years (Green, 2014), all these make weed management difficult and time consuming. According to the literature, at least 315 weed biotypes, 183 weed species including 110 dicots and 73 monocots worldwide have been reported to have acquired resistance to widely used herbicides (Vishnoi et al., 2009). Thus, there is an urgent need for novel herbicide discovery to overcome the weeds resistance. The 4-hydroxyphenylpyruvate dioxygenase (HPPD) inhibitor offers such solutions by bleaching herbicide mode of action.
HPPD, which was founded by Zeneca Group PLC in 1982, is nonheme Fe(II)-containing dioxygenases (Lee et al., 1997;Neidig et al., 2006). HPPD catalyzes the conversion of p-hydroxyphenylpyruvate (HPPA) to homogentisate (HGA) in aerobic metabolism. This reaction involves decarboxylation, substituent migration and aromatic oxygenation in a single catalytic cycle (Moran, 2005;Purpero and Moran, 2006). The transformation catalyzed by HPPD has both agricultural and therapeutic significance. In plants, HPPD is a key enzyme involving the biosynthesis of the prenylquinones, plastoquinone and tocopherols, and has been used for selective weed control since the early 1990s (Wu et al., 2002). Tocopherol and plastoquinone are produced by further transformation of HGA, and both of them are vital for the natural growth of plants (Yang et al., 2004). In the carotenoid biosynthesis pathway, plastoquinone is an essential cofactor for phytoene desaturase. HPPD inhibition will hinder HPPA-HGA conversion, which gives rise to the deficiency in isoprenoid redox cofactors such as plastoquinone and tocopherol, and finally causes growth inhibition, necrosis and death of treated plants (Borowski et al., 2004;Kovaleva and Lipscomb, 2008;Siehl et al., 2014). In human, the deficiency of tyrosine catabolism enzyme will lead to the tyrosine catabolism pathway suffocation (Raspail et al., 2011).
HPPD is a kind of new herbicide target enzyme, and several HPPD inhibited herbicides have been commercialized. Since the first launch of pyrazolynate by Sankyo in 1980, approximately 13 HPPD inhibitors have been commercialized (Meazzaa et al., 2002;Witschel, 2009). Several of class HPPD inhibitors are currently used as selective broad leaf herbicides including triketones, pyrazoles, isoxazoles, diketone nitriles and benzophenones, among them, the triketone herbicides have contributed to various commercialized HPPD inhibitors through structural modification (Figure S1), such as sulcotrione, mesotrione and benzobicylon (Mitchell et al., 2001;Sutton et al., 2002;Ahrens et al., 2013). These inhibitors also show numerous advantages, such as application security, high activity, low residual, broad-spectrum weeds control (including herbicideresistant weed biotypes), excellent crop selectivity and benign environmental effects (Beaudegnies et al., 2009;Woodyard et al., 2009). Many have been obtained potent inhibitors but none was as potent as sulcotrione. Because of emergence of sulcotrione resistant, therefore obtaining HPPD inhibitors with novel scaffolds is an urgent task for herbicide developers. The inhibitory mechanism has been found in previous studies (Lin et al., 2013;Silva et al., 2015). In HPPD, active site Fe(II) octahedral coordination sphere is accomplished by forming coordinating interaction with three protein ligand atoms and three water molecules. However, two coordinating water molecules are displaced by the 1,3-diketone moiety of the HPPD inhibitor in the enzyme-inhibitor complex. HPPD enzyme shares two His and a Glu residue as Fe(II) ligands and the distance to oxygen atom of His is 2.5 Å (Brownlee et al., 2004). Two other bidentate ligands are from the oxygen of HPPD inhibitors and the distance from oxygen to metal ions is between 1.9 and 2.4. The π-π stacking interaction also plays an important role in the binding model of complex (Ndikuryayo et al., 2017). Based on the above research, a virtual screening method was used to obtain novel HPPD inhibitor. In the present investigation, predictive 3D-pharmacophore model was built on basis of the known HPPD inhibitors reflecting the structure-activity relationship (SAR). Subsequently, the best pharmacophore model was used as a 3D query for searching four databases (Maybridge, Chembridge, ChemDiv, and Specs) to discern novel HPPD inhibitors and also utilized as a predictive program to estimate bioactivity of HPPD inhibitors. Further, molecular docking and molecular dynamics (MD) simulation were performed to identify the most potential HPPD inhibitors with strong ligand binding affinity. The obtained compounds were subsequently assayed the bioactivity in vitro to verify the inhibition on AtHPPD. The detailed screening workflow was shown in Figure 1.

Data Collection and Preparation
This training set, including 16 structurally diverse chemical compounds (Figure 2) from various literatures Cho et al., 2013;Xu et al., 2014;Lei et al., 2016) with wide activity were used to generate 3D-pharmacophore models. The molecular structures of the studied dataset compounds were depicted in SYBYL 6.9 program (SYBYL, Version 6.9 1 ). To validate the hypothesis, the test set 22 compounds obtained from literatures Dayan et al., 2009;Cho et al., 2013;Lei et al., 2016) were prepared using the same protocol (Figure 3).

3D-Pharmacophore Generation Using HypoGen
Based on spatial permutations of pharmacodynamic characteristic elements that play an important role in the activity, the key common chemical features were selected to create 3D-pharmacophore models with HypoGen. The common features of pharmacophore were generated in "Feature Mapping" module by identifying the important chemical features from the training set before establishing hypothesis. Hydrogen bond acceptor (HBA), hydrogen bond donor (HBD), hydrophobic features (HY), and ring aromatic (RA) were selected as key chemical properties to generate the pharmacophore model. Subsequently, pharmacophore models were computed using the "HypoGen" module of DS v3.5 (Catalyst, Version 4.10 2 ) and the top 10 hypotheses based on the fixed cost, null cost, and total cost values were saved. The best hypothesis (Hypo1) model was analyzed and selected according to the cost analysis and consideration of vital factors, namely configuration cost, root mean square deviation (rmsd) and correlation coefficient (r) of the training set between experimental and predicted values.

Validation and Evaluation of Pharmacophore Model
The best ligand-based pharmacophore model evidenced that as rmsd and configuration components are very low and the value of total cost would always be rather higher than null cost but close to the fixed cost. Regard to external validation, firstly the correlation of test set between the experimental and predicted activity was an important index. Secondly randomizing the data using Fischer's Randomization test was performed to obtain cross validation. All process was carried out with the "Ligand Pharmacophore Mapping" module.

Pharmacophore Based on Database Virtual Screening
The well-validated Hypo1 was further used to screen the database consisted of 111560 compounds belonging to four different chemical databases. The four databases were built by "Build 3D Database" protocol of DS v3.5. Subsequently, Chembridge, Maybridge, ChemDiv and Specs were selected as "Input Database" and Hypo1 was imported to "Input Pharmacophore" to discriminate potential HPPD inhibitor from database using "Search 3D Database" of DS v3.5.

Molecular Docking
1917 Compounds which were obtained according to the property of fit value and predictive bioactivity in the previous steps further were used to molecular docking. Docking simulation was carried out utilizing CDOCKER module of the DS v3.5 to explore the binding mode of compounds. The AtHPPD crystal structure (PDB ID: 1TFZ) was downloaded from the Protein Data Bank. The protein was prepared by removing the water, adding hydrogen and correcting the incomplete residues using "Clean Protein" tool in "Prepare Protein" module, then the protein were assigned potentials with CHARMm force field. The active site of protein was predicted and identified using "Edit binding site" module of DS v3.5 according to the native ligand and the radius was set to 10Å. The obtained receptor was used as the "Input Receptor" molecule parameter. All hit compounds subjected to first filtering processes were chosen as "Input Ligand" and docked into the active site of HPPD. The "Pose Cluster Radius" was defined as 0.5 Å for increasing the diversity of the docked poses. The Top Hits was set to 10, which means top the 10 conformations were saved for each ligand based on scoring and ranking by the negative value of CDOCKER energy. The remaining parameters were default. The best binding modes were determined by docking scores and also the comparison with available complex crystal structure of DAS869 with AtHPPD as reference.

Molecular Dynamics Simulations
In order to analyze the ligand-target interactions, further 9 hits compounds with the best docked poses in complex with AtHPPD were submitted to the MD simulation in Amber16 (Case et al., 2017). Meanwhile, mesotrione and 2-(aryloxyacetyl)cyclohexane-1,3-diones were used as positive control for the in silico simulations. The general AMBER force field gaff and ff14SB force field were employed for the ligand and protein, respectively (Wang et al., 2004;Hornak et al., 2006). The 3D structure coordinate files of candidate ligand and control compound were manually edited to match atom number and naming conventions consistent with pdb format for input into the "Antechamber" module of Amber16. The partial atomic charges of ligands were calculated using the AM1-BCC method (Jakalian et al., 2002). The force field of the Fe(II) treated in the "MCPB" module of Amber which was used to build nonbonded model. The nonbonded model with simple form and excellent transferability implemented in the metal center parameter builder (MCPB) tool was employed to treat Fe(II)-protein interaction (Peters et al., 2010;Li and Merz, 2014). The side chain model including Fe(II) coordination sphere with His205, His287, and Glu373 was first created in MCPB module. The geometry optimization and the atomic partial charges of side chain model were calculation through the restrained electrostatic potential (RESP) technique in Gaussian03 (Frisch et al., 2004). Then parameter information included bond, angle, torsion, improper, van der Waals, and electrostatic terms about the Fe(II) coordinating three residues were generated in the "MCPB" module of Amber. The charge neutralized and solvated progress was performed in the "LEaP" module of Amber16. The resulting structure was immersed into a TIP3P water box with 11,614 water molecules in a rectangular periodic box of 10 Å distance around the complex and eight sodium counter ions were added to maintain electro-neutrality of all systems. The energy minimization and equilibration protocol was carried out in the Sander program. First, a minimization with a tightly restrained protein with a force constant of 500 kcal mol −1 Å −2 was applied to all the atoms of the complex to relieve bad contacts in the surrounding solvent. Then, only the protein backbone atoms were fixed with a restraint force of 5.0 kcal mol −1 Å −2 to minimize the side chain and ligand. Finally, a minimization was performed for all atoms with no restraint. In each step, energy minimization was first performed using the steepest descent algorithm for 2,500 steps, and then the conjugated gradient algorithm for another 2,500 steps. The temperature was gradually raised from 0 to 298 K in the canonical (NVT) ensemble with Langevin thermostat to relax the location of the solvent molecules. Then short equilibration with 500 ps was performed to adjust the solvent density in the isothermal isobaric (NPT) ensemble with Monte Carlo barostat. Finally the system was equilibrated for 1 ns without any restraint in the NTP ensemble. All simulations were run with the PMEMD program without any restraints for 10 ns in NTP ensemble (298 K and a pressure of 1 atm) with a 2 fs time step. During the MD simulation, the particle mesh Ewald (PME) algorithm was used to deal with long-range electrostatic interactions, with a cut-off distance of 10Å (Darden et al., 1993;Essmann et al., 1995). And bond lengths involving hydrogen atom were constrained using the SHAKE algorithm (Ryckaert et al., 1977).
The binding free energy ( G bind ) was calculated with the molecular mechanics and Poisson-Boltzmann solvation area (MM/PBSA) methodology was applied based on stable MD trajectory. Entropy was omitted herein because we pay more attention to relative binding free energy for a series of very similar systems.
Where G complex , G receptor , and G ligand are the free energy of complex, receptor and ligand molecules, respectively. The free energy (G) was calculated based on an average over the extracted snapshots from the MD trajectories. Each term can be expressed as follows: Where E gas is gas-phase energy and can be decomposed into internal energy E int , van der Waals energies E vdw and Coulomb E ele ; G sol is salvation free energy; G GB is the polar solvation contribution calculated by solving the GB equation; G np is the nonpolar solvation contribution and was estimated by the solvent accessible surface area (SASA) determined using a water probe radius of 1.4 Å (Yang et al., 2011). The surface tension constant (γ ) was set to 0.0072 kcal mol −1 Å −2 (Sitkoff et al., 1994); TS is entropy term. Free energy decomposition was performed to obtain the key residues contribution to the binding energy of inhibitors. MM-PBSA was used to compute the interaction energies to each residue in Amber16 by considering molecular mechanics energies and solvation energies without considering the contribution of entropies. The binding interaction between residue-ligand pair comprises three terms: the van der Waals contribution ( E vdw ), the electrostatic contribution ( E ele ) and the solvation contribution ( G sol ).

AtHPPD Inhibition Bioassay in Vitro
Recombinant AtHPPD and homogentisate 1,2-dioxygenase (HGD) was derived from Guangfu Yang research team. HPPD inhibitory activity was determined by using method of coupled enzyme in previously published works (Wang et al., 2015). Assays were performed in 96-well plates at 30 • C using a UV/vis plate reader to monitor the generation of maleylacetoacetate at 318 nm. The IC 50 was then calculated based on the plot of the residue activity against different concentration of test compounds at certain concentrations of substrate by fitting the curves. Each experiment was three replicates and averaged. Compound 3881 was purchased from Innochem. Mesotrione was obtained from Wuhan Yuancheng Technology Co. Ltd. and was recrystallized before use.

HypoGen Model Generation for HPPD Inhibitors
The best quantitative pharmacophore models indicated that HBA, HBD, HA and RA features were important. The top 10 pharmacophore hypotheses were generated based on the structure and the activity values of the training set compounds. The values of 10 hypotheses such as total cost, r, and rmsd were statistically significant (  670). Also, the cost difference 122.85 between null and fixed cost was more than 70 bits, which illustrated that total cost was far from the null cost. Comparing with other hypotheses, Hypo1 with higher Cost, better r and low rmsd accounted for all the pharmacophore features and had good predictive ability. Therefore, Hypo1 selected as a best hypothesis was used for further analyses. The Hypo1 chemical features with its geometric parameters were shown in Figure 4A and matched heat map of the ten hypotheses from training set was displayed in Figure 4B. The results indicated that all the training compounds were well matched to the Hypo1 and the compounds with better activity and fitvalue tended to red.

Pharmacophore Model Hypo1 Validation
Verification hypothesis is one of the important processes in the formation of pharmacophore. Test set and Fisher's Randomization test were employed to confirm the quality of pharmacophore.
The test set was used to validate whether the best hypothesis was capable of prediction for the active compounds other than the training set molecules (Debnath, 2002). Hypo1 has predicted the biological activities of most of the test set compounds in the same activity scale with the r being 0.757 (Figure 5). The r for the training set given by Hypo1 was 0.952. This result suggested that the Hypo1 was not only fit for training set compounds but also for the external compounds. It was also verified legalization of Hypo1 to screen databases.   Fischer randomization calculations were performed on the training set to validate statistical robustness of Hypo1 (Sakkiah et al., 2010). 19 Random hypotheses were generated in 95% confidence level and compared with Hypo1, and it was found that the total cost of Hypo1 was lower than that of 19 random models (Figure 6), which meant that the Hypo1 was robust and stable. Meanwhile, statistical r of Hypo1 was far more superior to the randomly generated 19 hypotheses (Figure 7).This result clearly indicated that Hypo1 was not generated occasionally, and the relationship between the structures and bioactivity did exist in the training. Based on these validation results, the best validated Hypo1 as a 3D query was used to retrieve HPPD inhibitors with novel scaffold from four databases.

Database Screening Analysis
Virtual screening based on chemical databases is a fast and accurate approach to identify novel and potential drug (Gogoi et al., 2016;Guedes et al., 2016). The Hypo1 have used as a 3D query tool to screen the chemical databases like Maybridge (20565), Chembridge (40637), ChemDiv (12620), and Specs (37738). According the average fitvalue of the training set was approximately 6.6 and the compounds with a value greater than 6.6 were screened from four databases. As a result, 643, 392, 416, and 466 (a total of 1917) compounds from Maybridge, Chembridge, ChemDiv and Specs, were mapped upon all the features of Hypo1. Screening compounds had scored the HypoGen estimated activity value between 0.01 and 7.3 µM, and thus considered best for further studies.

Molecular Docking Analysis
Molecular docking is a perfect method for predicting the interaction between small molecules and the receptor binding cavity at the molecular level. In the current study, 3D structure of AtHPPD complexed 1TFZ with ligand DAS869 was selected as target protein. In the enzyme-inhibitor complex, the Fe(II) coordinated to the three amino acids and the 1,3-diketone moiety of the DAS869 inhibitor. The distances from the oxygen atoms to the Fe(II) were restrained to a range of 1.9-2.4 Å. This ensured the octahedral geometry and provided a strong ligand orientation and binding force (Yang et al., 2004). Based on the results of the above literature, the native ligand was redocked into corresponding binding pocket. Figure 8 showed that re-docked ligand and natural ligand share the same binding site, whose RMSD values of 0.55 were calculated, confirming the accuracy and feasibility of the docking method of CDOCKER. These 1917 hits obtained through pharmacophorebased screenings were submitted to molecular docking studies by DS v3.5. The distance from Fe(II) to two oxygen atoms of nine candidate ligand was the range of 1.9-2.4 Å. The selection was further obtained 9 compounds ( Table 2) based on their -CDOCKER energy, which picked out 2, 4, 1, and 2 FIGURE 6 | The difference in costs between HypoGen runs and the scrambled runs. The 95% confidence level was selected.
FIGURE 7 | The difference in correlation coefficient between HypoGen runs and the scrambled runs. The 95% confidence level was selected. molecules from Maybridge, Chembridge, ChemDiv and Specs, respectively.
Compared with the binding affinity the known inhibitors, the selected hits were ranked according to the interaction of the amino acid residues at the binding cavity. Most of the sorted molecules from training set and test set showed good interactions with the critical residues like phe360 and phe403. According to this binding modes and binding affinity, four hit compounds ( Table 2) from Chembridge were selected as the target inhibitors. The four compounds selected were structurally similar to the triketones and bind to AtHPPD in the same pose as the triketone. Compared with the training set of compounds, all potent HPPD inhibitors (compound 1-10) containing 2-benzoylethen-1-ol moiety showed better bioactivity, but the activity without the above structure was reduced, such as compound 11-16 (Figure 2) as well as the test set. The common chemical structure of representative HPPD inhibitors ( Figure S1) also contained the 2-benzoylethen-1-ol moiety, which could coordinate to the Fe(II) and interact with phe360 and phe403, as the essential pharmacophore for a significant inhibition of the HPPD activity. The four hit compounds included common subunit of triketone. Two oxygen atoms of hit compounds were responsible for coordinating with Fe(II). It is an effective way to obtain new HPPD inhibitors by lengthening the aryl moiety of triketone compounds according to the literatures (Wang et al., 2016). Aryl in the extended side chain which inserted a C-O between the triketone and side chains formed a more favorable sandwich π-π interaction with residues Phe360 and Phe403 of AtHPPD in the active pocket. Similarly, a C=C was inserted between the triketones and the aryl moiety of the four hit compounds. In addition, an extensive review on the triketone inhibitors indicated that modification of substituent on the benzoyl was a common way to increase inhibitory activity, but transformation of 3-hydroxycyclohex-2en-1-one was hardly seen. The 3-hydroxycyclohex-2-en-1-one part was replaced by 4-hydroxy-6-methyl-2H-pyran-2-one of the obtained four compounds and formed favorable π-π stacking interaction with Phe398, which effectively increased the stability of the ligand binding to AtHPPD. It was also first reported the interaction with Phe398. The predicted pK a for the obtained compounds were less than 6.0, and the weak acidity was favorable for plant uptake and conduction.
An analysis of the co-crystallized HPPD-inhibitor complex showed the interaction pattern between the inhibitor and the HPPD binding site (Figure 9). All selected 4 compounds from Chembridge were inserted well into the binding groove and showed metal-coordination binding to Fe(II). The C-terminal harbored a wide cavity exposed to the solvent that accommodated the Fe(II) cofactor. Ferrous ion coordination was fulfilled by the amino acids His205, His287, and Glu373, while the two remaining coordinating water were replaced by the oxygen atoms of inhibitors. When the distances of two oxygen atoms and the Fe(II) were compared (Figure 10), we found that distance from Fe(II) to two oxygen atoms of 4293 was 2.2 and 2.3Å and the distance of Fe(II) with oxygen atoms of compound 3885, compound 3881 and compound 3882 was 2.3 and 2.3Å. However, the distance of control 1 was 2.4 and 2.4Å. For control 2, the distance was 2.3 and 2.4Å. The decreased distances of the Fe(II) from oxygen atoms of four hit compounds may increase HPPD inhibition.
The detailed chemical interactions of the best hits were presented in Figure 10. The benzene ring of the hit compounds sandwiched by the phenyl of Phe360 and Phe403, forming hydrophobic interaction with the two residues. In addition, hit compounds generated another π-π stacking interactions with Phe398, whereas this interaction was disappeared in the binding of control compounds with HPPD. Therefore, the aromatic subunit should be introduced on the 1,3-dione part for design novel HPPD inhibitors in the future, which lead the π-π interaction with phe398 to increasing the binding affinity. All the interactions indicated that the selected compounds were potential HPPD inhibition.

MD Simulations Analysis
MD simulations were carried out to validate the dynamic interactions between ligands and receptors, and the MM/PBSA program was applied to calculate the binding free energy. To assess the dynamics stability of all the complexes during the MD simulations, root-mean-square-deviations (RMSD) was used to monitor entire MD simulation of each complex. As shown in Figure 11A, all the RMSD values of backbone atoms of protein was very smooth in the whole simulation process maintaining at around 1.5-2.0 Å. Figure 11B showed that the RMSD of the active site of compound 522 suffered a conformational change during the MD simulation after 9 ns. It could be seen that the averaged RMSD of all the heavy atoms of the ligand reached equilibrium at about last 1 ns, which indicated that the all the trajectories was stable during the course of simulation after 1 ns course of simulation in Figure 11C. Therefore, the last 1 ns trajectory was used to analyze the binding free energy and free energy decomposition of all complexes.
For each system, 200 snapshots of each complex were extracted from the last 1 ns stable MD trajectory and used to calculate binding free energy calculations and the results were shown in Table 3. The calculated G bind of the four compounds from Chembridge was higher 30 kcal/mol, however, binding free energy of other compounds from Maybridge, ChemDiv and Specs was lower than control and other compounds in Chembridge. The compound 3885 bound systems remained highest than other bound systems. The E vdw , E ele , and G SA energy were the favorable contribution to G bind , and the positive free energy ( G PB/GB ) displayed adverse effect for all the systems. The E ele made the greatest contribution to the binding free energy for all complexes. And electrostatic energy of control 2 and control 1 was lower than compounds from databases. Therefore, the electrostatic interaction should be increased in modification of the compound to strengthen the interactions between each other. The G PB/GB of the control 1 (91.38. kcal/mol) was higher than that of the control 2 (75.94 kcal/mol) and compound 3881 (59.52 kcal/mol). The probable causes were that the lengthening the aryl side chain of control 2 and compound 3881 reduced the polar solvation energy. It was found that compound 3881(−28.79 kcal/mol) displayed relative lower E vdw energy than those of control 2 (−25.91 kcal/mol). It should be noted that the sum of the nonpolar term ( E vdw + G SA = −32.51 kcal/mol) of 3881 were obviously stronger than that of control 2 with −28.87 kcal/mol. Therefore, these results indicated that nonpolar interactions played a crucial role in the increased binding affinities of 3881 system. The free energy decomposition was calculated to investigate the contribution of key residue for the binding process from the last 1 ns stable MD trajectory and the per-residue contribution for the binding of all systems was plotted in Figure 12. Residues Val207, His287, Phe360, and Phe403 made the greatest contribution to binding energy for control 2. And Residues Phe403 and Phe407 was important composition for binding energy of control 1. His205 was unfavorable for the binding energy of compound 118 and His287 had a negative effect on the binding free energy of compound 120. Residues Leu244, Phe398, and Phe403 had a more than 1 kcal /mol free energy contribution to the binding of compound 3881, 3882, 3885, and    4293 and phe403 made the biggest contribution to the binding between compound 3882 and protein with −2.46 kcal/mol. It can be seen that the residues Phe403 and Phe360 has a greater interaction with control 2, four compounds from Chembridge and two compounds from Maybridge than control 1 and other compounds, and it was found that these compounds included extended side chain of aryl that formed a more favorable sandwich π-π interaction. Furthermore, the contribution of residues Phe398 to the binding of four compounds from Chembridge has an obvious increase, indicating the importance of π-π interaction between Phe398 and candidate ligand 3881, 3882, 3885, and 4293. This made a conclusion that molecular dynamics result verified the interaction with Phe398, Phe403, and Phe360 in molecular docking, which may lead to an increase in hydrophobic interaction. During MD simulation, it is regrettable that hydrogen bond could not be found and these findings are consistent with the previous studies (Moran, 2014) that no hydrogen bonds or ionic interactions contributed to the complex.

Bioactivity Analysis
Compound 3881, identified in silico as promising AtHPPD inhibitor, was submitted to in vitro assays in order to confirm the biological activity. The commercial mesotrione was used as positive control. Figure 13 showed dose-dependent inhibitory activity and the IC 50 values were calculated. The IC 50 of mesotrione was 0.204 µM. The inhibition of compound 3881 at 5 µM was 65.89%, and IC 50 was 2.489 µM, showing good inhibitory activity against AtHPPD in vitro indicating that it may serve as a good template for the design of potent HPPD inhibitor. Furthermore, the activity of compound 3881 was better than compound 10 (IC 50 :11.2 µM) that was parent skeleton of mesotrione in the Figure 2. Compound 3881 was the parent skeleton of all the four hit compounds and exhibited inhibition to AtHPPD, according to analogous to the similarity property principle (i.e., similar chemical structures share similar biological activities) (Johnson and Maggiora, 1990), so it was inferred that the virtual screening derivative 3882 with iodine, 3885 with methoxyl, and 4293 with bromine and hydroxyl also may show good activity of inhibitor. Therefore, compound 3881 will be used as a template compound for structural modification in the subsequent work to obtain higher activity HPPD inhibitors.

CONCLUSIONS
In summary, a reasonable and hierarchical virtual screening workflow was successfully constructed to identify potential HPPD inhibitors. Specially, 3D-pharmacophore and molecular docking were applied to obtain 4 hits from 111560 compounds. Compound 3881(3-cinnamoyl-4-hydroxy-6methyl-2H-pyran-2-one) subjected to biological validation against AtHPPD, showing good inhibitory activity with IC 50 being 2.489 µM. Molecular docking result indicted that hit compounds generated the bidentate coordination with Fe(II) and formed the sandwiched π-π interaction of the benzene ring with Phe403 and Phe360. In addition, it is firstly reported π-π stacking interactions with Phe398. Therefore, the aromatic subunit should be introduced on the 1,3-dione section for design novel HPPD to improve the binding affinity. The MD simulation and MM/PBSA calculations confirmed interaction with Phe398 to make great contributions to the binding free energy, which suggested that the constructed model was reliable and viable. As far as known, these compounds have not been previously reported as HPPD inhibitors. Therefore, this result offered interesting templates for design of novel and more potent HPPD inhibitors.

AUTHOR CONTRIBUTIONS
YF and FY: constructed the workflow; Y-NS and K-HY: built the pharmacophore-based screening work; M-QL and H-FC: performed and analysis the molecular docking; J-ZL: developed the MD experiment; Y-NS: carried out the bioassay of HPPD; YF and FY: discussed and concluded the results; YF: completed the paper.

ACKNOWLEDGMENTS
We are very grateful to the National Nature Science Foundation of China (31772208), the Research Science Foundation in Technology Innovation of Harbin (2015RAYXJ010), and the Nature Science Foundation of Heilongjiang Province (ZD2017002).

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fchem. 2018.00014/full#supplementary-material Figure S1 | Structures of the most representative triketone HPPD inhibitors.