To Bind or Not to Bind? A Comprehensive Characterization of TIR1 and Auxins Using Consensus in Silico Approaches

1. Introduction

Agricultural production is fundamental from an economic and nutritional standpoint. It has been estimated that from 1961 to 2020, food production from agricultural sources has increased significantly, with current per capita values rounding almost 3200 kcal/day, at the end of the supply chain [1]. Yet, the annual State of Food Security and Nutrition in the World (SOFI) report published by the Food and Agriculture Organization (FAO) showed that more than 700 million people faced hunger in 2022 [2]. Thus, there is a pressing need to improve food security, which translates to the assurement of higher crop yields.

Some of the main problems that pose a significant challenge to safe, productive and sustainable agricultural practices are the herbicide resistance and toxicity. Firstly, there have been several documented cases of resistances of weeds to several chemical groups worldwide since the discovery of herbicide resistance in common groundsel (Senecio vulgaris) to triazine herbicides in 1968 [3,4]. Nowadays, herbicide resistance is globally widespread, affecting 96 crops across 72 countries. This phenomenon has been documented in 513 distinct instances, involving 267 weed species that have developed resilience against herbicides [5]. Secondly, despite the generally low mammalian toxicity of many newly developed herbicides, growing experimental evidence indicates that exposure to these chemicals can have detrimental effects on the development and/or reproduction of various mammalian species. Notably, 2,4-D (2,4-Dichlorophenoxyacetic acid) and its combination with 2,4,5-T (2,4,5-Trichlorophenoxyacetic acid) have been associated with reproductive issues and malformations in humans [6].

Auxins can be regarded as phytohormones, metabolically derived from L-Tryptophan, responsible for plant growth and development [7]. With Indole-3-acetic acid (IAA) being the most well-known example of such compounds, often deemed as the “master hormone” in plants [8]. Since the 1940s, synthetic derivatives from IAA have been developed and tested, with notable examples being naphthalene acids (1-naphthalene acetic acid); phenoxycarboxylic acids (MCPA, 2,4-D); benzoic acids (dicamba); pyridinecarboxylic acids (e.g., picloram), and quinolinecarboxylic acids (quinmerac, quinclorac); all of which showed strong responses in plants, inducing damage at high concentrations while eliciting higher stability when compared to IAA [9]. Auxin compounds reached peak commercial success as herbicides, just in 2000, synthetic auxins accounted for roughly 1.3 billion USD, approximately 9% of total sales that year [10].

The intricate auxin signaling pathway is orchestrated by at least three protein families: transport inhibitor response/auxin signaling F-box (TIR1/AFB) proteins, auxin/indole-3-acetic acid (Aux/IAA) transcriptional repressors, and auxin response factor (ARF) transcription factors (Yu 2013). At the heart of this signaling cascade lies the TIR1 which is part of a small family of F-box proteins comprising five additional AFB proteins (AFB1-AFB5). These six proteins collectively function as auxin receptors [11]. Auxin transport occurs within and between cells via PIN, ABCB, and AUX/LAX transporters. Then interacts with the SCFTIR1/AFB complex, which upon creation of the SCF-auxin-Aux/IAA complex, triggers the ubiquitination of Aux/IAA transcriptional repressors, leading to the activation of auxin response genes [12]. Therefore, auxin exerts its regulatory influence by directly binding to TIR1, promoting the interaction between SCFTIR1/AFB and Aux/IAAs. This interaction doesn't induce any conformational changes; instead, auxin acts as a molecular glue, facilitating the formation of a binding pocket for the Aux/IAA protein within the TIR1 structure [13].

Herein, we present an in silico profiling of TIR1 binding, coupled with a pipeline for de novo design based on machine learning. We hope this study serves as an initial step towards further optimization in synthetic auxin design, as these molecules can provide additional tools not just as herbicides, but as molecular guidelines for engineering of plant growth and development [14] .

2. Materials and Methods

2.1. De Novo Design of Putative Auxins and Molecular Decoys

The crystal structure of the TIR1/SCP complex of Arabidopsis thaliana was obtained from the Protein Data Bank. The PDBID 2P1P was selected as its values for resolution, R-free and R-work are very good, also it is co-crystallized with IAA. The protein was then preprocessed with PDBFixer to add hydrogens, IAA and the co-factor were removed prior to modelling.

The program OpenGrowth [15] was used for de novo generation of molecules within the binding site for IAA. As OpenGrowth uses Openbabel as backend for compound minimization and atom typing the version 2.4.1 was manually compiled. Parameters for OpenGrowth were implemented as described in Table S1 in the Supplementary Information. Minimization was carried out using the MMFF94 force field and the SMoG2016 scoring function was used for ligand ranking. In the end, more than 1200 ligands were generated with these settings. A heuristic threshold of -8.0 kcal/mol was applied to select candidates for virtual screening.

To establish a set of non-binders, molecular decoys were obtained with Deepcoy [16]; these were generated to match DUD-E descriptors in a 150:1 decoy-to-reference ratio following selection of top decoys using the provided DeepCoy scripts.

As proof of concept, two sets of molecular descriptors were obtained: Lipinski space, composed by molecular weight (MW), topological polar surface area (TPSA), cLogP, rotable bonds, number of H-bond acceptors, and donors. A second set was selected from the molecular descriptors available in Mordred [17]. Feature selection was first performed by discarding those with zero variance and with less than three different values (semi constant variables). Highly correlated descriptors were also excluded, i.e., those with a Kendall correlation coefficient higher than 0.9; the final set contained 395 of the original 1613 variables. Then a subset of features was selected using the Boruta algorithm [18], these are listed in Table S2*, in the Supporting information. Finally, these chemical spaces were visualized via principal component analysis (PCA) and t-SNE using DataWarrior [19].

2.2. Machine Learning

To evaluate the overall pertinence of the generated compounds, we trained different classifiers as proof of concept. These included logistic regression (LR), supporting vector machines (SVC), random forest (RF), and extreme gradient boosting classifiers (xgboost classifier; XGBC). The comparison was performed with Scikit-learn (v.1.2.2) [20] and XGboost [21] Python modules. Each model was trained and validated by ten-fold cross-validation, grid search optimization was used to tune their respective hyperparameters [22]. Each compound was represented as a concatenated vector of molecular descriptors selected with Boruta; then a second iteration was made adding topological information: extended connectivity fingerprints (ECFP4), WHIM and GETAWAY descriptors, calculated with RDKit (v.2022.9.5) [23]. Given the numerical features differences, each representation was standardized, subtracting the mean and scaling to unit variance using Scikit-learn.

Classification performance was evaluated with several measures, including: validation curves, learning curves, ROC curves, Matthews correlation coefficient and detection error trade-off curves.

2.3. Molecular Modelling of Auxins

To better understand molecular recognition by TIR1, we conducted a comprehensive characterization with structure-based methodologies, including mixed-solvent molecular dynamics, conventional molecular dynamics and enhanced sampling with coarse metadynamics. Detailed information is provided in the following subsections.

2.3.1. Mixed Solvent Molecular Dynamics

In order to gain additional information on the plasticity and overall presence of hotspots, within TIR1, molecular dynamics with mixed solvents was used. For this, the MXMD script [24] was modified to include the probes included in Scheme 1.

A modified version of the Ghanakota et al. protocol was used [25]; briefly, the protocol involves the creation of replica systems with a volume/volume concentration of 95/5% (water/probe). MXMD uses pre-equilibrated boxes of the selected cosolvent probe which is then incorporated in a cubic box of water buffering the whole system by 15 Å. These systems were then submitted to default relaxation protocol and 15 ns of equilibration time. Production times were set to 20 ns, as it provides a good compromise between sampling and computing time [26]. Each MXMD run consisted of 15 simulations; three for each probe, accounting for an overall sampling of 300 ns. MXMD simulations were carried out with PDBIDs 2P1M, 2P1P and 2P1N; which are TIR1 structure without ligands, just IHP cofactor, bound to IAA and bound to 2,4-D; respectively.

Hotspots were then evaluated by the following score:

$$Score=\sum _p^{probes}\sum _s^{spots}\sum _g^{grids}Z(p,s,g)$$

(1)

Where Zxyz is the occupancy of functional groups of the probes, converted to Z-scores using the following equation:

$$Z_{xyz}=\frac{X_{xyz}-<X>}{\sigma }$$

(2)

2.3.2. Assessment of Pocket Solvation and Its Role in Auxin Recognition

Molecular recognition is often grounded on the “desolvation paradigm”, which considers that water molecules need to be displaced by ligands during binding events [27]. Only recently, this notion has been challenged, with the discovery of protein families such as bromodomains, which possess a tightly bound water network, serving as enthalpy control for ligand binding [28]. Plus, it has been shown that robust description of binding modes is more consistent when ligands are properly solvated [29,30].

Thus, we conducted a series of simulations under different starting conditions, using the following PDBIs 2P1N, 2P1O, 2P1P and 2P1Q. Structures were prepared with pdbfixer in order to add missing atoms and/or residues, and then submitted to the protein preparation wizard in Maestro. Hydrogen bonds were optimized using PropKa, finally two replicas of complexes were obtained: one preserving water molecules in the vicinity of co-crystal ligands and one without them. These structures were then submitted to the System Builder utility to buffer them in orthorhombic boxes of TIP3P water and NaCl at 0.15 M.

TIR1-ligand complexes were then submitted to molecular dynamics simulations using Desmond (v.6.1) [31]. Systems were parametrized with the OPLS_2005 force field, simulation, integration was carried out with the RESPA algorithm using an inner timestep of 4 fs and 8 fs for far timestep. Hydrogen mass repartition scheme was implemented to allow the increase in timesteps [32]. The prepared system was relaxed using Desmond’s default protocol which includes Brownian dynamics and an equilibration period to attain NVT and NPT conditions. Target values were 300 K for temperature and 1 atm for pressure. Following system relaxation the complex was equilibrated for 50 ns under the NPT ensemble, using MTTK barostat and Nosé-Hoover chain thermostat. Production runs of 500 ns for both replicas were obtained and analyzed for protein-ligand contacts, RMSD, RMSF and interaction transition networks.

Finally, we tested pocket solvation using grand-canonical Monte Carlo (GCMC), as it is a well-established method to achieve reasonable results [33,34,35,36]. As proof of concept, we selected the IAA/Tryptophan pair, each complex prepared as described above. Solvation protocol made use of the GCMC plugin in Desmond, defining the MC region in a buffer of 6 Å around ligands; solvent parameters were default values for TIP3P water, similar to those described elsewhere [37].

Simulations were carried out in five stages, beginning with GCMC/MD with restraints on heavy atoms. Consisting in 10 000 MC moves per cycle and a total of 500 ps of MD. This stage was followed by GCMC/MD with 100,000 moves per cycle, 1 ns of MD and no restraints. Systems were then relaxed with small runs of Brownian dynamics and Langevin dynamics for equilibration of volume and pressure; integration timestep was kept at 2fs and 4fs using the RESPA integrator. An additional equilibration of 10 ns using MTTK barostat and Nosé-Hoover thermostat was carried out, this step also used hydrogen mass repartition to increase timestep to 4fs and 8fs. Finally, production runs of 500 ns were made for each ligand keeping conditions from the previous stage.

2.3.3. Molecular Dynamics

We devised a series of comparisons between known ligands and chemically related “decoys''. Selected compounds included: IAA/Tryptophan; α-naphthylacetic acid/β-naphthylacetic acid; Diclophenac/Rinskor and Fluroxypyr/Quinclorac. By decoy, we intended to have a low affinity/non-binder for TIR1, based on experimental data [38,39]. As a first step, compounds were docked to PDBID 2P1P using PLANTS software, as it is a well-established method with high accuracy and consistency in the identification of putative binding modes [40,41,42,43]. Ligands were prepared using SPORES [44], to obtain proper input for PLANTS. The binding site was defined in a sphere with a radius of 15 Å, centered with IAA’s crystallographic coordinates. Additional parameters were kept at default values.

The protocol was the same as the generated compounds. The top ranking pose was selected for molecular dynamics simulation. Trajectory analysis was done under MDTraj (v.1.98) [45]. TIR1-auxin contacts were used to define the microstates to construct the networks. Interaction fractions were obtained with Maestro utilities and concurrence maps were calculated using ContactMap explorer (v.0.7.0) [46].

2.3.4. Coarse Metadynamics

To evaluate the metastability of binders and decoys, well-tempered metadynamics was used. This is a common method used to assess binding mode stability and overall affinity [47,48,49,50,51]. Simulation parameters were kept the same as from classical simulations. Collective variables (CVs) are selected degrees of freedom which are related to the desired reaction coordinate; in this case CVs were defined using distances between centers of mass, for the protein residues 435-437 and ligand’s heavy atoms; second CV was WHIM’s second principal moment for ligand’s heavy atoms. Gaussian height was set to 0.08 kcal/mol, with a width of 0.25 Å. The bias factor was set to 5.0 and Gaussians were deposited every 250 steps.

Additionally, we probed the role of the so-called engagement niche, a region composed by residues K410, S440, G441, A464, P465 [52]. For this case, compounds were docked at a centroid of these residues in PDBID 2P1M, CVs were kept the same. Both metadynamics simulations had production times of 50 ns.

3. Results & Discussion

3.1. De novo Design of Putative Auxins and Decoys

Using principal component analysis, we made a projection of the chemical spaces described in the methods section, loading values can be found in Table S3, in the supporting information. Figure 1 shows the chemical space delimited by the selected descriptors with Boruta. It can be seen that there is good superposition between generated binders and decoys in both principal component analysis (PCA) and t-distributed stochastic neighbor embedding (t-SNE). This suggests that classes may not be linearly separable, a similar behaviour is seen for the Lipinksi space (Figure S1). Such a trait is positive, as it shows that decoys are indeed robust due to sharing similar physicochemical spaces, which in turn may lead to effective discrimination of its latent space with a good value as screening tools.

Plus, this lack of clear separation underscores the challenge of distinguishing these classes. As it has been shown that overreliance on molecular properties can lead to negative effects such as molecular inflation [53,54].

With this, we concluded that a good baseline was established and proceeded to test the pertinence of these compounds using machine learning algorithms.

3.2. Machine Learning

Classification and prediction of biological activity based on latent spaces has been attempted since early quantitative structure activity relationships (QSAR). Still, after several decades and unexpected results it became clear that QSAR are hardly linear in nature [55] with the added complexity of highly similar compounds with highly variable biological activity; an issue now defined as activity cliffs [56] or even activity ridges, when said concept is found consistently in certain scaffolds and not just compounds [57].

With the rise of the Rule of Five [58], the existence of latent spaces for drug discovery gained traction once more. As of today, the notion of a druglikeness profile has been challenged [59]. Herein, we do not intend to find a set of empirical rules to characterize putative auxins, but rather a good starting point for chemical space exploration and expansion. We aim to provide a setting stone of sorts, to aid rational design of auxins and related agents.

To this end, we used the generated compounds to train commonly reported classifiers, in order to verify chemical space coverage and the overall significance of structural information within our data set. For classification problems several methods are used in conjunction to evaluate the accuracy and consistency of predictions. Hereunder we selected a combination of four measurements that convey a general picture of both data and model’s performance:

a) Validations curves: As these help to visualize the trend of variance and suggest the presence of bias; either underfitting or overfitting.

b) Learning curves: Which convey information on the learning rate as a function of n in data; i.e., if the provided examples are enough for effective classification or if the addition of data can benefit scoring.

c) ROC curve & Matthews correlation coefficient (MCC): These provide a general description of the predictive power and can be linked to the confusion matrix. The inclusion of both is complementary, as MCC can take values between -1 and 1, it is possible to assess if a model is indeed providing a significant difference when applied or if its success relies more on chance. Even when the area under the curve of the ROC plot would suggest otherwise.

d) Detection error trade-off: Similar to the ROC curve, this metric gives information on the success of classification tasks. During classification, there is always a trade-off between the rate of falsely classified values; therefore this plot conveys the trend of said trade-off; an ideal classifier would keep a low rate of either false positives and negatives, resulting in a curve downward and to the left, an inverse of sorts to an expected ROC curve.

The plots for the first three metrics for all classifiers can be found in the supplementary information (Figures S2-S9). ROC curve and Error trade-off for the models with the best parameters are shown in Figure 2.

As expected, models trained solely with physicochemical data showed low performance. Based on area under ROC curves, algorithms based on linear methods performed worst. Plus, it is quite significant that classifiers were barely superior to random chance. When analyzing the detection error trade-off, it also gives an almost linear trend. This suggests that any improvement in the discrimination of either false positives or negatives inherently results in an increase of misclassification of the other. This underscores that proper discerning between auxins and non-binders is no easy task, even random forest and supporting vector classifier had poor performance, which are often regarded as robust, for chemical information applications due to their non-linear approach [60,61,62,63,64]. In addition, these classifiers showed a strong tendency to overfitting (Figures S4 and S6 in the supporting information).

Even if these results seem negative, such is not the case, as this verifies that molecular decoys are robust indeed. A common use case of molecular decoys is for structure-based design; i.e., decoys should be molecules that share chemical space with actives, with little to none pharmacophoric similarity. The quest for such compounds resulted in the construction of datasets such as DEKOIS [65] or most notably DUDE-Z [66]. Only recently, use cases have been extended to ligand-based methods and scoring function construction [67,68]. Hence, for our case at the very least, we have confirmed that decoys do share the chemical space from selected descriptors, but also a very representative one for small molecules (Ro5 space or Lipinski space).

To confirm the second half of the provided definition, topological and structural features shall be more conclusive in this classification test. Indeed addition of these proved very significant, as the area under ROC curves increased roughly by 0.2 in all cases. Random forest classifier had the lowest AUC value, consistent with its MCC and error trade-off. Following this logic, it would seem that logistic regression was the top performer; yet, its MCC was 0.53. For this test non-linear methods were more consistent and performed well, at first glance SVC would seem the superior choice; then again when comparing MCC values suggests that XGBC would be preferred.

Generally speaking, our goal was achieved, as the generated compounds do serve as a reasonable sample for further analysis,while conveying enough information on true binders. Still, there is no denying that the rate of false negatives remains quite high even at notable trade-off values and thus major incidence of false positives. To some extent, this result was to be expected, as no computational method is completely accurate and reliable on its own. From here, we analyze the results of structure-based methods to provide further information on what characterizes a true binder of TIR1.

3.3. Molecular Modelling of Auxins

3.3.1. Mixed Solvents Molecular Dynamics

Structural alignment suggests little variation between conformations of TIR1; nonetheless, in a dynamic context substantial differences may arise [69,70]. Hence, to attain a comprehensive description and knowledge on these factors, we conducted molecular dynamics with mixed solvents. During drug discovery campaigns, molecular design can follow a fragment-based approach. A coarse description of the technique is that protein structures are crystallized iteratively, in the presence of probe scaffolds or leads. Thus, creating an empirical mapping on which fragments bind more effectively and where does this occur [71].

A computational equivalent of this can be found in molecular dynamics with mixed solvents, where a rather discrete amount of polar solvent molecules is introduced to the simulation cell. Following an equilibration period, the probe molecules will interact with the protein, and therefore occupancy of these can be evaluated and further compared. A notable advantage of the computational variation is that calculation of occupancy can be made from different probes from independent simulations. Which in turn increases descriptive power and can be used to delimit allosteric or even cryptic pockets within proteins [72].

We conducted 45 independent simulations accounting for 900 ns; using three main conformations of the auxin binding pocket. Figure 3 shows the computed surfaces from this sampling. It can be seen that the auxin binding pocket had the highest scoring, indeed. However, significant changes in the neighboring region to the engagement niche. Additional changes included the orientation of F79, 82 and 380.

In general, the conformational variations in TIR1 pocket changed the shape of the druggable regions. Overall, the ring shape is suggestive of the known mechanism of auxins, as these stabilize the protein-protein interface. Considering that said shape is notable in 2P1N, its worth considering the link between ligand polarity and pocket dynamics. As for individual probes, for 2P1M the probes with major occupancy were N-methylacetamide and acetohydroxamic acid. In contrast, 2P1P and 2P1N showed a major occupancy of druglike-probes, particularly imidazole.

3.3.2. Assessment of Pocket Solvation and Its Role in Auxin Recognition

As stated in the methods section, only recently the role of pocket solvation has been critically assessed within the context of small molecule binding. Based on crystallographic structures it does seem that TIR1 has several water molecules that could participate in auxin binding [73]. To determine the contribution to recognition and binding we compared protein-ligand interaction under different starting conditions of solvation (Figure 4).

Notable differences arise when comparing interaction profiles, major changes are the contact with R436, which is mostly driven by water bridges; suggesting a change in ligand orientation and directionality. Similarly, it seems that water interactions are significant within the engagement niche. H78 also shows subtle differences, it is difficult to assess if these are the result of solvation or merely the nature of numerical simulation. Another subtle effect is the change of hydrophobic contacts to stacking interactions, it can be argued that due to the lipophilicity of auxins these are prone to the latter. However, considering the cases of NAA and 2,4-D (2P1O and 2P1N, respectively) interaction distribution and persistence did change significantly due to the solvation environment. Figures S11 and S12 show some examples of this (vide infra).

For further proof, we constructed transition networks based on protein-ligand interactions. Briefly, the networks encode microstate information which can be related to metastability and with enough data it can confer kinetic information [74]. Thus, it is very important to delimit which information is used to describe said microstates; hence, an additional criterion was introduced, networks were constructed with and without the inclusion of water bridges in their definition. Figure 5 shows the obtained transition networks, overall it can be seen that just as previously suggested, the inclusion of water-ligand interactions has a significant impact. metastability was present in all three simulations with crystal waters; yet, this trait is further accentuated by the introduction of water bridges to microstate definition. Such behaviour is more noticeable in the synthetic auxins, which provides further support to the previous argument on the role of ligand polarity.

In stark contrast, networks without preservation of crystal water were highly variable. For the case of IAA (2P1Q) are reversal pattern is observed, as it shows some degree of metastability and a densely populated network. NAA on the other hand showed highly transient interactions, some states are indeed revisited but no apparent clustering arises under these conditions. Omission of water bridges does result in a more complex network but the transient nature of contacts persists. Finally, for 2,4-D metastability is clearly shown, suggesting a local minima or basin. Overall the inclusion of water mediated interactions did improve the representation and gave a more robust description of basins and metastability.

It can be argued that such observations are to be expected due to the complex nature of binding and the balance of enthalpic and entropic contributions. Nonetheless, our aim is not to strictly evaluate the role of water interactions, but their general contribution to discern between binders and non-binders. To this end, we conducted GCMC/MD simulations with IAA and TRP. For these simulations PDBID 2P1P was used, comparison was between a common MD protocol and GCMC/MD, using similar analyses described above. Figure 6 shows protein-ligand interactions retrieved from IAA and TRP. It can be seen that for the most part the GCMC/MD protocol was successful in establishing a water network were the main bridge between ligand and R436 is present. This is a very positive result as it highlights the robustness of the method and its placement of solvent molecules.

On the other hand, the interaction profile remains consistent and similar to those in the previous simulations. Providing further proof on the significance of water bridges in microstate definition. For IAA perhaps a notable change is the dominant stacking interaction with F79, instead of a hydrophobic contact. Due to this, we made a structural alignment and binding site superposition of TIR1 crystal structures; i.e., PDBIDs 2P1M, 2P1P, 2P1Q, 2P1O and 2P1N. Results of this alignment are shown on Table 1.

Based on these alignments, F79 and 82 are consistently flexible between structures. A notably higher deviation is shown precisely between 2P1P and 2P1Q; which possibly accounts for this transition between hydrophobic and stacking interactions. This together with the observed changes in molecular dynamics with mixed solvents, reinforce the importance of stacking interactions within the binding pocket.

Other residues showing changes between TIR1 conformations are H78, L439 and A464. With the latter being more displaced when 2,4-D is bound to TIR1. This trait is noteworthy as A464 belongs to the engagement niche. Perhaps the observed variation for L439 can also be linked to the presence of water bridges. Unfortunately, this was only observed for IAA so a clear reason remains elusive. Nevertheless, this may be due to the presence of nitrogen in the indole scaffold, a feature absent in other binders under study. Still previous studies have suggested that the presence of this heteroatom is paramount, for potency and overall affinity of IAA [75].

In contrast, one of the most consistent residues was R436; which in turn suggests that the water bridge between it and auxins can be regarded as “canonical” and not anecdotal behaviour. Plus, considering the complete absence of said bridge in conventional protocols is very interesting. A similar case was obtained in 2P1O where the bridge was not formed under standard conditions. Of note, simulations of both 2P1Q and 2P1N were able to form this interaction. The fact these subtle variations within the binding site of TIR1 introduce such divergent results is remarkable. Also, this stands as a cautionary tale for the well known limitations of structure base methods, while highlighting that structure selection is no trivial matter [76,77]. Thus, it seems that 2P1P is prone to the loss of ligand-solvent interactions; whereas 2P1Q and 2P1N may yield more reliable and consistent results. For further assessment, we analyzed root mean square fluctuation of auxins in these sets of simulations (Figs S10-S13).

Considering IAA and TRP in 2P1P, RMSF plots indicate higher fluctuation when water molecules are correctly placed in the pocket. Under standard protocol, the fluctuations are significantly suppressed, more so for IAA’s case. For both ligands, higher fluctuation was observed in the nitrogen’s side of the scaffold. This is worth mentioning as this part of the molecule faces in the vicinity of F79. When comparing to NAA and 2,4-D interesting patterns were found, for the former the fluctuation profile was very consistent, with respect to its fit in the binding site. The main difference being found on the carboxylic acid moiety, which may be explained by the loss of the water bridge with R436. For 2,4-D, fluctuation profile changes but in contrast to IAA, TRP and NAA, higher values were observed under standard conditions. Due to this, and considering that atoms showing this behaviour were both chlorines and the carbon atoms of this side of the phenyl ring; this is yet another case where the solvation environment is the main driving force.

As for ligand’s fit, fluctuation was very consistent and almost unaffected. In all cases the highest value was found in the charged oxygen, for TRP this is to be expected due to the amino group. Beyond that, only 2,4-D showed differences in fluctuation for this atom between protocols.

3.3.3. Molecular Dynamics

Considering prior results, we used an integrative approach for the analyses of known binders and non-binders. In an attempt to facilitate the correct assessment of binding; we propose the use of protein-ligand interaction fractions; concurrence maps, transition networks and ligand-protein interactions simultaneously. Perhaps it would seem that this combination conveys redundant information. Still, let us present two cases as proof: diclofenac (DCF) and quinclorac (QCL). Both are arylic carboxilic acids, with quinclorac being an auxin herbicide; however the molecular target of quinclorac is not TIR1, but rather F-box protein ABP5 [38].

DCF is an non-steroidal anti inflammatory agent (NSAID) and thus, not an auxin per se; but it has been shown that exposure to this drug can cause effects similar to those of IAA [39]. Figure 7 present the consensus analyses of molecular dynamics simulation of DCF after docking in TIR1. Indeed, it can be seen that the interaction profile is similar to that of synthetic auxins; however, there is no ionic nor stacking interactions with H78. This may be attributed to the ring’s orientation due to hindrance of the chlorine atoms. Another significant observation is the diminished presence of the water bridge with R436.

When analyzing contact concurrence and interaction network, these suggest that DCF remains in the general vicinity of the main pocket. This, plus the observed basins and metastability can be taken as an indication of true binding. Nevertheless, the diminished contact with the water bridge was also observed in the ligand-protein contacts.

These observations provide a strong case in favor of the proposed protocol, still a criterion related to binding affinity is needed to assess the relative potency of a given ligand.

On the other hand QCL has no experimental affinity for TIR1. Figure 8 shows the results of the consensus analyses from MD. Based on the interaction profile, all contacts are significantly weaker for quinclorac. Moreover, the water bridge with R436 is non-existent; perhaps it is due to this that the interaction network shows a highly transient behaviour. No basins were observed and this is consistent with the concurrence plot. It suggests that quinclorac remains near hydrophobic residues; i.e., F82 and 380. Due to its planarity the presence of stacking interactions with these resides may be enough to stabilize QCL in the pocket. Beyond that, any significant contact is barely present; the main anchor seemed to be R489, as confirmed by ligand-protein interactions. Yet, this interaction seems to be lost in favor of hydrogen bonding with H78.

Additional analyses for other auxins can be found in the supplementary material, Figures S14-S19. Akin to simulations in section 3.3.2, IAA showed a highly stable binding. Stacking interactions with F82 and 380 were observed. Transition network showed only one highly populated basin, with significant metastability. Of note, is the persistence of the water bridge with R436. Interestingly, as shown in the previous section TRP also showed a highly populated basin; still, the interaction profile was inconsistent with that of auxins, just as QCL. Moreover, it must be stated that based on concurrences and ligand protein interactions, TRP showed a very stable binding mode. Even after 500 ns this non-binder was able to remain in the pocket. This may be attributed to the main anchors being present for this ligand, with a strong presence of thr R436 bridge. This also goes to show that even within “long” simulations non-binders can be falsely identified as candidates due to pharmacophore similarity. Just like DCF this is yet another case where affinity would provide a decisive criterion for assessment and pruning.

As for NAA, both isomers had similar concurrence maps and ligand-protein interactions. Despite this, the interaction profile of the α-napthyl acetic acid showed stacking interactions with H78, F79 & 82; traits absent for the β isomer. The main difference however, was found in transition networks. With the one of NAA showing a clear basin with significant metastability. Whereas 2-NAA showed a network more in line to the one observed for QCL. Experimental evidence confirms that indeed, NAA has the higher potency.

Finally, we discuss the results for fluroxypyr (FXY) and rinskor (RSK); both auxin herbicides belonging to the picolinate subclass [78]. Similar to QCL, these show preference to ABF proteins; still fluroxypyr is able to bind significantly to TIR1. Interestingly, FXY seems to be a binder whose main anchors are H78 and R489. Moreover, this arginine exhibits an analogue behaviour to that of R436. This, added to its general positioning, makes us hypothesize that auxin orientation initiates with K410 just outside the engagement niche; with R489 serving as the pulling mechanism toward the auxin site. A premise supported by the transition network, where some minor basins seemed to be found with shared microstates. Based on the interaction profile and the presence of additional water bridges, this can be the origin of said transitions. Concurrence plot also seems to favor this, since FXY stays in the general vicinity of H78 during the simulation; making close contact with D487 only during the last segment of the run.

As for RSK, the interaction profile was dominated by hydrophobic and stacking interactions. Notably, water bridges were also present, again with R489. A crucial detail however, is that RSK only showed inverse binding; i.e, a general orientation toward the engagement niche and R489 and the phenyl rings in the vicinity of H78. Uzunova et al., found similar arrangements for TRP using tomographic docking [52]. Perhaps this can be attributed to the absence of the carboxylic acid moiety, esterified in RSK. Of note, both RSK and related picolinate Arylex are readily metabolized to carboxylates, it has been suggested that their high potency is due to this activation step [79]. In our case, the preservation of the ester was deliberate; precisely to test the interaction profile and prior observations. Transition network shows a putative basin with transient behaviour; still as the simulation continues the microstates are novel and the basin no longer explored. In this regard the concurrence map suggests that overall RSK remains very stable during most of the run. Therefore, this is yet another case where water bridges may be the reason behind this pattern.

3.3.4. Coarse Metadynamics

Finally, we conducted a series of coarse metadynamics simulations in order to have an additional criterion for cases where consensus analyses could prove ambiguous. Briefly, metadynamics belongs to enhanced sampling methods, for this case a bias is introduced as specific degrees of freedom (CVs) which are subject to Gaussian potentials. These potentials help to fill the minima defined by CVs; in the context of small molecule binding a coarse approach has been proposed elsewhere [48]. While it is true that metadynamics has been successfully applied to determine binding mode stability, a major challenge is CV selection.

Herein, we used distance between centers of mass as it has proven to be a very robust choice in previous studies. Now, for ligand’s degrees of freedom, we settled on the WHIM descriptor. Due to its conceptualization and the overall positive results found during machine learning. Hence, WHIM’s second principal moment of the ligands was used. Commonly, metadynamics simulations are made from docking poses; thus, in order to test the engagement niche we carried out docking in the vicinity of this region and used metadynamics to “force” the ligands towards the binding pocket. Figure 9 shows free energy surfaces (FES) for DCF.

This compound only showed one basin in each metadynamics run. Distance and affinity values for these basins were close. Nonetheless, when binding poses within the basins were selected an interesting picture emerged (Figure 10). When DCF starts within the pocket and is pulled outward; it gets stuck in the vicinity of F351 and R489. In contrast, when DCF starts from the engagement niche, it is able to enter the auxin pocket and adopts a canonical orientation toward H78, F82 and 380; an indicative of true binding. It is true that in both arrangements the DCF is enclosed by several water molecules; still this could be due to the noisy behaviour of short metadynamics. Now, taking all prior results and the computed affinity value (~ -19 kcal/mol) confirms the role of DCF as a true binder. Nevertheless, the significance of this value remains relative to that of other auxins calculated with the same protocol, discussed below.

Recovered FES for QCL are shown in Figure 11, for this compound, a clear unbinding event took place. Additionally, this was the only basin found, when visual inspection of the binding mode was done (Figure 12) QCL is outside the neighboring region of the engagement niche. Contrary to DCF, the synthetic auxin was not able to effectively bind to any region during its exit. When QCL starts from the engagement niche, based on COM distance values it would seem that binding did occur. Again, when the basin is inspected it would seem that a true binding was indeed achieved, as QCL is interacting with H78; however, this is only via hydrogen bonding, there is no interaction with F79, 82 or 380 nor any contact or presence of any arginine residue. Even if the suggested affinity seems rather high (-32 kcal/mol) this binding mode is not characteristic of any auxin under study. Also, this value can be the result of low sampling of the observed basin. As stated, even affinity estimates from metadynamics are qualitative in nature, more so when short runs are produced, given enough time this value is likely to diminish towards a more representative one.

Now for the additional auxins and non-binders; their FES diagrams and representative binding poses from basins can be found in the supplementary information (Figures S20-S31). Herein we just provide some details and highlights, but we think that results are self-explanatory. IAA, serves to establish a baseline for affinity values. Akin to experimental reports, it showed the highest value of all true binders (-24 kcal/mol) when initiating from within the binding pocket. Followed closely by NAA isomers (-22 kcal/mol), also consistent with experimental affinity which is usually one order of magnitude below IAA [14,38,80]. For these cases metadynamics could not distinguish between isomers just with affinity estimates; nonetheless, only NAA showed true binding modes in both metadynamics runs.

TRP was correctly classified as a non-binder, even when a quasi-canonical arrangement was found; the affinity estimate was very low (-16 kcal/mol). RSK showed even lower values (-14 kcal/mol), notably this was consistent between runs, and also its tendency for remaining inversely placed within the pocket.

The case for FXY is a rather interesting one, as an unbinding event occurred and it remained tightly bound to the engagement niche surface. Given the previous results we cannot provide an explanation for this discrepancy. Indeed, this case presses for further study while also highlighting the need for additional experimental information.

4. Conclusions

Herbicide development is an endeavor of high interest, from both an economical and chemical standpoint. Currently, there are worldwide efforts to replace commonly used agents such as glyphosate or 2,4-D. Herein, we presented a computational characterization of TIR1, one of the molecular targets associated with auxin pathways. With the aid of structure based methods we analyzed a series of compounds in silico, we found that water plays an important role in binding and were able to characterize significant protein-ligand interactions.

Using de novo design, we proposed a library of putative binders and molecular decoys. Chemical space showed a good superposition between them, a trait further supported by machine learning results. Herein, non-linear methods were more consistent and performed reasonably, while conveying enough information on true binders.However,current results also showed that current data makes correct classification of false negatives difficult. This behaviour was recurrent in this and other studies.

Hence, we conducted structure-based characterization of the TIR1 pocket. In order to determine major traits that could provide a robust protocol for true binder assessment. We found that inclusion of water mediated interactions gave a more robust description of local minima and metastability. Also, we found a persistent water bridge between ligands and R436; which to the best of our knowledge has not been described as significant.

We also probed the engagement niche, confirming its role in molecular recognition and the overall plasticity of the auxin binding pocket. Changes in ligand polarity make a significant impact on pocket dynamics and the orientation of hydrophobic residues involved in stacking interactions.

Using coarse metadynamics, we were able to further assess the affinity of known auxins and TIR1. Overall, the collective variables used in this study provide a robust description of ligand degrees of freedom and provide yet another criterion for virtual screening campaigns.

We are confident that the database and proposed analyses in this work provide a valuable starting point for further developments such as molecular generators or virtual screening campaigns.