Computational study of the interactions among structural analogues of acyl homoserine lactones (AHLs) and the Agrobacterium tumefaciens TraR binding site

Dataset

The chemical structures of 34 AHL analogues were taken together with their half maximal inhibitory concentration (IC₅₀) biological activity from Geske et al. (Table 1)^11,12. These structures were optimized using the theory of functional density with the functional X3LYP, together with the DGDZVP calculation base^13,14, using the Gaussian09 program¹⁵ (as an open-source alternative, GAMESS could be used¹⁶).

Table 1. Structures and IC₅₀ values of AHL analogues.

Group	2D Structure	ID	Substituent (X)	IC50 (µM)
AHL		4	-	0.69
		A3	-	0.77
		A4	-	1.05
Other analogues		5	-	0.83
		6	-	1.12
		A9 n=1	-	0.61
		A11 n=3	-	0.83
		A12 n=4	-	3.49
		D6	-	0.57
		D17	-	1.40
		E12		4.2
PHL		8	4-Br	4.73
		C10	4-I	1.25
		C13	3-NO3	2.25
		C20	4-CF3	0.81
		E6	3-SCH3	8.4
		E8	3-SCF3	4.7
		E11	3,4-Cl	2.4
POHL		D15	4-OCF3	0.46
		E15	4-CH3	1.8
		E16	3-CH3	1.1
		E17	4-OCH3	4.3
		E19	4-F	2.7
		E20	4-Cl	0.62
		E21	4-Br	0.51
		E22	4-I	0.44
		E23	4-CF3	1.6
		E24	4-NO2	0.29
PPHL		B7	4-Br	0.92
		E30	4-Cl	3.3
		E32	4-I	2.1
		E34	4-SCH3	1.1
		E35	4-CF3	0.21
		E37	4-NO2	1.6

AHL, acyl homoserine lactone; PHL, phenylacetyl homoserine lactone; POHL, phenoxyacetyl homoserine lactone; PPHL, phenylpropionyl homoserine lactone. Structures and IC₅₀ values for each analogue were obtained from Geske et al. 2007¹¹.

Molecular docking

Molecular docking calculations were performed to predict the binding modes among the 34 analogues of AHL and the TraR receptor. The 3D structure of the TraR receptor and its native co-crystallized ligand 3OC8HSL were downloaded from the Protein Data Bank, with the identifier 1L3L. The complex in PDB format was visualized and prepared using the structure preparation tool available in the Sybyl X 2.1.1 package¹⁷ and the Amber force field was used to assign the partial atomic charges of the protein. As an alternative, the open-source program, Avogadro could be used¹⁸. The native ligand and all water molecules were removed from the complex, hydrogen atoms were added and side chain amides and imidazoles were fixed (protonated), assuming a physiological pH of 7. The structure of extracted native ligand was optimized with the same conditions as the 34 analogues, mentioned above.

Once the structure of the native ligand and the TraR receptor was obtained, the docking protocol was initiated in the molecular simulation program AutoDock 4.2¹⁹. The preparation of both the protein and ligand was performed using the graphical user interface of AutoDock, called AutoDockTools²⁰. This preparation involved the addition of hydrogen and the assignment of Gasteiger charges²¹. After this, the potential maps for each type of atom in the ligand were calculated, executing the grid parameter using the Autogrid program incorporated in AutoDock. In this study, the grid box was delimited, centered on the TraR binding site, with dimensions 50x50x50 Å and spacing between points of 0.375 Å. Finally, for the conformational search, the Lamarckian genetic algorithm was chosen and for the validation 500 conformations were registered, with the remaining parameters used as per the default setting. With these parameters, docking validation was carried out by re-coupling the native ligand to the binding site of the TraR receptor.

Once the docking protocol was validated, and its ability to reproduce the interactions that promote protein-ligand binding confirmed for this particular case, TraR-3OC₈HSL, the docking calculations for the 34 AHL analogues and the receptor in question were continued following the previous protocol, with the only difference being that 200 conformations were recorded for each analogue, of which the best were chosen according to their binding energy values and their position within the binding site of the TraR protein. The interactions of the selected protein-ligand complexes were analyzed on the Protein-Ligand Interaction Profiler webserver²².

Construction of 3D-QSAR models

Two 3D-QSAR models were constructed from structural alignments. The first model was constructed by taking the structure of the ligand, with the best inhibitory activity, reported by Geske et al. 2007¹¹, as a template. The other model was constructed based on the 34 best-selected poses of each analogue-TraR molecular docking²³, which generated an automatic alignment²⁴. The comparative molecular field analysis (CoMFA)²⁵ descriptor was obtained by using the QSAR tool implemented in Sybyl X 2.1.1. Steric and electrostatic field energies were calculated by using the Lennard-Jones and Coulomb potentials, respectively, with a dielectric constant 1/r based on the distance of all interactions. To calculate the steric and electrostatic energies, we used the force field of Sybyl X.2.1.1, Tripos²⁶. As an alternative, the open-source program, Open3DQSAR could be used²⁷.

For the evaluation of the 3D-QSAR models obtained, the cross-validation technique was used with the ‘leave-one-out’ method. With this validation, the predictive power of the models was determined and, in this way, the best model was determined. An exploratory analysis of the correlation between molecular fields and biological activity was also performed using PLS regression; thus, the correlation coefficients (R²) of each model were obtained using the SAMPLS algorithm of Sybyl X 2.1.1. Finally, the cross-validation was done, giving the values of Q² and R²_CV. According to these values both models were compared, and it was deduced which had the best predictability.

Docking validation

In order to guarantee the reliability of the protocol used for molecular docking among analogues of AHL and the TraR receptor, validation of the search algorithm and scoring function combination was performed by re-docking the native ligand to the binding site of the TraR receptor using the molecular docking approach. From the large number of conformations generated, the conformations that best reproduced the position and interactions of the native ligand within the binding site were chosen. This conformation was compared with the co-crystallized structure of the same native ligand, as is depicted in Figure 1. In the alignment of the conformations of the native ligand in (green color) and the resulting docking (pink color), it can be seen that the coupled structure located at the TraR binding site adopts a conformation almost identical to the co-crystallized structure of the native ligand and in this way, both the lactone ring and the acyl side chain is located suitably at the receptor binding site.

Figure 1. Alignment between the ligand present in the co-crystallized complex in green and the coupled ligand in pink, using the Autodock 4.2 program.

The mean quadratic root deviation (RMSD) value calculated by AutoDock was used to confirm whether the chosen conformation was correct or not. The RMSD value of the chosen conformation is 0.56 Å; generally, an RMSD value of less than 1.5 or 2 Å is considered to be optimal²⁸. In Figure 2, the interactions established in the conformation resulting from the molecular docking are visualized. In this ligand-receptor complex, the majority of hydrogen bond interactions that contribute to the binding of the native ligand and TraR are reproduced, as well as the correct folding of the TraR receptor²⁹. These interactions include hydrogen bonding between the Trp 57 residue and the carbonyl group of the lactone ring, hydrogen bonding between Asp 70 and the amino group and a bond between the second carbonyl group and Thr 129, with the exception of the bond between the Tyr 53 residue and the first carbonyl group located in the side acyl chain. Hydrophobic residues also contribute to the stable binding of the native ligand with TraR by hydrophobic interactions with residues Leu 40, Tyr 53, Tyr 61, Val 72, Trp 85. The interactions with Thr 51, Phe 62 and Ile 102 do not occur in said conformation.

Figure 2. Interactions present in the protein-native ligand complex.

Modes of binding among analogues and TraR

The structural interactions of this group of molecules were classified into five groups according to the differences present in the acyl chain reported by Geske et al.^11,12. This classification was made taking into account the importance of the influence of the variation of the acyl chain on the interactions that form the 34 analogues of the AHL within the TraR binding site. In addition to the classification by groups, the individual structure of each analogue was divided into two regions, A and B, as is shown in Figure 3, the first consisting of the lactone ring and the carbonyl group and the second region comprising the acyl chain, which is the one that varies between the different analogues. In this way, the analysis of the interactions is facilitated since each analog tends to form mostly the same interactions in the region A, unlike those formed in region B, which tend to be different³⁰.

Figure 3. Sketch of the 34 acyl homoserine lactones analogues.

Region A is formed by the lactone ring and the amide function and region B has variable side chains.

For each analogue-TraR docking, two or more clusters were generated. In these clusters, the different ways of binding can be identified. From the analysis of the clusters of each analogue-TraR docking, the best 34 conformations were chosen using the orientations that reported the lowest binding energy as the criterion of choice. The analogue-TraR complexes chosen, in comparison to 3OC₈HSL, are, in general, similarly oriented at the TraR receptor binding site. The structure of these analogues was constructed based on the structure of homoserine lactone, with variations in the acyl chain of the native ligand. Table 2 shows the IC₅₀ values and the free energies calculated using AutoDock for all 34 ligands.

Ligands 4, D6, C10, E22 and E35 have the highest values of inhibitory activity (IC₅₀), so the analysis of their structure-activity relationship for each group was taken as precedents. Analogues 4, C10, E22 and E35 form two hydrogen bonds with the TraR receptor; one between the carbonyl group of the lactone ring and the amino group in the Trp 57 residue and the second between the carbonyl group of the first carbon of the acyl chain and the residue Thr 129. Ligand D6 also forms these two hydrogen bonds, but also forms another with the Asp 70 residue and the carbonyl group of the acyl chain. The native ligand forms the same three hydrogen bonds mentioned that contribute to the bond at the TraR protein binding site. Similarly, the remaining ligands form these bonds but with the absence of one or two of them and differences in bond distances. The interactions of ligands D6 and E35 are detailed in Figure 4.

Figure 4. Interactions among ligands D6, E35 and TraR.

Hydrogen bonds are shown in blue, hydrophobic interactions in gray, π-π stacking in green and halogen bonds in cyan.

On the other hand, the possible hydrophobic interactions formed by these ligands in the hydrophobic pocket of TraR involve both the lactone ring and the acyl chain. Structurally, ligand 4 is classified within the group of non-native AHLs with an acyl chain of seven carbon atoms; one less than the native ligand 3OC₈HSL. This side chain allows the formation of eight possible hydrophobic interactions with Ala38, Leu40, Gln58, Phe62 and two interactions with each residue of Tyr61 and 63. In addition, in the lactone ring are the interactions with Val72 and Trip85, for a total of ten interactions. The D6 ligand has in its structure a side chain a carbonyl group on carbon 4 and a toluyl group. This structural modification places it within the group of other analogues, and between this chain and the protein it is possible to form 11 hydrophobic interactions; eight with the same residues as ligand 4 and also those formed with Ala49, Thr51 and Ala168. The C10 ligand has in its side chain a phenylacetanoyl group and E22 a phenoxyacetanoyl group, which places them in the phenylacetyl homoserine lactones and phenoxyacetyl homoserine lactones groups, respectively. Both form six hydrophobic interactions; four in the side chain with Leu40 and Tyr53 with different distances of connection, and two in the lactone ring with Val 72 and Trp85. Finally, of the 34 analogues, ligand E35 reported the highest inhibitory activity against TraR. It is classified in the phenylpropionyl homoserine lactones group since it has the phenylpropionyl group in its side chain and its possible hydrophobic interactions are the same as those mentioned for C10 and E22.

The differences among the possible hydrophobic interactions mentioned above for the different ligands representing each group and the interactions of the native ligand in the hydrophobic pocket of TraR are notable. However, these are not the only differences; some analogue-TraR complexes reported other possible interactions, such as π-π stacking, which are favored in ligands that in their structure possess an aromatic group, as can be seen in the complexes of ligands D6, C10, E22 and E35. The formation of this interaction occurs between the benzene rings of these molecules and the Tyr 61 residue. This type of interaction is also reported in other ligands, but with the amino acid Tyr53.

On the other hand, analogues that have halogen atoms in their side chains can form other types of interactions called halogen bonds. In this study, ligands C10, D15, E21, E23, B7, and E35 have halogen bond interactions with the Gln 58 residue and the E23 ligand with Trp 57. Studies of the characteristics of this type of halogen bond establish that they are strongly directional and the distance between the halogen and the nearby electronegative atom should ideally be less than the sum of its Van der Waals radii. Moreover, although the fluorine atom is a halogen, it does not have the capacity to form these interactions due to its great electronegativity³¹. According to the above, the halogen bonds of ligands D15, E23 and E35 are discarded as they are not possible since the interaction occurs between a fluorine atom and an oxygen atom. The sum of the van der Waals radii for the remaining ligands C10, E21 and B7 are calculated as follows: in the C10 ligand where its iodine atom and the Gln 58 oxygen interact, this value is 3.50 Å; in the E21 ligand where its bromine atom and the Gln 58 oxygen interact, the sum of the radii is 3.37 Å; and in ligand B7 where its bromine atom and oxygen of Gln 58 interact, the sum of the radii is 3.37 Å. The respective distances between the atoms bound for each ligand-receptor complex are 3.90 Å , 2.14 Å and 2.06 Å. Comparing these values, it follows that the only halogen bonds that are shorter than the value of the sum of the radii of the bound atoms are those of ligands E21 and B7 with the amino acid Gln 58.

It is quite clear that the addition of certain groups such as phenyl and halogenated substituents can generate other interactions, which could confer changes in the increase or decrease in the inhibitory activity of these analogues against the TraR receptor. These extra interactions (halogen bond, π-π stacking), formed by the analogues and the TraR receptor, are not present when TraR is in complex with the native ligand.

Analysis of the 3D-QSAR structure-activity quantitative relationship

For ligand-based alignment, ligand E35 was used as a template since it was the ligand with the highest inhibitory activity against TraR and for the receptor-based alignment, the best poses generated in the molecular docking, that is, those that reported the lowest binding energy and the correct orientation in the active site of the receptor, were used. Figure 5A and 5B show the superposition of all ligands for the ligand-based model (Model 2) and receptor-based model (Model 1), respectively.

Figure 5. Structural alignment.

A) Model 1 and B) Model 2.

The correlation analysis between CoMFA values and in vitro biological activities of the 34 analogues was performed by PLS and the evaluation of the predictive capacity of the models generated through the leave-one-out cross-validation method. Q² values of 0.457 and R² are obtained with cross validation of 0.856 for the Model 2, and a Q² of 0.697 and R² with cross validation of 0.781 for the Model 1. In effect, Model 1 shows a lower correlation so that their statistical confidence is lower than in model 2. Table 3 shows the other results of both PLS and cross-validation analyzes.

Figure 6 shows correlation of models 1 and 2, it can be seen that in model 2 there is a better pattern of adjustment of the points along the trend line, which translates into the prediction of biological activity values closer to those reported experimentally, unlike the graph for Model 2 where the residual dispersion is greater. From these comparisons, the Model 2 is statistically profiled as the one with the best correlation and the greatest predictive power. Even so, the Model 1 should not be discarded altogether since it also takes into account the structural similarity between the analogues and the native ligand; the position and orientation taken by the ligands at the receptor binding site according to the conformational freedom allowed by the amino acid residues that comprise it.

Figure 6. Predicted vs experimental activity scatter plot obtained by partial least squares analysis using comparative molecular field analysis (CoMFA).

Through the interpretation of the contour maps obtained by CoMFA, it is possible to know which structural characteristics of the ligands correlate with the changes in their biological activity values. These contours are represented by colors; steric contour maps are represented in green (region favorable for bulky groups) and yellow (region not favorable for bulky groups), while those of red (region favorable for electronegative substituents) and blue (region favorable for electropositive substituents) correspond to electrostatic contours³².

In Figure 7 the contour maps of the CoMFA models calculated are shown, taking as an illustrative example those generated for ligand D15, which records one of the highest values of inhibitory activity against TraR in the phenoxyacetyl homoserine lactones group. In these contours, it can be observed that the sterically favored regions are in positions close to the fluorine atoms found in the benzene ring substituent in the ligand structure in both Figure 7A and 7B. This is because this region has high electronegativity, which favors interactions with residues that have partially positive charges in their structures, while the sterically unfavorable regions, in yellow, spread around the phenyl group. However, on most of the hydrogen atoms in the phenyl ring, the positively charged regions that favor the electrostatic field are identified, as can be seen in more detail in Figure 7B. In the same way, carbonyl and ether groups with a negative charge are seen as electrostatically favorable regions due to the negative potential of oxygen atoms. In addition, the existence of halogens promotes the formation of halogen bonds with residues containing oxygen atoms in their structure. All these groups mentioned are identified as regions that would contribute to the antagonistic activity of ligand D15 against TraR.

Figure 7. 3D comparative molecular field analysis (CoMFA) contour maps for ligand 5.

A) receptor-based model and B) ligand-based model. Steric field: green: sterically favorable; yellow: sterically unfavorable. Electrostatic field: red: negative potential favorable; blue: negative potential unfavorable.

Underlying data

Figshare: Receptor_Alim_2019_Pose_E_D_1.mdb (1).zip. https://doi.org/10.6084/m9.figshare.10255739.v1²³.

This project contains the following underlying data:

Data are available under the terms of the Creative Commons Attribution 4.0 International license (CC-BY 4.0).