Insect RDL Receptor Models for Virtual Screening: Impact of the Template Conformational State in Pentameric Ligand-Gated Ion Channels.

The RDL receptor is one of the most relevant protein targets for insecticide molecules. It belongs to the pentameric ligand-gated ion channel (pLGIC) family. Given that the experimental structures of pLGICs are difficult to obtain, homology modeling has been extensively used for these proteins, particularly for the RDL receptor. However, no detailed assessments of the usefulness of homology models for virtual screening (VS) have been carried out for pLGICs. The aim of this study was to evaluate which are the determinant factors for a good VS performance using RDL homology models, specially analyzing the impact of the template conformational state. Fifteen RDL homology models were obtained based on different pLGIC templates representing the closed, open, and desensitized states. A retrospective VS process was performed on each model, and their performance in the prioritization of active ligands was assessed. In addition, the three best-performing models among each of the conformations were subjected to molecular dynamics simulations (MDS) in complex with a representative active ligand. The models showed variations in their VS performance parameters that were related to the structural properties of the binding site. VS performance tended to improve in more constricted binding cavities. The best performance was obtained with a model based on a template in the closed conformation. MDS confirmed that the closed model was the one that best represented the interactions with an active ligand. These results imply that different templates should be evaluated and the structural variations between their channel conformational states should be specially examined, providing guidelines for the application of homology modeling for VS in other proteins of the pLGIC family.

Introduction

Pentameric ligand-gated ion channels (pLGICs) constitute a large family of transmembrane receptors widely expressed in animals, from insects to mammals, as well as in a few bacterial and archaeal species.[1] They mediate fast synaptic transmission in the central and peripheral nervous system, allosterically converting the binding of a neurotransmitter in their extracellular domain (ECD) to the opening of an ionotropic pore in the transmembrane domain (TMD).[2] The pLGIC family include γ-aminobutyric acid receptors (GABAA-Rs), the nicotinic acetylcholine receptor (nAChR), glycine receptors (GlyR), the serotonin 5-HT3 receptor, the eukaryotic glutamate receptor (GluCl), and the bacterial Erwinia chrysanthemi and Gloeobacter violaceus ligand-gated ion channels (ELIC and GLIC, respectively).[3] These proteins share an evolutionarily conserved cylinder-shaped architecture, in which five identical or different subunits are arranged around a central five-fold axis. Each subunit consists of a large ECD containing the agonist binding site, a TMD with four α-helices (labeled from M1 to M4), which encloses the channel pore, and an intracellular domain (ICD) between M3 and M4. Besides the binding of their specific neurotransmitters, pLGICs are important targets for many pharmacological drugs and, in the case of invertebrates, insecticides.[4]

GABAA-Rs are the specific target for γ-aminobutyric acid (GABA), the major inhibitory neurotransmitter in both vertebrate and invertebrate central nervous systems (CNSs).[5] The binding of GABA selectively regulates the gating of chloride anions through the channel pore in neurological synapses.[6] In vertebrates, a diversity of GABAA-Rs subtypes exist due to the assembly of different combinations of 19 subunit isoforms. The most abundant GABAA-R class in the human CNS is a heteropentamer with two α subunits, two β subunits, and one γ subunit.[7] On the contrary, in insects, the most representative class of GABAA-R is a homopentamer, formed by the so-called RDL subunit. RDL homopentamers are related to their vertebrate counterparts in that their activation opens a chloride channel, but they show significant differences in the pharmacological profile.[8] They are widely expressed in the insect CNS, developing essential functions such as learning, as well as visual and olfactory processing, so their blockage would lead to the insect death. Therefore, RDL receptors constitute a relevant target for insecticides.[8,9] Different sites in this receptor are known to bind several families of insecticidal compounds.[10] The noncompetitive antagonist type IA (NCA-IA) site, located inside the transmembrane pore (Figure ), is the target for many insecticides that block the channel, ranging from natural compounds to synthetic chemicals, including picrotoxinin, polychlorocycloalkanes (e.g., lindane and toxaphene), cyclodienes (e.g., dieldrin and α-endosulfan), and the phenylpyrazole fipronil.[10−12] The NCA-IA site is the most studied RDL binding site and the one with the largest amount of known active ligands.

Figure 1

Overall structure of an RDL homology model (template: PDB ID: 3RHW) (A) Side view of the homopentamer, with each of the subunits in a different color. ECD: extracellular domain. TMD: transmembrane domain. (B) Bottom view of the receptor from the intracellular region. The TMD of each subunit is composed of four α-helices, named M1 to M4. The M2 α-helices (shown in orange) surround the channel pore. (C) Sequence of the M2 residues that constitute the NCA-IA binding site, with their corresponding index numbers (above). The asterisks below the sequence indicate the residues that face the channel pore and have been reported as the key residues for the interaction with insecticides.[11,12,21] (D) Detailed view of the key residues of the NCA-IA site with their corresponding index numbers. Only two subunits of the RDL receptor are shown.

Because of the difficulties related to the expression and crystallization of large-size membrane proteins, the determination of the 3D structures of pLGICs has remained elusive. Indeed, there are only 23 unique structures deposited in the Protein Data Bank that belong to this family.[13] Considering the lack of experimentally obtained structures for the entire pLGIC family, there has been a huge development of homology models for these proteins.[14] Specifically, the 3D structure of the insect RDL homopentamer has not been determined yet. For that reason, several studies have developed RDL homology models to study their interactions with agonists in the orthosteric site[15] and with insecticide molecules in allosteric sites. The first models were obtained by using the bacterial ELIC and GLIC as templates.[16,17] Later models were based on the Caenorhabditis elegans glutamate-gated chloride channel (GluCl),[17−23] the human GABAA-R β3 homopentamer,[24−26] the human Gly-R α3 homopentamer,[27,28] and the human GABAA-R α1β3γ2 heteropentamer.[29] With these models, ligand–protein docking studies and molecular dynamics simulations (MDS) have been conducted, shedding light on the interactions between RDL and insecticides, as well as the effects of some specific mutations on the affinity for these ligands. Also, recently, a structure-based virtual screening (VS) process of potential antagonists binding to an intersubunit site was performed in an RDL homology model based on the GABAA-R β3 homopentamer.[30]

Despite the considerable usefulness of homology modeling, this methodology presents some caveats in the case of GABAA-Rs and pLGICs in general. In particular, pLGICs are highly dynamic, undergoing discrete allosteric transitions between at least three conformational states.[3] These conformations are characterized by a specific channel gating state. In the resting state, the channel is closed by hydrophobic residues in the extracellular end of the TMD, while in the activated form, the channel is open, allowing an ion current flow.[31] After prolonged exposure to agonists, the receptors become desensitized, a state characterized by a closed channel, although in a different conformation from that of the resting state.[32] The above-mentioned homology models of the RDL homopentamer were obtained from templates whose structures had been determined in different conformational states. As an example, the C. elegans GluCl[33] presents an open conformation, while the GABAA-R β3 homopentamer is in a desensitized state.[34] However, there has been no detailed investigation of the impact that these differences could exert on molecular docking and VS performance.

In order to study any binding site structure in detail and to perform molecular docking or VS using homology models, it is necessary to assess which of the available structures are the most suitable templates and to clarify if several templates have to be considered to study the consequences of protein motion on the site of interest.[35] An increasingly employed strategy to assess homology models consists of using molecular docking to evaluate if the binding site can identify known active compounds among decoys. A model that displays high enrichment of known ligands can be considered as a good representative of the receptor structure and suitable for VS. Moreover, the VS performance of structurally different homology models of the same protein can be correlated with properties of interest, contributing to the understanding of the structural features of the system.[36−38] This retrospective VS approach has been extensively used for the assessment of G protein-coupled receptors (GPCRs) homology models, but no intensive studies have been found to apply it in pLGICs.[36−39]

Considering that the RDL homopentamer NCA-IA site is located inside the transmembrane pore, we hypothesize that the conformational state of the pLGIC which is used as a template for obtaining a homology model of this protein could be relevant for the performance of VS. To evaluate this, we obtained 15 homology models of the RDL receptor, based on pLGIC templates in different conformational states. A set of known active ligands, combined with decoys, were docked to the binding site of each model, and the retrospective VS performance was assessed. Then, we correlated these results with the structural properties of the channel to analyze which are the determinant factors for the correct identification of active ligands. Finally, we performed MDS with three of the best performing models in complex with the docked pose of insecticide fipronil, each one in a different conformational state, to assess the possible differences in the interactions with this representative ligand. Our work aims not only to develop RDL homology models suitable for VS of new insecticides but also to provide insights and methodological guidelines that could be useful for the application and assessment of homology models for VS in other proteins of the pLGIC family.

Results

Construction of RDL Homology Models

Homology models of the insect Musca domestica RDL homopentamer were generated based on 15 pLGIC templates whose structures had been determined in different conformational states of the channel: four in the open state, five in the closed state, and six in the desensitized state (Table ). The sequence alignments were facilitated by the high structural conservation of the pLGIC family, specially in the TMD. Some manual adjustments were made in the alignments of the ECDs to avoid the placing of gaps inside the α-helix or β-sheet regions. The target sequence of the M. domestica RDL subunit is shown in Figure S1, and the final alignments between the target and each of the template sequences are shown in Figure S2. The 15 templates present a similar global sequence identity (Seq Id) with the target, of around 40–45%, although for the two GLIC templates the percentage is lower (22.4%). In the case of the TMD, the Seq Id is slightly higher for all the templates, while in the M2 α-helix region, where the binding site under study is located, the templates cover a wider range of Seq Ids, going from 31.8 to 70.8%.

Table 1

Experimental Structures of the pLGICs Selected as Templates for the Homology Modeling of the RDL Receptor

			Seq Id (%)
protein	PDB codea	conformational state of the channelb	globalc	TMDd	M2e
GLIC	3EAM	open[46]	22.4	23	31.8
α1 glycine receptor	3JAD	closed[47]	44.6	57.1	66.7
α1 glycine receptor	3JAE	open[47]	44.6	57.1	66.7
α1 glycine receptor	3JAF	desensitized[47]	44.6	57.1	66.7
glutamate-gated chloride channel	3RHW	open[33]	40.5	46.1	58.3
β3 GABAA receptor	4COF	desensitized[34]	45.9	55.3	58.3
glutamate-gated chloride channel	4TNV	closed[31]	40.5	46.1	58.3
α3 glycine receptor	5CFB	closed[40]	44.7	55.5	70.8
α3 glycine receptor	5TIO	desensitized[41]	44.7	55.5	70.8
α3 glycine receptor	5VDH	desensitized[42]	44.7	55.5	70.8
chimeric α1 GABAA-R/ELIC	6CDU	desensitized[43]	37.0	50.4	62.5
GLIC	6HY9	closed[44]	22.4	23	31.8
α1 glycine receptor	6UD3	open[45]	44.6	57.1	66.7
α1 glycine receptor	6UBS	closed[45]	44.6	57.1	66.7
α1 glycine receptor	6UBT	desensitized[45]	44.6	57.1	66.7

Code of the protein used as the template in the Protein Data Bank.

Conformational state of the channel assigned to the protein structure (see references).

Seq Id between the template and the whole target sequence (including both the ECD and the TMD).

Seq Id between the template and the target TMD, as defined in Figure S1.

Seq Id between the template M2 α-helix and the target M2 α-helix, as defined in Figure S1.

As a measure of the structural quality of the models, QMEAN6 scores were calculated and are displayed in Table S1. This score is a linear combination of six structural descriptors, and a higher value (ranging between 0 and 1) reflects a stronger reliability for the model.[48] For comparison, the respective template structures were also evaluated by QMEAN6. In most cases, the homology model scores are comparable to those of the templates from which they were built. Additionally, PROCHECK Ramachandran plots were obtained for the models, demonstrating the good quality of the backbone geometries. In all cases, more than 90% of the residues were located in the most favored regions of the plot, while less than 0.7% were located in disallowed regions. The root mean square deviations (rmsds) between each model and its template were also calculated. Considering the whole proteins, rmsd values were very low, going from 0.2 to 0.47 Å, whereas when taking into account only the M2 region, the values were even lower, from 0.04 to 0.26 Å. This structural similarity allows us to assign the RDL homology models to the same conformational state as that of the template they came from.

To characterize the common structural features of the closed, open, and desensitized states at the NCA-IA binding site, Figure shows the pore diameter profile, obtained using the HOLE[49] software, of three representative RDL homology models. Closed models present a cylindrical arrangement that becomes progressively funnel-shaped in the open and desensitized structures, in which the channel tapers toward the intracellular end. Moreover, in closed models, the pore is narrower, with a prominent constriction at 9′Leu. In the open structures, the pore diameter profile shows an expanded pathway for ion permeation, without substantial constriction points. Finally, in the desensitized state, there is an important constriction in the intracellular end of M2, at the level of -2′Pro. However, it is important to consider that the diameter profiles of different models, even if they are assigned to the same conformational state, are not identical. As can be seen in Figure S3, they present variable values of pore diameter at key residues of the binding site. These differences are more pronounced among open and closed models, in comparison to the desensitized ones, that share more similar profiles. In the case of the closed state, all the models share the constriction at 9′Leu, a particularly relevant feature because it constitutes the gate that controls ion permeation.

Figure 2

Pore diameter profiles of representative RDL homology models in the closed (RDLC, blue), open (RDLO, red), and desensitized (RDLD, green) conformational states. The closed model was based on template 6UBS, the open model on 3RHW, and the desensitized model on 6UBT. The locations of the binding site key residues that face the channel lumen are indicated on the right of the plot.

Molecular Docking and Retrospective VS Performance

Molecular docking screens of active ligands and decoys were carried out to evaluate the VS performance of our RDL homology models in different channel conformational states. The set of active ligands included 42 compounds that, when appropriate, were docked in different tautomeric, chiral, and ring conformational forms, giving a total of 194 variants. On the other hand, the property-matched decoy set included 1734 molecules that were selected after filtering an initial larger set obtained from DUD-E to match as much as possible the active ligands physicochemical properties (see the Computational Methods section). The decoys were subjected to the same variant generation protocol used for the active ligand set, giving a total of 5729 variants. For the subsequent assessment of VS performance, only the best scored variant of each ligand was considered.

Following the molecular docking runs, the prioritization of active ligands over decoys was assessed for each of the homology models by the receiver operating characteristic (ROC) curve. The area under the ROC curve (AUC) was calculated as a measure of the general selection of active ligands over decoys across the entire screen. In addition, the Boltzmann-enhanced discrimination of the ROC (BEDROC) score was also calculated to evaluate the early enrichment of active ligands. This parameter is specially relevant because the ultimate aim of a VS protocol is to select a subset of the top scoring compounds for further experimental validation.

Considering that the molecular docking performance is strongly dependent on the sampling algorithm and the scoring function used,[50] three docking tools were used: AutoDock Vina, 2Vinardo, and LeDock. Figure compares the retrospective VS performance parameters obtained using AutoDock Vina for each of the RDL homology models. In terms of the AUC, almost all the models performed better than random selection, presenting AUC values above 0.5, except for the models based on templates 3JAE and 3JAF. The AUC results ranged from 0.46 to a maximum of 0.76 for the model based on template 6UBS, a closed structure for the α1-glycine receptor. In the case of BEDROC (α = 20), there was a wider variation among the models, considering that five of them performed worse than random selection, with BEDROC (α = 20) values below 0.06 (see the Computational Methods section). BEDROC (α = 20) results for the other 10 models ranged from 0.07 to a maximum of 0.30, also corresponding to the model based on the 6UBS template. With regard to BEDROC with α = 100, a more challenging parameter given that it focuses on the top 1.6% of the ranking, only four models performed better than random selection. The model based on the 6UBS template was also the one that resulted in the highest value of this metric.

Figure 3

Retrospective VS performance parameters obtained using AutoDock Vina for RDL homology models. The bars represent the AUC value (A), the BEDROC (α = 20) score (B), and the BEDROC (α = 100) score (C). PDB codes of the templates from which the homology models were obtained are indicated in the x axes (for references, see Table ). Bar colors represent the conformational state of the models (blue: closed, red: open, and green: desensitized). For each plot, the bars were sorted from the left to the right following the improvement of the corresponding parameter. The values of the parameters for each model are indicated above the bars.

It is apparent from the plots in Figure that homology models presenting an open conformational state of the channel tended to perform worse than closed and desensitized ones. This trend becomes more evident when the early enrichment parameter BEDROC is evaluated, given that all the open models performed worse than random selection. Taking together the AUC and BEDROC results, the models based on templates 6UBS and 6UBT can be considered as the best performing ones. The former presents a closed channel, while the latter corresponds to a desensitized structure. The model based on the 3RHW template was the best performing among the open models, although with poor results.

The 2Vinardo performance parameters for each homology model are presented in Figure S4. It can be seen that the AUC values are consistently lower than the values obtained using AutoDock Vina. The same trend was observed in terms of BEDROC, where only two models resulted in values above random expectation. In addition, Figure S5 depicts the performance of each model that was obtained using LeDock. The AUC values were also lower than AutoDock Vina results in all the models, while BEDROC values were comparable with α = 20 and higher with α = 100. Although the performances of 2Vinardo and LeDock were, in general terms, worse than the results obtained using AutoDock Vina, it is important to highlight that the outcomes of the three scoring functions agree in two main aspects. First, in the three cases, the results show a general trend that favors closed and desensitized models over the open structures. This indicates that the three scoring functions showed differences in performance that were determined by the conformational states of the models. In addition, in all cases, the best performing model was the one based on the 6UBS (closed) template. For subsequent analyses of the properties that explain the variations of VS performance in our RDL homology models, the AutoDock Vina results were used, considering the better global performance of this scoring function in the system under study.

Relation between Retrospective VS Performance and Target/Template Seq Id

To analyze which are the determinant factors for a successful VS campaign on RDL homology models, we first evaluated the possible correlations of the AUC and BEDROC with the target/template Seq Ids. This evaluation was conducted by considering the whole sequences of the target and each of the templates, taking only the TMD sequences or including exclusively the residues of the M2 region, where the NCA-IA binding site is located.

There were no statistically significant correlations (p > 0.05) between target/template Seq Ids and the homology models retrospective VS performance parameters (Figure S6). Regarding the global Seq Id, homology models with similar Seq Ids, of around 40–45%, resulted in a wide range of AUC and BEDROC values. Moreover, the results of the two models with the lowest global Seq Id fell in the same range as the rest of the structures. This lack of correlation was also observed when only the TMD sequences were considered. In the case of the M2 region, even though the Seq Ids between the target and most of the templates are higher, this variable did not explain the differences in the VS performance.

Relation between Retrospective VS Performance and Structural Properties of the Channel

Considering the lack of correlation between RDL homology models VS performance and Seq Ids, as well as the observed trend that favors closed and desensitized models over the open ones, we now turn to analyze the impact of some structural properties of the binding cavity on the prioritization of active ligands in the docking screens. First, a correlation analysis was conducted to evaluate the relation of the AUC and BEDROC with the solvent-accessible (SA) area and SA volume of the binding site in each of the models (Figure ). It can be observed in the scatter plots that the AUC and BEDROC values tend to decrease as the size of the cavity increases, both in the SA area and in the SA volume. Pearson correlation coefficients (r) were found to be significant (p < 0.05) and showed negative values in all cases. The strongest correlation (r = −0.84) was the one between the AUC and SA area. No clear differences were observed between groups of models corresponding to each of the three conformational states. However, it has to be considered that, as previously mentioned, the sizes of the channel at the binding site are not identical among models assigned to the same state (Figure S3). For this reason, the relation of the screening performance and the pore diameter at the key residues was further analyzed.

Figure 4

Correlation analysis between VS performance and SA area and volume of the binding cavity, as obtained from CASTp 3.0.[51] To account for the screening performance, both AUC (upper panels) and BEDROC α = 20 (lower panels) values were included. In panels A and C, the SA area is represented, while B and D panels refer to the SA volume. Each circle represents a RDL homology model based on a different template, and the colors indicate its conformational state (blue: closed, red: open, and green: desensitized). The Pearson correlation coefficient (r) is indicated for each case in the upper left corner of the plot.

Figure shows the relation between the VS performance and the pore diameter at the four residues that face the channel pore in the RDL NCA-IA binding site: −2′Pro, 2′Ala, 6′Thr, and 9′Leu. The scatter plots not only reflect the common structural features of each conformational state, such as the constriction in −2′Pro for desensitized models and in 9′Leu for closed channels, but also exhibit the differences in pore diameter among models assigned to the same state. The most interesting aspect of these plots is that the inverse correlation found between the screening performance and the size of the cavity can be mainly attributed to the spatial arrangement of some specific residues. In particular, remarkable negative correlations were found between the AUC and BEDROC with the pore diameter at 2′Ala (r = −0.85 and r = −0.74, respectively). Pearson coefficients were also high for the correlations between the pore diameter at −2′Pro and screening performance parameters. On the contrary, the correlations of the AUC and BEDROC with the pore diameter at 6′Thr and 9′Leu were not significant. The hydrophobic interactions with −2′Pro and 2′Ala seem to be specially relevant for the identification of active ligands in the docking screens.

Figure 5

Correlation analysis between the VS performance parameters (AUC and BEDROC α = 20) and the pore diameter at four residues located in the NCA-IA binding site: −2′Pro (panels A and E), 2′Ala (panels B and F), 6′Thr (panels C and G), and 9′Leu (panels D and H). Each circle represents an RDL homology model based on a different template, and the colors indicate its conformational state (blue: closed, red: open, and green: desensitized). The Pearson correlation coefficient (r) is indicated for each case in the upper left corner of the plot.

Molecular Dynamics Simulations

Three homology models of the RDL receptor, based on the templates in the closed (6UBS, RDLC), open (3RHW, RDLO), and desensitized (6UBT, RDLD) states, were subjected to MDS to assess their dynamic behavior in the presence of fipronil, a canonical NCA-IA insecticide. The initial coordinates were taken from the docking poses obtained using AutoDock Vina. These three models were selected because they were the best performing among models corresponding to each of the conformational states. It is important to highlight that in the three cases, fipronil was ranked among the top 10 active ligands in the docking screens. Each system was embedded in a hydrated 1-palmitoyl-2-oleylphosphatidylcholine (POPC) lipid bilayer, equilibrated, and simulated for 150 ns.

The rmsd from the initial structures was plotted versus the simulation time to evaluate the stability of the models. The protein alpha-carbon (αC) rmsd showed a convergence value of 3–3.5 Å (Figure A). This value becomes lower when only the αC units of the binding sites are considered, with the closed template showing rmsd values for the binding site considerably lower than those of the open and the desensitized ones (Figure B). Figure C shows the superposition of fipronil in several snapshots taken along the trajectories. Fipronil remained bound to its binding site in the three models along the simulation time, with less orientation changes in the closed and desensitized states in comparison to the open model. The convergence of the rmsd for the αC and the binding stability observed for fipronil led us to conclude that the RDL homology models in complex with the insecticide were stable along the simulated trajectories. The pore diameter profiles were also calculated along the MDS (Figure S7). Some differences in the diameter profiles were observed between the homology models and the initial structures used for MDS. These variations are due to the movements of the amino acid side chains that occurred during the equilibration stage. It can be seen that in the closed model, the pore maintained almost the same diameter profile of the initial structure, while the open and desensitized models showed marked changes, with their pore diameters tending to decrease. These changes in the pore diameter, particularly in the open model, may be related both to the instability of the protein structure and to the effect of the interactions between the channel pore residues and fipronil.

Figure 6

(A) rmsd vs the simulation time for the protein αC of the three simulated RDL homology models, in the closed (RDLC), open (RDLO), and desensitized (RDLD) conformational states. (B) rmsd vs simulation time for the binding site αC. (C) Superposition of the fipronil molecule orientation in snapshots taken at 10 ns (red), 25 ns (green), 50 ns (blue), 75 ns (magenta), 100 ns (orange), 125 ns (black), and 150 ns (cyan) for the three simulated models in the closed (left), open (middle), and desensitized (right) states.

To compare the binding modes of fipronil in the three RDL homology models, its interactions with the binding site residues were assessed in the initial docking poses, as well as along the MD trajectories. Previous experimental and computational studies have reported that the key residues in the NCA-IA site that interact with fipronil are −2′Pro, 2′Ala, 6′Thr, and 9′Leu.[11,18,21] −2′Pro, 2′Ala, and 9′Leu present mainly hydrophobic interactions, while 6′Thr has been reported to be very relevant for the binding of fipronil because it forms hydrogen bonds (H-bonds) with the insecticide molecule.[11,18,22] The comparison of the docking poses shows that the fipronil orientation is different according to the conformational state of the channel. It presented a tilted orientation in the open and desensitized models, while it was vertically erected in the closed model (Figure S8). Fipronil formed several of the expected hydrophobic interactions with RDL residues, specially in the desensitized and closed models. However, none of the docking poses showed the expected H-bonds with 6′Thr. On the contrary, the formation of these H-bonds could be observed in the MD trajectories.

The analysis of the H-bonds between fipronil and 6′Thr along the MDS is shown in Figure . In the open RDL model, the H-bonds were formed intermittently and none of these interactions remained stable in the final part of the trajectory (Figure A). The desensitized model presented even more unstable H-bonds, being mostly absent along the whole simulation time. Lastly, the closed model showed the most persistent H-bonds, with two stable interactions in the first 50 ns and one that remains until the end of the trajectory. The occupancy plot (Figure B) indicates that the more stable H-bond is formed with the oxygen atom of fipronil acting as the acceptor and the 6′Thr side chain as the donor. The other, less stable, H-bond is formed with the NH2 group of fipronil as the donor and the 6′Thr side chain as the acceptor. Figure D shows an MD simulation snapshot depicting the two H-bonds formed between fipronil atoms and two 6′Thr residues from adjacent subunits.

Figure 7

Analysis of the H-bonds between fipronil and 6′Thr in RDL homology models along the MD trajectories. (A) Number of H-bonds formed between fipronil atoms and the 6′Thr side chain in each of the three RDL homology models along the simulation time. (B) Occupancy percentage for each of the fipronil atoms that form H-bonds with 6′Thr −OH. The atoms of fipronil included in this plot are marked in red in the molecular structure depicted in (C). (D) Representative snapshot of RDL6UBS MD simulation showing in detail two H-bonds: one with fipronil O acting as the donor and the 6′Thr side chain as the acceptor and the other one with fipronil NH2 as the donor and the 6′Thr side chain as the acceptor. The 6′Thr residues belong to two different adjacent subunits. H-bonds are depicted in detail with green thick lines.

Persistent H-bonds between fipronil and 6′Thr are expected to contribute to the binding energy. Calculations of the interaction energies between RDL models and fipronil were carried out using gmx energy for the last 100 ns of the trajectories. The average interaction energy for fipronil in RDLO was −168.1 ± 2.3 kJ/mol, in RDLD was −158.7 ± 5.0 kJ/mol, and in RDLC was −175.6 ± 1.2 kJ/mol. Therefore, MD simulations showed that the binding of fipronil resulted to be more energetically favorable in the closed model than that in the open and desensitized ones. The contributions of M2 residues to the interaction energies were also calculated, considering both the electrostatic (calculated by the Coulomb potential) and the van der Waals (calculated by the Lennard-Jones potential) forces (Figure ). The plot clearly indicates the relevant contribution of the electrostatic forces in the interaction of fipronil with 6′Thr in the closed model. As expected from the analysis of the H-bonds, this energy is lower in the open model, while it is insignificant in the desensitized state. The interaction energies with −2′Pro and 2′Ala are similar among the three models and mostly represented by van der Waals forces. In the case of 9′Leu, the energies are considerably higher in the open and closed models than that in the desensitized one, given that the huge pore diameter at this residue in this conformational state may have impeded a stronger interaction with fipronil. In agreement with previous studies,[18] some residues that do not face the channel pore, such as 5′Val, also showed relevant van der Waals interactions with fipronil.

Figure 8

Average interaction energies between fipronil and the M2 residues in the binding site. The contribution of electrostatic (calculated by the Coulomb potential) and van der Waals [calculated by the Lennard-Jones (LJ) potential] forces are shown for each of the RDL homology models.

Discussion and Conclusions

Homology modeling approaches have been intensively used for gaining knowledge about the pLGIC structures, in order to elucidate physiologically relevant aspects, such as their gating mechanism[52] and their interactions with ligands.[53] Particularly, in the case of GABAA receptors, efforts have been made to obtain good-quality homology models,[14] both for vertebrate subunits and for the insect RDL homopentamer. The most common assumption in homology modeling is that the template with the highest Seq Id to the target should be selected because it would result in the best obtainable model. However, in a previous work, the evaluation of different templates has proved to be useful for achieving successful VS in the allosteric sites of the vertebrate GABAA receptor.[54] Indeed, one of the most challenging aspects in the homology modeling of GABAA receptors is the highly dynamic behavior of this ion channel. The conformational state of the template on which the models are based becomes specially relevant in the case of the insecticide NCA-IA binding site in the RDL receptor, given that it is located inside the channel pore.

In this work, the retrospective VS performance of RDL homology models based on different pLGIC templates was assessed. We evaluated which were the determinant factors for the correct identification of active ligands, by correlating the obtained performance parameters with some relevant properties of our models. The most interesting finding was that the RDL homology models showed differences in VS performance according to their conformational states. With the three scoring functions used, closed and desensitized models tended to present better results than the open ones. Moreover, the three scoring functions agreed on the fact that the best performing model corresponded to a closed structure, based on the 6UBS template.

A note of caution is due here since even the best performing model presented rather modest enrichment metrics. About this, it is important to consider that AUC values of around 0.75 and BEDROC scores of around 0.2–0.3 are still acceptable. These results do not interfere with the aim of this work, considering that the purpose was not to adjust a methodology in order to obtain the best achievable metrics for VS performance in RDL homology models but to analyze the relations between relevant features of the models and their retrospective VS performance, specially focusing on the template conformational state.

The correlations between the VS performance parameters and target–template Seq Id were not statistically significant, even when only the binding site sequences were considered. Although this observation is contrary to the rule of thumb for template selection in homology modeling, it is not entirely surprising, considering that previous studies with GPCRs arrived to similar conclusions, demonstrating that the target/template Seq Id is not always a good predictor of VS performance.[36−38] Also, this finding may be somewhat limited by the fact that almost all the templates used in this work presented a very similar global Seq Id to the target, of about 40–45%, and templates with higher identities were not available. Therefore, a thorough exploration with a wider spectrum of template Seq Ids with the RDL receptor was not possible to achieve.

On the other hand, significant linear correlations were found between some structural properties of the channel and VS performance parameters. First, negative correlations were observed between the SA area and volume of the binding site cavity and the screening performances of the RDL models, meaning that molecular docking was more effective at the active ligands scoring in models with a more constricted cavity. Besides, to analyze if this relation could be assigned to the spatial location of some specific residues, the correlations of the ROC AUC and BEDROC with the pore diameter were also evaluated. Significant negative correlations were found between these performance parameters and the pore diameter at −2′Pro and 2′Ala. Thus, AUC and BEDROC values were particularly higher in RDL models with more constricted pore diameters at the level of these residues. These results can explain the bad performance of the closed and desensitized models based on templates 3JAD, 5CFB, and 3JAF (see Figure ) since these models present wide binding site cavities, with the highest pore diameters at the key residues among their corresponding groups of models (see Figure S3).

The strongest correlations between VS performance and pore diameter were found at 2′Ala. This result points to the key role of this residue in the recognition of active ligands. This is in agreement with a large body of experimental evidence that has demonstrated the relevance of this residue for the hydrophobic interactions between RDL and NCA-IA ligands, given that a naturally occurring mutation at this amino acid leads to the development of resistance for the action of insecticides.[55,56]

Several previous studies have used templates in the open conformational state to build RDL homology models for performing molecular docking and MD simulations with NCA-IA insecticides.[16,19,22,23] This choice may be based on the fact that under physiological conditions, the open state is the one that facilitates the channel blocking by these antagonists.[57,58] In our docking screens, however, a closed model (based on template 6UBS) was the one that obtained the highest AUC and BEDROC values, while all the open structures tended to perform worse. This result can be explained by considering that a more constricted binding cavity, such as the one that the 6UBS-based model presents, may be the most representative structure of the blocked channel, which is the final state produced by the binding of NCA-I insecticides. Therefore, molecular docking could be more effectively describing the interactions with the active ligands in this type of structures than that in the models with wider cavities. What is somewhat surprising is that in general terms, the desensitized models resulted in good VS performance. Although not physiologically meaningful, the good performance of this state can be explained by its structural features. The desensitized state is characterized by a marked constriction in −2′Pro, and it is also partly constricted at 2′Ala. Thus, it satisfies the determinant features for successful VS that were observed in the correlations with the pore diameter. Nonetheless, this analysis on the desensitized state can be enriched in the light of our MD simulations.

MD simulations were performed for three complexes between RDL and the insecticide fipronil, each one representing the closed, open, and desensitized states. Fipronil was the active ligand selected for the simulations because there is a large body of previous evidence about its interactions with the RDL receptor.[19,22,59] None of the docking poses for fipronil, that were used as the initial structures for MD simulations, presented the expected H-bonds with 6′Thr.[19,60] This may be caused by the fact that H-bonds can be more difficult to identify than hydrophobic contacts in protein–ligand docking using AutoDock Vina.[61] However, the MD trajectories showed that these H-bonds did form in the closed model, while they were unstable in the open model and did not appear in the desensitized state. This shows that although the desensitized models resulted in acceptable VS performance parameters, when the best performing model among this group was subjected to MD simulations, it was not able to reproduce the expected interactions with a representative ligand.

It is important to highlight that the H-bonds with 6′Thr influenced the interaction energies between each model and fipronil, resulting in a more favorable energy for the closed model. In a previous work, Zheng et al. performed molecular docking and MDS of fipronil in complex with an RDL homology model in the open conformational state.[19] They found that a H-bond with 6′Thr was formed in the docking pose, but it disappeared during the subsequent simulation. The comparison of this previous result with our findings supports the conclusion that the results of computational studies on the interactions between NCA-IA insecticides and RDL homology models may vary according to the conformational state of the template on which the models are based. Moreover, a closed template may be the most suitable and representative conformational state, not only for performing VS but also for carrying out MDS to analyze the interactions with known ligands.

Although the present work is a case study about the insect RDL receptor and the results are particularly valid for the NCA-I insecticide binding site, our methodology and results have relevance for future computational studies in GABAA receptors in general and other pLGICs as well, given that these proteins share a common global architecture. Particularly, we expect our study to be useful for studying binding sites located in the TMD of pLGICs, which can be remarkably affected by the channel conformational state. Our results suggest some guidelines for the pLGIC homology modeling for VS. The Seq Id between the template and the target may not be a good predictor of VS performance. Therefore, different templates should be evaluated, and the structural variations between their conformational states should be examined. Moreover, considering the dynamic behavior of these ion channels, templates in the conformational state that better represent the final conformation produced by ligand binding should be prioritized for building the target structures.

Computational Methods

Selection of Templates and Homology Modeling

The target sequence of the M. domestica RDL subunit was retrieved from UniProtKB (accession number: Q75NA5) (Figure S1). The sequence of the predicted cytoplasmic loop of the ICD of the RDL subunit was eliminated due to the lack of suitable templates for modeling it.[19] Fifteen homopentameric pLGICs, whose structures are available in the Protein Data Bank, were selected as templates. This selection aimed to include templates in different conformational states (closed, open, and desensitized), regardless of their Seq Id with the target protein (see Table ). The sequence alignments between the target and each of the templates were performed with the align2d() function of Modeller v9.24.[62] For each template, an iterative alignment–modeling–evaluation approach was performed to avoid the placing of gaps inside α-helices or β-sheets.[63] With the final alignments, 400 models were generated in Modeller, with a slow refinement level. The quality of the models was evaluated following a consensus approach,[64] combining the Modeller built-in objective (molpdf) and discrete optimized protein energy functions, the ProSA z-score[65] and the QMEAN6 normalized score.[66] The percentages of residues located in the disallowed regions of the Ramachandran plot, obtained from PROCHECK,[67] were also considered. The best-quality model for each template was selected for further studies.

Active Ligand and Decoy Sets

A set of known active ligands (IC50 < 1 μM) for the NCA-IA site of the M. domestica RDL receptor was built, including 42 compounds retrieved from ChEMBL[68] and collected from the literature.[69−73] The selection of active ligands aimed to include compounds with a reported IC50 value for the displacement of [3H]-EBOB, a specific radioligand probe for the RDL NCA-IA site.[12] The list of active compounds is shown in Table S2. 50 property-matched decoys were obtained for each of these ligands by using the DUD-E approach.[74] To reduce the artificial enrichment, efforts were made to match as much as possible the physicochemical properties of these decoys with the determined physicochemical properties for the active ligand set.[75] Therefore, eight molecular properties (number of heavy atoms, number of rings, net charge, log P, number of H-bond acceptors, number of H-bond donors, number of rotatable bonds, and polar surface area) were calculated using RDKit[76] for the active ligand set and for the initial decoy set. The initial decoy set was filtered based on the properties that showed the most notable differences when compared with the active ligand set properties (Figure S9A). In this way, decoys with log P values above or below the limits of log P distribution for the active ligands were excluded. The same procedure was carried out for the number of heavy atoms. Also, decoys with the net charge different from zero were eliminated. In this way, a filtered set of decoys was obtained, with more similar physicochemical properties with respect to the active ligands (Figure S9-B). This set of decoys, containing 1734 compounds (approximately 40 decoys per active ligand), was used for molecular docking screens.

Preparation of Ligand and Receptor Structures

To treat active ligands and decoys in the same way, they were first depicted in SMILES format. The active ligands that were not in the ChEMBL database were sketched from the closest ChEMBL entry. Then, Gypsum-DL 1.1.7[77] software was used to obtain the geometry-optimized 3D representations of the compounds. This software also allowed us to obtain different tautomeric, chiral, cis/trans isomeric, and ring conformational forms for the decoys and the active ligands. This is an important consideration to take into account in the structure-based drug discovery programs, given that in experimental assays, chemical compounds are often not tested in pure forms but rather in mixes containing all possible isomers.[78] In all the receptor homology models, the rotamers of Thr-6′ (see Figure ) were refined using UCSF Chimera.[79]

Molecular Docking

Molecular docking calculations were performed using AutoDock Vina 1.1.2,[80] 2Vinardo,[81] and LeDock v1.0,[82] using default parameters. Active ligands and decoys were considered flexible, whereas the receptors were held rigid. The search space for docking calculations in each of the homology models was set around the centroid of −2′, 2′, 6′, and 9′ residues of the five M2 helices (see Figure ), with size dimensions of 20 × 20 × 20 Å. It was checked that this box size contained the key residues of the binding site in all the RDL homology models. For molecular docking with AutoDock Vina, both receptor and ligand structures were prepared using AutoDock Tools,[83] by assigning them partial Gasteiger charges and combining the nonpolar hydrogen atoms. In the case of 2Vinardo, the structures were prepared using Open Babel.[84] For LeDock, the receptors were prepared using LePro, and the ligands were converted to the mol2 format using Open Babel.

Retrospective VS Performance Metrics

The enrichment of active ligands over decoys was analyzed for each of the RDL homology models by calculating the ROC curve with the R package “enrichvs”.[85] Then, the AUC was calculated. The AUC value is indicative of the global recovery of active ligands in the docking screen and takes values from 0 to 1. An ideal VS process, ranking all actives higher than decoys, will have an AUC value of 1, whereas if all decoys are ranked higher than actives, it will yield an AUC of 0. An AUC of 0.5 represents a VS process that randomly ranks actives and decoys. The BEDROC was also calculated using the “enrichvs” package. The BEDROC parameter emphasizes the early recognition of active ligands.[86] Its calculation involves an α coefficient that determines the weight given to the top-ranked compounds. In this work, two α coefficients were used: 20 and 100. With α = 20, 8% of the top-ranked molecules will account for 80% of the score, while with α = 100, the 80% of the BEDROC value will be determined by the top 1.6% of the ranking.[86] The BEDROC also ranges from 0 to 1, but the value it takes in the case of random ranking depends on a variety of factors, such as the α coefficient, the number of active ligands, and the number of decoys.[36] With the active and decoy sets used in this work, random ranking would yield a BEDROC score of 0.06 for α = 20 and a score of 0.03 for α = 100. The standard error of AUC and BEDROC scores were calculated using 2000 bootstrap resamplings of the VS rankings through the R package “boot”. Plots and other statistical calculations were performed in Python 3.7.6.

Three different RDL homology models were selected for MDS, according to their retrospective VS performance, each one in a different conformational state. The obtained docked structure of each fipronil/RDL complex was used for MDS using the GROMACS 2018.8 package with the all atom SLIPIDS force field[87] for the phospholipids and the AMBER99SB force field for the protein.[88] The construction of fipronil units to be used in MDS was performed using the protocol described before.[89] The optimized structure and restrained electrostatic potential charges were obtained as before using the Gaussian 03 package,[90] units were obtained using the GAFF force field,[91] and parameter files were translated to be used in GROMACS with ACPYPE.[92]

A POPC bilayer patch with 128 lipid molecules equilibrated at 293 K was obtained from SLIPIDS online resources (http://www.fos.su.se/~sasha/SLipids/Downloads.html) and replicated four times to create a 512-lipid bilayer patch. Water molecules were added to increase the size in the z axis for the patch to be able to accommodate the RDL model. The final bilayer system consisted of 512 POPC molecules with approximately 90,000 water molecules and a simulation box of 132 × 132 × 170 Å. This bilayer was equilibrated for 350 ns at 300 K in order to reach a surface tension value of zero.

Each of the three fipronil–RDL complexes was inserted into a pre-equilibrated POPC lipid bilayer using the InflateGRO script[93] after 26 iterations of scaling down by 0.95 and using a strong position-restraining force (100 kcal/mol/Å2) on the protein heavy atoms. To electrically neutralize the system, chloride ions were inserted randomly into the solvent, according to the charge of the model used. The total size was ∼470,000 atoms for each system. All the systems were initially minimized for 50,000 steps and then equilibrated for 100 ps in an NVT ensemble at 300 K, followed by 1 ns at NPT. In both cases, a restraint of 10 kcal/mol/Å2 fixing the protein αC was applied. Surface tension reached zero in all systems after equilibration. Finally, these equilibrated structures were used for MD production runs of 150 ns, in which the restraints on the protein αC were removed. The temperature was kept at 300 K using the Nose–Hoover thermostat,[94] and the Parrinello–Rahman barostat[95] was used to maintain the pressure at 1.01 bar. A simulation time step of 2 fs was set. All bonds were constrained using the LINCS algorithm. Cutoff distances were set to 1 nm for Lennard-Jones and electrostatic interactions. The particle-mesh Ewald method[96] was used for long-range electrostatic interactions.