In Silico Design and Validation of a Novel HPPD‐Inhibiting Herbicide Candidate Based on Benzofuran and Arylthioacetic Acid Scaffolds
Abstract
Inhibition of 4‐hydroxyphenylpyruvate dioxygenase (HPPD) is a well‐established strategy for weed control, yet the emergence of resistance underscores the need for more potent inhibitors. In this study, datasets of benzofuran analogues and arylthioacetic acid–derived triketones were analyzed using multivariate image analysis of quantitative structure–activity relationships (MIA‐QSARs) to guide the identification of promising candidates. Molecular docking was conducted on HPPD, including Co(II) substitution to explore metal–ligand interactions, and molecular dynamics simulations, parametrized with MCPB.py, evaluated the stability and interaction patterns of the complexes. Fourteen candidates were proposed, six of which exhibited higher predicted activity and improved performance relative to mesotrione. Among them, candidate P1 emerged as the most promising, reproducing key interactions of mesotrione while displaying lower binding energy and stable convergence during dynamics. These results demonstrate P1 as a novel HPPD inhibitor and highlight the utility of combining MIA‐QSAR, docking, and MCPB.py–parametrized dynamics in rational metalloenzyme inhibitor design.
Article type: Research Article
Keywords: MIA‐QSAR, molecular docking, molecular dynamics, HPPD, benzofuran
Affiliations: Department of Chemistry Institute of Natural Sciences Federal University of Lavras Lavras Brazil
License: © 2025 The Author(s). Chemistry & Biodiversity published by Wiley‐VHCA AG. CC BY 4.0 This is an open access article under the terms of the http://creativecommons.org/licenses/by/4.0/ License, which permits use, distribution and reproduction in any medium, provided the original work is properly cited.
Article links: DOI: 10.1002/cbdv.202503221 | PubMed: 41247843 | PMC: PMC12761364
Relevance: Moderate: mentioned 3+ times in text
Full text: PDF (3.2 MB)
Introduction
The hydroxyphenylpyruvate dioxygenase (HPPD) active site is a key target in the search for safe and selective herbicides [ref. 1, ref. 2, ref. 3, ref. 4]. Triketone compounds, designed to structurally compete with the natural substrate hydroxyphenylpyruvate, have been developed following the observation of the allelopathic effects of leptospermone on weeds surrounding the tree Melaleuca citrina, with sulcotrione being the first commercially recognized and patented compound in 1991 [ref. 5, ref. 6, ref. 7].
From a biological perspective, HPPD catalyzes the reaction that produces homogentisate [ref. 8]. This compound, as a precursor in the biosynthesis of plastoquinone and α‐tocopherol, is essential for plants because it participates directly in carotenoid synthesis, reactive oxygen species scavenging, and electron transport in the photosynthetic chain [ref. 1, ref. 7, ref. 9, ref. 10]. Consequently, its depletion causes severe plant damage, primarily characterized by chlorophyll destruction, which has led to the classification of HPPD inhibitors as bleaching herbicides [ref. 3, ref. 5, ref. 11, ref. 12].
However, the growing resistance of weeds to commercial compounds has driven research toward the development of new structures capable of inhibiting HPPD and, consequently, homogentisate synthesis [ref. 13]. Among the previously susceptible weed species in which HPPD inhibitors were once effective, seven have already been reported to have developed resistance [ref. 14]. This finding justifies the design of new, optimized molecules to overcome current limitations, such as the need for higher application rates and the use of undesirable combinations. Thus, the use of quantitative structure–activity relationship (QSAR) methods emerges as a promising strategy to guide the design of new candidate molecules, which, within a defined set of analogous and diverse structures, may exhibit improved inhibitory activity by exploring combinations of substituents reported in the literature [ref. 15].
In this context, a multivariate image analysis of QSAR (MIA‐QSAR) approach is particularly attractive due to the ease of descriptor generation and its low computational cost. In addition to its simplicity, enhancements that account for Van der Waals radii and color diversity in atom representations allow periodic properties to be considered for each element type [ref. 16], providing a more robust analysis than older binary descriptor methods, which could not incorporate additional information [ref. 17, ref. 18]. Beyond model calibration, MIA plots complement this approach by visually representing the statistical results, thereby facilitating the interpretation of group performance and supporting the design of new molecules with enhanced activity [ref. 18].
The requirements for an effective HPPD inhibitor are as previously described: a conserved chelation region, which mimics the α‐keto acid substrate moiety to enable competitive inhibition, and a binding site on an aromatic moiety, which exhibits higher variability among the set of possible interactions [ref. 19, ref. 20]. Among all these possible compounds, different triketone HPPD hybrids are able to be combined for QSAR‐based optimization, due to their ketonic congeneric moiety, as shown in Figure 1.

Despite the stacking interactions promoted by the aromatic ligand moieties, electronic effects primarily govern the enhancement of these compounds’ bioactivity. The introduction of electron‐withdrawing groups on the aromatic ring is frequently recommended, as they increase the α‐carbonyl pKa, thereby promoting the formation of the enolic tautomer, which represents the active inhibitory form [ref. 21, ref. 22, ref. 23]. In addition, such substituents reduce the electron density of the aromatic rings, thereby facilitating stronger π–π interactions with phenylalanine residues [ref. 20, ref. 24, ref. 25].
Understanding these key interactions has encouraged researchers to explore the influence of different aromatic ring systems, aiming to optimize π–π stacking interactions among diverse structural frameworks. Given the well‐established bioactive properties of benzofuran and arylthioacetic acid derivatives, some studies introduced structural modifications at the triketone moiety to identify new hybrid triketone herbicides [ref. 4, ref. 26, ref. 27, ref. 28, ref. 29]. As a result, several active analog inhibitors were discovered in vitro, providing a valuable foundation for subsequent biological and computational studies.
In the present study, both compound sets were analyzed using the MIA‐QSAR approach, and new candidate molecules were proposed based on the obtained results. Rigid molecular docking was then performed to investigate potential enzyme–ligand interactions, followed by molecular dynamics simulations to validate the stability of these formed complexes.
Methods
MIA‐QSAR Procedure
Two datasets comprising 59 congeners derived from benzofuran and arylthioacetic acid [ref. 4, ref. 28] were merged and analyzed using the MIA‐QSAR method. In the Y block, the inhibitory activity values, measured by in vitro assays in the original works, were converted to pK i values (–log K i) [ref. 30, ref. 31, ref. 32] (Table 1). The selection of both scaffold molecules was guided by specific methodological requirements. First, the biological assay procedures reported in the original studies were comparable, and the reference compound exhibited similar activity values across both datasets, ensuring consistency in the bioactivity measurements. Second, the previously demonstrated inhibitory potential of these compounds confirmed their feasibility as biologically active candidates. Finally, the presence of a congeneric moiety within the molecular structures was considered essential for MIA‐QSAR modeling, allowing proper alignment and descriptor correlation.
TABLE 1: Triketone derivatives used in the calibration step of quantitative structure–activity relationship (QSAR) modeling and their corresponding pK i values (K i in mol/L).
The chemical structures for the X‐block (descriptors) were sketched in GaussView, ensuring consistent alignment across all compounds [ref. 33]. Shadowing effects were disabled, and the structures were saved in MOL2 format. Pictorial representations of each compound were obtained as pixel values (432 × 300 pixels), ranging from 0 to 765 on the RGB color scale. These values were then modified by incorporating periodic properties such as electronegativity (ε), Van der Waals radius (rVdW), and the ratio rVdW/ε, and paired with inhibitory activities using Chemoface software [ref. 34]. The correlation between descriptors and responses was modeled using partial least squares (PLS) regression, and William´s plots were employed for outlier identification. Parameters from cross‐validation and Y‐randomization, such as r 2, q 2, c r 2 p, and their respective root mean square error (RMSE) values, were used to assess model quality [ref. 35, ref. 36, ref. 37, ref. 38].
External validation was performed using a bootstrapping methodology, in which 25% of the samples (totaling 14 samples for each cycle) were randomly separated in each of 10 modeling iterations, in which the 14 removed samples, based on which iteration, can be found in the Supporting Information. Validation parameters at this stage included r 2 pred, Avg. r 2 m, Δr 2 m (Roy’s parameters) [ref. 27], and the concordance correlation coefficient (CCC) [ref. 39, ref. 40]. VIP and b‐plots, related to variance and PLS regression coefficients, were generated and used to propose new potential structures, whose activities were predicted with the previously developed models [ref. 18]. Finally, a synthetic route, inspired by the original dataset article, was proposed for the most active molecule obtained to illustrate the feasibility of its synthesis.
Molecular Docking Protocol
To evaluate the correlation between the inhibitory activity of the proposed molecules and their interaction energies within the enzyme cavity, molecular docking was performed with the HPPD enzyme from Arabidopsis thaliana (PDB ID: 5YWG; EC 1.13.11.27) [ref. 41]. The complex contained a cobalt (II) cofactor and the ligand mesotrione, which was also used in the QSAR dataset as the standard compound. For the removal of water molecules and ligand separation, Discovery Studio software was employed [ref. 42]. The ligand, saved in MOL2 format, was processed by adding hydrogens, calculating partial charges with the Gasteiger force field, and correcting atom hybridization and bond types.
The redocking procedure, as protocol validation, was carried out with AutoDockTools [ref. 43]. The enzyme (apo form, in PDB format) had polar hydrogens and partial charges added using the Gasteiger method. Cofactor parameters were set according to the user manual, and the file was saved in PDBQT format. Using a text editor, the charge of cobalt was manually adjusted to +2 to reflect its oxidation state.
The ligand, also converted to PDBQT format after identifying active torsions, was docked with a Genetic Algorithm search method, using a rigid protein structure. The binding cavity was defined based on the crystallized ligand position, with the grid box centered at coordinates (24.434, ‐7.429, ‐29.328) for the x, y, and z axes, respectively. To encompass the entire cavity region, 40 grid points were used in each dimension with a grid spacing of 0.375 Å.
The docking parameters included an initial population of 150 individuals, a mutation rate of 0.02, a crossover rate of 0.8, and a maximum of 27,000 generations. A total of 100 independent runs were conducted. The quality of the docking poses—and thus the reliability of the protocol—was evaluated using the Docking Accuracy (DA) equation, based on the root mean square deviation (RMSD) relative to the crystallized reference structure [ref. 44]. This analysis identifies poses with the smallest deviations from the experimental conformation, with DA values closer to 1 indicating better performance. Equation 1 expresses this relationship, where f corresponds to the fraction of poses within a defined RMSD range, and l and h denote the lower and upper RMSD limits, respectively (l < h) [ref. 44]. The adopted l and h values were 2 and 3 Å, consistent with those commonly used in the literature [ref. 45].
For the proposed ligands, the same preparation and docking protocol was applied, except for an additional structural minimization performed using the DREIDING force field in Discovery Studio software.
Molecular Dynamics
For molecular dynamics simulations, the MCPB.py tool in AmberTools24 was employed [ref. 46, ref. 47, ref. 48]. Based on the optimized mesotrione–HPPD complex, a hybrid approach was used, in which only the amino acid residues directly bound to the metal ion were included [ref. 49]. For the bonded‐system parametrization in Gaussian16, the cobalt (II) ion was coordinated with histidine and glutamate residues, and the standardized ligand was used to preserve octahedral geometry [ref. 50]. Calculations were performed at the B3LYP/6‐31G level, after which ligand–metal bonds were removed and the system was converted to GROMACS parameters [ref. 51].
Four simulation steps were performed in GROMACS: (i) energy minimization using the steepest descent method; (ii) NVT equilibration (–DPOSRES, 100 ps, dt = 2 fs) to stabilize temperature with positional restraints; (iii) NPT equilibration, under the same conditions, to adjust system density and pressure; and (iv) production dynamics without restraints, consisting of 100 ns simulated over 50 000 000 steps (dt = 2 fs). From the resulting trajectories, the RMSD of the ligand within the enzyme cavity was calculated, as well as root mean square fluctuation (RMSF) values for all residues. B‐factors were also used as a pictographic tool to visualize structural flexibility in both the enzyme and ligands [ref. 52].
Physicochemical Properties and ADMET Analysis
For the human safety evaluation, aimed at ensuring low health risks of the proposed compounds, ADMET predictions were performed using the ADMETlab webserver [ref. 53]. The two most promising proposed molecules, in comparison with mesotrione, were assessed based on commonly employed ADMET parameters reported in previous studies: solubility (log S), absorption (Caco‐2), distribution (BBB and PPB), and toxicity, including hepatotoxicity, skin sensitization, and carcinogenicity [ref. 2, ref. 54, ref. 55, ref. 56, ref. 57]. Finally, the log P value for each proposed molecule was predicted with the Molinspiration web server, a reliable predictor for this parameter [ref. 58].
Results and Discussions
MIA‐QSAR Procedure
The MIA‐QSAR method applied in this study is advantageous because it reduces both operational and computational costs while providing valuable interpretative insights through MIA plots. The resulting models typically exhibit high accuracy and predictive power, further supported by docking analyses that confirm the correlation between inhibitory activity and the binding energies of the predicted poses within the biological target.
The dataset was monitored to identify potential outliers. Of the total samples, 59 were deemed suitable for analysis, as they had defined inhibitory activity values (K i). The most suitable models were defined as those showing the lowest RMSE in the cross‐validation procedure. Within this subset, two samples (B19 and B110) showed significant influence in the Williams plot, as evaluated by the threshold of 2.5 for p in the studentized residual analysis (see Supporting Information). The decision to remove these two samples was based on improvements in the statistical parameters adopted and on the increased information content in both the X and Y blocks.
In the developed models, the PLS method demonstrated consistent performance across the three variables used as descriptors for block X. Performance evaluation included the goodness of fit (r 2), the correlation between dependent variables (r 2 y‐rand and c r 2 p), the predictive ability for the calibration set (q 2), and the accuracy in predicting the behavior of external samples (r 2 pred, r 2 m, CCC). Acceptable values were obtained for all three adopted properties, as summarized in Table 2. Among the models generated for each periodic property, the best‐performing model was defined as the one with the highest average values above 0.5. Accordingly, Model 7 (from 10 bootstrap runs) was identified as the most suitable. Its MIA plots, shown in Figure 2, allow a clear distinction of the groups through visualization of the superimposed sample images.
TABLE 2: Statistical parameters to multivariate image analysis of quantitative structure–activity relationship (MIA‐QSAR) validation obtained from the bootstrapping procedure.
| Parameter | rvdW | ε | rvdW/ε | Average | Std. Dev. | Cut‐off |
|---|---|---|---|---|---|---|
| PLS comp. | 6.60 | 5.30 | 7.80 | 6.5667 | 1.2503 | — |
| RMSEC | 0.1725 | 0.1802 | 0.1681 | 0.1736 | 0.0061 | — |
| r2 | 0.9022 | 0.8932 | 0.9057 | 0.9004 | 0.0064 | ≥ 0.6 [ref. 38, ref. 59 ] |
| RMSEy‐rand | 0.3772 | 0.3886 | 0.3695 | 0.3784 | 0.0096 | — |
| r2y‐rand | 0.5268 | 0.5017 | 0.5487 | 0.5258 | 0.0235 | [r 2 y‐rand << r 2] [ref. 38] |
| cr2p | 0.5788 | 0.5908 | 0.5679 | 0.5792 | 0.0115 | ≥ 0.5 [ref. 35] |
| RMSECV | 0.3663 | 0.3634 | 0.3594 | 0.3630 | 0.0035 | — |
| q2 | 0.5819 | 0.5868 | 0.5939 | 0.5875 | 0.0060 | ≥ 0.5 [ref. 38, ref. 39, ref. 59] |
| RMSEP | 0.2555 | 0.2399 | 0.2731 | 0.2561 | 0.0166 | — |
| r2pred | 0.7315 | 0.7676 | 0.7098 | 0.7363 | 0.0292 | ≥ 0.5 [ref. 39] |
| Avg. r 2 m | 0.6006 | 0.6573 | 0.5697 | 0.6092 | 0.0445 | ≥ 0.5 [ref. 39, ref. 40] |
| Δr 2 m | 0.1384 | 0.1196 | 0.1279 | 0.1286 | 0.0094 | < 0.2 |
| CCC | 0.8326 | 0.8600 | 0.8157 | 0.8361 | 0.0223 | ≥ 0.8 [ref. 60] |

The overall analysis of the MIA plots revealed significant variance and strong correlations between substituents and inhibitory activity, particularly at the R3 position. The –OCH3 and –Cl substituents were identified as the most effective contributors to increased activity, consistent with previous reports highlighting the superior performance of electron‐withdrawing groups directly attached to the aromatic ring adjacent to the triketone carbonyl. In contrast, alkyl substitutions at the R1 position produced undesirable effects, as also noted in earlier studies, due to an unfavorable increase in steric hindrance. Other positions showed beneficial, though less pronounced, effects. Consequently, –Cl and –F at R2 and –CH3 at R4 were also considered in the design of new structures, as summarized in Table 3.
TABLE 3: Proposed structures based on the multivariate image analysis (MIA) plot from quantitative structure–activity relationship (QSAR) analysis.
As shown in Table 3, compounds P1 and P2 had the highest predicted activity values, followed by P8 and P9. Moreover, 6 of the 14 proposed molecules displayed higher predicted activity than the reference compound mesotrione (pKi = 7.699), with 9 of them surpassing the activity of all other compounds. Based on the original Wang et al. report [ref. 4], the synthetic route for P1, the molecule with the highest predicted activity, was developed to provide guidance for its future synthesis. The proposed route is illustrated in Figure 3.

Molecular Docking Protocol
For the docking methodology, the selected enzymatic complex was used in the redocking procedure, where the crystallographic ligand served as the reference in activity assays. To evaluate the obtained poses, the Docking Accuracy (DA) equation was applied as described in the methodology section. In the first range, defined by an RMSD of up to 2 Å from the crystallized pose, 83 poses were found, while the remaining 17 were located within the second range, between 2 and 3 Å. Given that all poses were distributed within these predefined ranges and considering the high DA value obtained (DA = 0.915), the adopted docking methodology was considered suitable.
Among all compounds, representative samples were selected based on their activity performance and interaction energy profiles. These included: the best compounds from each original dataset (A7 and B9), the least active sample (B1), the best proposed compound (P1), and the crystallized ligand (Mesotrione), along with its best redocked pose (69th pose). Their binding energy components are summarized in Table 4. The free binding energy was calculated as the sum of all interaction energies, excluding the internal energy. In this context, torsional and internal terms describe ligand intramolecular contributions, while electrostatic and van der Waals terms represent intermolecular interactions with the enzymatic cavity.
TABLE 4: Decomposition of binding free energy values obtained from AutoDock 4 for the most favorable pose of each compound (kcal/mol).
| Samples | ΔG bind. | Electro. | VdW | Internal | Torsional |
|---|---|---|---|---|---|
| B9 | −7.99 | −2.24 | −7.55 | −0.41 | 1.79 |
| A7 | −9.12 | −1.68 | −8.33 | −0.63 | 0.89 |
| B1 | −8.56 | −2.32 | −7.73 | −0.05 | 1.49 |
| P1 | −9.57 | −1.48 | −9.58 | −1.39 | 1.49 |
| Mesotrione Cryst. | −6.36 | −2.81 | −3.97 | −1.06 | 1.49 |
| Mesotrione – 69th | −6.56 | −2.05 | −6.01 | −2.15 | 1.49 |
Regarding the binding free energy calculated with AutoDock 4, the P1 compound was identified as the most stable ligand inside the enzymatic cavity. Although B9 exhibited the highest experimental activity within the calibration set, its binding free energy was less favorable than that of both B1 and A7.
To better understand ligand accommodation within the enzyme’s active site, key molecular interactions were analyzed. Metal chelation appears to be the primary factor underlying the ligand’s competitive inhibitory behavior—preventing substrate binding and, consequently, catalytic activity. In addition, π–π stacking interactions between the ligand’s aromatic moiety and phenylalanine residues (PHE321 and PHE356 in the 5YWG structure) support an effective binding mode consistent with inhibition [ref. 25, ref. 61]. Considering that the docking protocol treats the protein as a rigid system, any structural variations in ligand fitting mainly reflect differences in how individual inhibitor moieties interact with the surrounding amino acid residues. Figures 5 and 6 depict the three‐dimensional orientations of the ligand conformations and highlight the specific interactions established within the binding cavity.


From Figure 4, while B1 and B9 exhibited greater deviations in their aromatic moieties, the remaining compounds showed better structural alignment. Notably, analysis of Figure 6 indicates that all presented compounds were able to chelate the cobaltous ion. Regarding π–π stacking interactions, only B9 failed to establish an interaction with the PHE356 residue, although the expected interaction with PHE321 was preserved. The compound A7 displayed the most similar interactions to the crystallized ligand, maintaining both aromatic and LYS353 alkyl interactions. However, despite adopting a different binding mode, P1 exhibited the lowest binding energy among the compounds. As shown in Table 4, this is likely due to a strong Van der Waals contribution, reflected in the alkyl interactions (Pi‐alkyl: HIS248, PHE321; alkyl: LYS353, LEU308, and LEU359).

To assess the correlation between the QSAR‐predicted activities and the binding free energies of the 14 proposed compounds, a plot was constructed (Figure 6). The X‐axis represents the binding free energy obtained from the most favorable docking poses, while the Y‐axis corresponds to the three activity values derived from each periodic property used in the MIA‐QSAR procedure, as denoted in Table 3. For each sample, the error bars on Y were calculated as the standard deviation of the three measurements, and a linear regression was performed.
The linear fit yielded a slope of –0.25 ± 0.03 and an intercept of 5.6 ± 0.2. The correlation coefficients indicate a strong agreement between predicted activity and binding affinity (R2 = 0.884; adjusted R2 = 0.874). Compounds P1 and P2 occupy nearly identical positions on the plot, corresponding to the highest predicted activity and the lowest binding free energy. Additional analyses were performed to ensure the reliability of the results: outlier analysis gave studentized residual as 2.7084, and the lack‐of‐fit test yielded F = 1.177 with p = 0.345, indicating no significant deviation from linearity and confirming that the linear model adequately represents the system.
Molecular Dynamics
Molecular dynamics simulations were performed to evaluate the binding stability of each ligand within the enzyme cavity. Three ligands were selected for comparison: the reference ligand mesotrione, P1 (the top candidate from docking and QSAR analyses), and P2 (the second‐best candidate). For each ligand, the initial pose corresponded to the lowest‐energy conformation obtained from docking.
The analysis considered ligand RMSD convergence within the cavity, conservation of key residues via RMSF, and structural flexibility of both enzyme and ligand through B‐factors. As shown in Figure 7, all ligands converged within the first 20,000 ps, with fluctuations confined to the 1–3 Å range. No deviations exceeded 4 Å over the 100 000 ps simulation, although P2 showed slightly higher deviations between 80 000 ps and the end of the simulation.

From the generated trajectories, mesotrione and P1 exhibited the most similar RMSD behavior, with the P1 curve showing slightly smaller deviations until convergence at 80 000 ps.
RMSF analysis (Figure 8) identified five key enzyme residues: His198, His280, and Glu366, which coordinate the cobalt (II) ion, and Phe353 and Phe396 (previously Phe321 and Phe356 in the docking analysis, as no loop modeling was performed), which are involved in π–π stacking interactions with the ligand’s aromatic moiety. Due to the absence of coil regions in the crystallized 5YWG structure, some modeled residues displayed slight positional differences compared to the docking results. Across all monitored residues and ligand types, no deviations greater than 2 Å were observed. For the metal‐coordinating residues, the mean RMSF was 0.5 Å, whereas the phenylalanine residues exhibited greater flexibility over time. Specifically, Phe353 was more flexible in the presence of mesotrione (RMSF = 1.105 Å), while Phe396 showed increased flexibility with P2 (RMSF = 1.641 Å).

Although the small deviations observed for these residues during the dynamic analysis suggest structural conservation of the enzyme cavity in the presence of ligands, some regions of the enzyme exhibited higher deviations. To identify these flexible regions, B‐factors were obtained, and colored structural representations were generated using DiscoveryStudio. The B‐factor is an index used to describe the inherent flexibility of atoms in a crystal lattice; considering thermal vibrations, it reflects atomic displacements around equilibrium positions throughout the trajectory and is strongly correlated with RMSF values [ref. 52].
Given the similar findings across all molecular dynamics simulations, the macromolecular structures were superimposed. The coloring scale represents structural flexibility, with blue indicating the most conserved regions and red indicating the most flexible areas, as determined from B‐factors and atomic fluctuations around equilibrium positions (Figure 9).

In the superimposed structural representations, unfolded regions—such as loops and terminal segments—showed larger deviations and were highlighted as white/red residues, indicating a lack of convergence. In contrast, well‐defined secondary structure elements exhibited smaller deviations, reflecting their stability over time. These observations, consistent with previous HPPD studies, suggest that despite the presence of different ligands in the enzyme cavity, the overall effects of exogenous molecules on enzyme flexibility were similar [ref. 62].
Physicochemical Properties and ADMET Analysis
To assess the potential risks associated with the P1 molecule compared to the standard compound mesotrione, ADMET parameter predictions were performed. Predictions were conducted using the ADMETlab web server, chosen for its high accuracy, precision, and coverage [ref. 63]. The obtained results are summarized in Table 5.
TABLE 5: ADMET properties calculated for P1 and mesotrione.
| Parameter | P1 | Mesotrione |
|---|---|---|
| Solubility | −4.246 | −3.135 |
| Caco‐2 | −4.710 | −4.582 |
| BBB | 0.035 | 0.005 |
| PPB | 97.088% | 58.991% |
| Hepatoxicity | 0.596 | 0.583 |
| Skin sensibilization | 0.396 | 0.755 |
| Carcinogenicity | 0.581 | 0.579 |
Given the output values for each evaluated property, BBB penetration was assessed by the predicted probability of crossing the Blood‐Brain Barrier (range: 0–1), PPB by the Plasma Protein Binding index, and toxicity, skin sensitization, and carcinogenicity by the probability of being toxic or sensitizing (0 = false; 1 = true). Comparing P1 with mesotrione, it is evident that mesotrione exhibits higher solubility. Conversely, parameters such as Caco‐2 permeability, BBB penetration, hepatotoxicity, and carcinogenicity show similar performance for both compounds. Specifically, both molecules display low predicted intestinal absorption (with ‐5.15 log units as the lower threshold) and limited BBB permeability, along with a moderate probability of hepatic and cardiac toxicity. The most pronounced differences between the two compounds were observed for PPB and skin sensitization. P1 demonstrated a higher capacity for plasma protein binding, whereas mesotrione exhibited a greater likelihood of inducing allergenic or sensitizing effects.
Regarding lipophilicity, the tenth proposed compound displayed the lowest log P value (0.34), compared with 0.42 for mesotrione (Figure 10). In contrast, the most active compound, P1, presented a log P value of 3.23, indicating considerably higher lipophilicity.

Conclusion
The dataset used was suitable for applying the MIA‐QSAR method, providing strong evidence that P1 is the most promising herbicidal candidate among the studied and proposed compounds. Its potential was further supported by docking and molecular dynamics simulations.
In docking analyses, the high R2 value between predicted activities and binding energies indicates not only a strong correlation between these parameters but also the reliability of the original in vitro activity data, the close link between compound interactions and biological activity, and the effectiveness of the computational methods employed. During molecular dynamics simulations, P1 and mesotrione exhibited similar behavior and convergence patterns, suggesting that P1 could effectively inhibit HPPD, analogous to mesotrione. ADMET analysis further indicated comparable safety profiles, with P1 showing even greater predicted safety, although appropriate management practices remain necessary.
This study provides an initial framework for designing new selective and potent herbicides. The higher predicted activity of P1 compared to the commercial standard could translate into improved inhibition and potentially allow reduced application rates. Overall, these findings underscore the value of this computational approach for advancing environmentally sustainable and safe herbicidal design.
Funding
This work was supported by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior (CAPES, funding code 001), Conselho Nacional de Desenvolvimento Científico e Tecnológico (CNPq, grant number 306830/2021‐3), and Fundação de Amparo à Pesquisa do Estado de Minas Gerais (FAPEMIG).
Conflicts of Interest
The authors declare no conflicts of interest.
Supplementary Materials
References
- Synthesis and Herbicidal Activity of Triketone‐Aminopyridines as Potent p‐Hydroxyphenylpyruvate Dioxygenase Inhibitors,”. Journal of Agricultural and Food Chemistry, 2021. [DOI | PubMed]
- Discovery of Novel HPPD Inhibitors Based on a Combination Strategy of Pharmacophore, Consensus Docking and Molecular Dynamics,”. Journal of Molecular Liquids, 2022. [DOI]
- Design, Synthesis and Biological Activity of Novel Triketone‐Containing Quinoxaline as HPPD Inhibitor,”. Pest Management Science, 2022. [DOI | PubMed]
- Discovery of Novel Benzofuran Scaffold as 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitors,”. Pest Management Science, 2021. [DOI | PubMed]
- 4‐Hydroxyphenylpyruvate Dioxygenase (HPPD)‐inhibiting Herbicides: Past, Present, and Future,”. Weed Technology (, 2023. [DOI]
- Synthesis and Herbicidal Evaluation of Triketone‐Containing Quinazoline‐2,4‐diones,”. Journal of Agricultural and Food Chemistry, 2014. [DOI | PubMed]
- Advances in Understanding the Toxicity of 4‐Hydroxyphenylpyruvate Dioxygenase‐Inhibiting Herbicides,”. Journal of Agricultural and Food Chemistry, 2024. [DOI | PubMed]
- 4‐Hydroxyphenylpyruvate Dioxygenase,”. Archives of Biochemistry and Biophysics, 2005. [DOI | PubMed]
- Diversity and Origin of Carotenoid Biosynthesis: Its History of Coevolution Towards Plant Photosynthesis,”. New Phytologist, 2021. [DOI]
- Formation of α‐Tocopherol Hydroperoxide and α‐Tocopheroxyl Radical: Relevance for Photooxidative Stress in Arabidopsis,”. Scientific Reports, 2020. [DOI | PubMed]
- True Oxygen Reduction Capacity During Photosynthetic Electron Transfer in Thylakoids and Intact Leaves,”. Plant Physiology, 2022. [DOI | PubMed]
- Broad 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitor Herbicide Tolerance in Soybean with an Optimized Enzyme and Expression Cassette,”. Plant Physiology, 2014. [PubMed]
- Metabolism of the 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitor, Mesotrione, in Multiple‐Herbicide‐Resistant Palmer Amaranth (Amaranthus palmeri),”. Journal of Agricultural and Food Chemistry, 2024. [DOI | PubMed]
- WeedScience.org. (. 2025
- Study of Two Combined Series of Triketones With HPPD Inhibitory Activity by Molecular Modelling,”. SAR and QSAR in Environmental Research, 2023. [DOI | PubMed]
- Coloured Chemical Image‐based Models for the Prediction of Soil Sorption of Herbicides,”. RSC Advances, 2015. [DOI]
- MIA‐QSAR: A Simple 2D Image‐based Approach for Quantitative Structure–activity Relationship Analysis,”. Journal of Molecular Structure, 2005. [DOI]
- MIA‐plot: A Graphical Tool for Viewing Descriptor Contributions in MIA‐QSAR,”. RSC Advances, 2016. [DOI]
- Survey on the Recent Advances in 4‐Hydroxyphenylpyruvate Dioxygenase (HPPD) Inhibition by Diketone and Triketone Derivatives and Congeneric Compounds: Structural Analysis of HPPD/Inhibitor Complexes and Structure–Activity Relationship Considerations,”. Journal of Agricultural and Food Chemistry, 2022. [DOI | PubMed]
- 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitors: From Chemical Biology to Agrochemicals,”. Journal of Agricultural and Food Chemistry, 2017. [DOI | PubMed]
- Aqueous pKa Prediction for Tautomerizable Compounds Using Equilibrium Bond Lengths,”. Communications Chemistry, 2020. [DOI | PubMed]
- NMR and Computational Studies on Tautomerism of 3‐hydroxy‐2‐(2‐thienylcarbonyl)cyclohex‐2‐en‐1‐one,”. Journal of Molecular Structure, 2014. [DOI]
- Predominance of the Triketo Tautomer in Acyldipivaloylmethanes in Solution and the Solid state,”. Journal of Molecular Structure, 2014. [DOI]
- Structural Features of Two Pyridyl Compounds of 1,5‐Bis‐(2′‐pyridyl)pentane‐1,3,5‐trione and a New Salt of Doubly Protonated Hydroxyterpyridinium,”. Journal of Chemical Crystallography, 2020. [DOI]
- Design, Synthesis and Herbicidal Activity of Novel Quinazoline‐2,4‐Diones as 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitors,”. Pest Management Science, 2015. [DOI | PubMed]
- Benzofuran Derivatives and Their Anti‐tubercular, Anti‐bacterial Activities,”. European Journal of Medicinal Chemistry, 2019. [DOI | PubMed]
- Synthesis of Novel Five‐membered Nitrogen‐containing Heterocyclic Compounds From Derivatives of Arylsulfonyl‐ and Arylthioacetic and ‐Propionic Acids,”. Chemistry of Heterocyclic Compounds, 2001. [DOI]
- Discovery of Novel Arylthioacetic Acid Derivatives as 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitors,”. Pest Management Science, 2020. [DOI | PubMed]
- Synthesis and Herbicidal Activities of Aryloxyacetic Acid Derivatives as HPPD Inhibitors,”. Beilstein Journal of Organic Chemistry, 2020. [DOI | PubMed]
- MIA‐QSAR: A Simple 2D Image‐based Approach for Quantitative Structure–activity Relationship Analysis,”. Journal of Molecular Structure, 2005. [DOI]
- aug‐MIA‐QSAR Modeling of Antimicrobial Activities and Design of Multi‐target Anilide Derivatives,”. Journal of Microbiological Methods, 2013. [DOI | PubMed]
- The MIA‐QSAR Method for the Prediction of Bioactivities of Possible Acetylcholinesterase Inhibitors,”. Archiv Der Pharmazie, 2012. [DOI | PubMed]
- GaussView Version 6 (. 2019
- Chemoface: a Novel Free User‐Friendly Interface for Chemometrics,”. Journal of the Brazilian Chemical Society, 2012
- Prediction Reliability of QSAR Models: An Overview of Various Validation Tools,”. Archives of Toxicology (, 2022. [DOI | PubMed]
- 36 K. Héberger , A. Rácz , and D. Bajusz , “Which Performance Parameters Are Best Suited to Assess the Predictive Ability of Models?” in Advances in QSAR Modeling: Applications in Pharmaceutical, Chemical, Food, Agricultural and Environmental Sciences, K. Roy (Springer, 2017), 89–104, 10.1007/978-3-319-56850-8_3.
- Real External Predictivity of QSAR Models: How To Evaluate It? Comparison of Different Validation Criteria and Proposal of Using the Concordance Correlation Coefficient,”. Journal of Chemical Information and Modeling, 2011. [DOI | PubMed]
- Validation of QSAR Models–Strategies and Importance,”. International Journal of Drug Design and Discovery, 2011
- Further Exploring R Metrics for Validation of QSPR Models,”. Chemometrics and Intelligent Laboratory Systems, 2011. [DOI]
- Real External Predictivity of QSAR Models. Part 2. New Intercomparable Thresholds for Different Validation Criteria and the Need for Scatter Plot Inspection,”. Journal of Chemical Information and Modeling, 2012. [DOI | PubMed]
- Molecular Insights Into the Mechanism of 4‐Hydroxyphenylpyruvate Dioxygenase Inhibition: Enzyme Kinetics, X‐Ray Crystallography and Computational Simulations,”. FEBS Journal, 2019. [DOI | PubMed]
- 42 BIOVIA Discovery Studio , BIOVIA Discovery Studio, Release 2024 (Dassault Systèmes, 2024).
- AutoDock4 and AutoDockTools4: Automated Docking With Selective Receptor Flexibility,”. Journal of Computational Chemistry, 2009. [DOI | PubMed]
- 44 G. Bitencourt‐Ferreira and W. F. de Azevedo , “How Docking Programs Work,” Docking Screens for Drug Discovery, ed. W. Filgueira de Azevedo Jr. (Springer New York, 2019), 35–50, 10.1007/978-1-4939-9752-7_3.
- An Automated Strategy for Binding‐Pose Selection and Docking Assessment in Structure‐Based Drug Design,”. Journal of Chemical Information and Modeling, 2016. [DOI | PubMed]
- AmberTools,”. Journal of Chemical Information and Modeling, 2023. [DOI | PubMed]
- MCPB.Py: A Python Based Metal Center Parameter Builder,”. Journal of Chemical Information and Modeling, 2016. [DOI | PubMed]
- Parameterization of a Dioxygen Binding Metal Site Using the MCPB.py Program,”. Methods in Molecular Biology, 2021. [DOI | PubMed]
- Building and Evaluating Force Field Parameters Around Transition Metal‐including Active Site in 4‐Hydroxyphenylpyruvate Dioxygenase (HPPD),”. Journal of Computational Chemistry, 2022. [DOI]
- 50 M. J. Frisch , G. W. Trucks , H. B. Schlegel , et al., Gaussian 16 Revision C.01 (Gaussian Inc., 2016).
- GROMACS: High Performance Molecular Simulations Through Multi‐level Parallelism From Laptops to Supercomputers,”. SoftwareX, 2015. [DOI]
- B‐Factor Rescaling for Protein Crystal Structure Analyses,”. Crystals, 2024. [DOI]
- ADMETlab 3.0: An Updated Comprehensive Online ADMET Prediction Platform Enhanced With Broader Coverage, Improved Performance, API Functionality and Decision ort,”. Nucleic Acids Research, 2024. [DOI | PubMed]
- Identification of 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitors by Virtual Screening, Molecular Docking, Molecular Dynamic Simulation,”. Journal of the Science of Food and Agriculture, 2023. [DOI | PubMed]
- The Novel 4‐Hydroxyphenylpyruvate Dioxygenase Inhibitors In Vivo and In Silico Approach: 3D‐QSAR Analysis, Molecular Docking, Bioassay and Molecular Dynamics,”. Arabian Journal of Chemistry, 2022. [DOI]
- Integrated Virtual Screening and Validation Toward Potential HPPD Inhibition Herbicide,”. Journal of Agricultural and Food Chemistry, 2024. [DOI | PubMed]
- Design, Synthesis, and Bioactivity of Novel Ester‐Substituted Cyclohexenone Derivatives as Safeners,”. 2023. [DOI]
- A Comparison of Theoretical Methods of Calculation of Partition Coefficients for Selected Drugs,”. Acta Poloniae Pharmaceutica, 2006. [PubMed]
- LINCS: A Linear Constraint Solver for Molecular Simulations,”. Journal of Computational Chemistry, 1997. [DOI]
- A Molecular Dynamics Method for Simulations in the Canonical Ensemble,”. Molecular Physics, 1984. [DOI]
- Canonical Dynamics: Equilibrium Phase‐space Distributions,”. Physical Review A, 1985. [DOI]
- Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method,”. Journal of Applied Physics, 1981. [DOI]
- Canonical Sampling Through Velocity Rescaling,”. Journal of Chemical Physics, 2007. [DOI | PubMed]
