Multiple conformational states in retrospective virtual screening - homology models vs. crystal structures: beta-2 adrenergic receptor case study.

BACKGROUND: Distinguishing active from inactive compounds is one of the crucial problems of molecular docking, especially in the context of virtual screening experiments. The randomization of poses and the natural flexibility of the protein make this discrimination even harder. Some of the recent approaches to post-docking analysis use an ensemble of receptor models to mimic this naturally occurring conformational diversity. However, the optimal number of receptor conformations is yet to be determined. In this study, we compare the results of a retrospective screening of beta-2 adrenergic receptor ligands performed on both the ensemble of receptor conformations extracted from ten available crystal structures and an equal number of homology models. Additional analysis was also performed for homology models with up to 20 receptor conformations considered.
RESULTS: The docking results were encoded into the Structural Interaction Fingerprints and were automatically analyzed by support vector machine. The use of homology models in such virtual screening application was proved to be superior in comparison to crystal structures. Additionally, increasing the number of receptor conformational states led to enhanced effectiveness of active vs. inactive compounds discrimination.
CONCLUSIONS: For virtual screening purposes, the use of homology models was found to be most beneficial, even in the presence of crystallographic data regarding the conformational space of the receptor. The results also showed that increasing the number of receptors considered improves the effectiveness of identifying active compounds by machine learning methods. Graphical abstractComparison of machine learning results obtained for various number of beta-2 AR homology models and crystal structures.

Background

G protein-coupled receptors (GPCRs) constitute a large superfamily of signaling proteins that share a common topology of 7 transmembrane (7TM) helices and transduce signals across the cell membrane. Because GPCRs are responsible for most of a cell’s communication with its environment, their malfunctions are associated with various disease states, mainly those related to the central nervous system (CNS). For this reason, GPCRs are a very important target base for drugs [1-3].

The beta-2 adrenergic (B2AR) receptor, the subject of this case study, is representative of the class A GPCRs and is involved in mediating the relaxation of smooth muscle, glycogenolysis and glucogenesis in the liver and regulation of the metabolism of cells in skeletal muscle. Β2AR is also responsible for increased cardiac output, facilitation of the release of neurotransmitters, and regulation of various other physiological processes [4-7]. B2AR is also one of the most studied 7TM structures; it was first crystalized in 2007 [8], and as of November 2014, 16 crystals with a variety of structurally and functionally unique ligands are available via the Protein Data Bank (PDB), making this receptor a strong base for in silico structural studies.

Our previous study applying Machine Learning (ML) to post-docking analysis used Structural Interaction Fingerprint (SIFt) profiles created upon three different crystalline conformations of receptors [9,10]. It showed the applicability of this approach to ligand-protein complexes evaluation for Virtual Screening (VS). In addition to the issue of the applicability of crystal structures in VS, this study also investigates the influence of the number of conformations used in per-ligand interaction profiles, for both crystal structures and homology models, on retrospective screening performance. The VS setup consisted of four groups of compounds: active, inactive, DUD (Directory of Useful Decoys) decoys, and random ZINC subsets; three sets of experiments were prepared to discriminate between active and inactive compounds from each of the decoy collections.

The support vector machine (SVM) was the classification algorithm chosen and the VS performance was measured with the Matthews Correlation Coefficient (MCC).

Results and discussion

Crystal structures vs. homology models

Because the number of crystal templates used for homology models construction would affect the clarity of the presented results, the comparison of homology models and crystal structures is shown for the templates providing the best (M2R) and the worst (D3R) results (in terms of the discrimination between actives and true inactives) – Figure 1; the outcomes for the remaining templates are available in the Additional files section (Additional file 1: Figure S1). The use of vast numbers of templates for homology modeling follows the protocols used in previously published data and ensures maximum VS performance [11]. Due to a limited number of available crystal structures, the maximum number of receptor conformations in this comparison is restricted to 10 (starting from 3).

Figure 1

Comparison of MCC values obtained in the ML-based experiments of docking results to homology models built on M2R and D3R template and crystal structures for discrimination between a) actives/true inactives, b) actives/DUDs, and c) actives/ZINC. The figure presents the MCC values obtained for homology models of beta-2 adrenergic receptor (the best and the worst template) and for crystal structures of this receptor in experiments distinguishing the following class of compounds: (a) actives/true inactives, (b) actives/DUDs and (c) actives/ZINC.

The results of retrospective VS (Figure 1) show that homology model-based screening significantly outperforms experiments conducted for the collection of crystal structures, with MCC improvement of 0.4 for the best set of conformations. In addition, all types of classification experiments (actives/true inactives, actives/DUDs, and actives/ZINC cmds) confirm this dependency. The MCC spreads for different templates were of little significance: variation between the best and the worst performing template ranged from 0.1 for actives/true inactives experiments to less than 0.05 for the other two VS scenarios.

For homology models, MCC values obtained for actives/true inactives discrimination were the lowest (~0.5 – 0.55). However, for actives/DUDs and actives/ZINC cmds classifications, MCC values exceeded 0.8, with a slight preference towards actives/ZINC experiments.

On the other hand, studies performed for crystal structures resulted in MCC of 0.2 for actives/true inactives (this best result was obtained for the SIFt profile composed of 8 receptor conformations), 0.47 for actives/DUDs (6 conformations) and 0.55 for actives/ZINC (9 conformations).

The obtained results show that the conformational flexibility provided by homology models allows for better accommodation of diverse ligands and therefore better screening performance in this interaction-centric type of experiments. Because crystal structures are limited in terms of chemical space of co-crystalized ligands, they are not yet able to provide a sufficient conformational landscape for efficient identification of active compounds.

Influence of the number of considered conformations on screening performance for homology models

Due to the substantial amount of data, a detailed analysis was conducted only for the best performing set of SVM parameters, in terms of MCC value, and also for the best and the worst template only (M2R and D3R, respectively); full datasets are included within the Additional files (Additional file 2: Figure S2). The MCC values related to different numbers of considered conformations are presented in Figure 2.

Figure 2

MCC values obtained for various numbers of models included in the SIFt profile for models built on M2R and D3R templates for a) actives/true inactives, b) actives/DUDs, and c) actives/ZINC discrimination. The figure presents the MCC values obtained for various numbers of models included in the SIFt profile for homology models constructed on M2R and D3R templates in the form of the heat map for a) actives/true inactives, b) actives/DUDs, c) actives/ZINC cmds discrimination. Red frames indicate the step when the best model in terms of the AUROC value was added to the profile.

In addition, differential graphs illustrating MCC change after including subsequent models (adding one-by-one- from 3 to 20 forming at the end 20-models-based profile) were prepared (Figure 3, Additional file 3: Figure S3). In each case, the addition of the model that was characterized by the highest area under the ROC curve (AUROC) at the stage of models evaluation was highlighted – in Figure 2 by red frame, in Figure 3 by the application of brighter colour.

Figure 3

Difference in MCC caused by the inclusion of additional receptors in the profile for a) actives/true inactives, b) actives/DUDs, c) actives/ZINC cmds discrimination. The figure presents the changes in MCC obtained after the inclusion of additional receptors in the SIFt profile for homology models constructed on M2R and D3R templates. The introduction of the best model in terms of the AUROC values is indicated by lighter colour.

The general outcome emerging from the results (Figure 2; Figure 3) aligns with the results obtained for the comparison of crystal structures with homology models. An increased number of conformations included in a SIFt profile leads to an improvement of VS performance, however, for some isolated cases, the contribution of subsequent models may be negative. The number of models providing the highest MCC was 20 in the majority of cases (as shown in Table 1), and the worst performing set of conformations was 3 for all but two sets of models. The improvement of MCC was not linear; however, its values were noticeably lower for a low number of conformations considered.

Table 1

The optimal and the worst number of models included in the SIFt profiles in terms of classification effectiveness

Template	Actives vs. true inactives		Actives vs. DUD		Actives vs. ZINC
	Optimal number of receptors	The number of receptors with min MCC	Optimal number of receptors	The number of receptors with min MCC	Optimal number of receptors	The number of receptors with min MCC
5-HT1BR	20	4	18	4	8	3
5-HT2BR	18	3	18	3	19	3
A2AR	19	3	20	3	20	3
Beta1R	17	10	20	6	20	3
CXCR4R	10	4	15	3	19	6
D3R	20	3	20	3	15	3
H1R	20	5	20	3	20	4
M2R	20	3	20	3	20	3
M3R	20	3	20	3	18	3
Crystal structures	8	3	6	3	9	3

MCC fluctuations occurring in the actives/true inactives classification stage of the experiment were the highest out of all three scenarios. This scenario also had several situations where additional considered conformations lowered the screening performance. On the other hand, filtering actives against DUD and random ZINC decoys led to a clear dependency between MCC and the number of models: the higher the number of models included in the profile, the higher the MCC values.

The impact of receptors bearing the highest AUROC values during the model selection step (conformation 4 and 19 for D3R and M2R templates, respectively) proves that the performance of individual homology models has little influence on the obtained results and, in some cases (conformation 19 based on M2R in screening against ZINC subset – Figure 2c), can even lower VS performance.

Although the MCC changes induced by including new conformations to the ligand profiles seem to be negligible, the cumulative effect for VS experiments leads to a significant improvement of screening performance by up to 20% (Additional file 4: Figure S4). Interestingly, the absolute values of MCC difference oscillated at approximately 0.1, regardless of the scheme of the experiment.

Influence of the number of considered conformations on screening performance for crystal structures

The results show that the experiments using multiple crystalline conformations are significantly more prone to screening performance fluctuations (Figure 4). The amplitude of these fluctuations ranges from a 0.3 improvement to a 0.4 decrease in terms of MCC value. This variation of the results obtained for crystal structures is connected with the specificity of the individual crystallized with the proteins. The drop in MCC, observed after adding the last two conformations (3P0G [12] and 3PDS [13]), is a consequence of the crystals being agonist bound. The conformation of an activated GPCR results in limited ligand accessible volume, significantly reducing the quality of docking results. Following this lead, the influence of activation state of the crystal template on classification efficiency was also examined for beta-2 homology models constructed on activated and deactivated M2R template. The direction of changes when adding subsequent SIFts to the profile was preserved but the results were slightly (~2-3%) better for models constructed on deactivated template (Additional file 5: Figure S5).

Figure 4

Changes in MCC after the inclusion of additional crystals in the profile. The figure presents the changes in MCC values caused by the addition of subsequent crystals to the profile (adding one-by-one- from 3 to 10 forming at the end 10-crystals-based profile) in experiments where docking was performed to crystal structures of beta-2 adrenergic receptor.

Conclusions

In this study we compared the performance of the collection of crystal structures with corresponding sets of homology models in retrospective VS experiments designed to consider multiple conformations of a target receptor. The results demonstrated that the bundle of homology models significantly outperformed the crystalline-based approach in terms of MCC, which agrees with results from previous reports [14]. The main reason behind this difference in screening effectiveness is the limited conformational space of the crystal structures, which is a consequence of adaptation to the co-crystalized ligands, thus biasing the conformation of the complex. Shallow conformational landscapes of the crystal structures of the receptors are also caused by low structural diversity of the crystalline ligands, limiting the possible spatial orientations of residues.

The second component of this research investigated the effect of increasing the number of considered conformations. The conclusion emerging from all schemes of experiments (screening active compounds against truly inactive, DUD, and ZINC decoys) is that high coverage of the conformational space of the receptor models leads to more effective screening. A probable reason behind this observation is that the inclusion of more conformations into a docking protocol neglects the fluctuations of docking poses and provides a more coherent binding mode for a given ligand, therefore enabling a clearer discrimination between active and inactive compounds. Extending the population of conformations would most likely increase the MCC up to the limit defined by the number of compounds that were not docked into any receptor model, yet the increasing computational cost of such tests may render the results not worth the effort. Although there is no actual boundary for the number of conformations to include, the results shown here prove that three models/crystals leave sufficient space to improve VS performance.

Methods

To maintain the coherence of the ligand data, the compounds of known activity (divided into sets ‘actives’ and ‘true inactives’) were extracted using a strict protocol. All structures with verified activity towards the B2AR were selected from the ChEMBL database [15]. Only those compounds whose activity was quantified in Ki or IC50 (with the assumption that Ki = IC50/2) and that were tested on human cloned, rat cloned, or native receptors were taken into account. A compound was considered active when the Ki value assigned to it was lower than 100 nM, and the compound was considered inactive when this activity parameter was higher than 1000 nM. The compounds were clustered with Canvas [16], and the number of clusters was set to approximately 30% of the total number of compounds from a particular group. Cluster centroids were used for the primary evaluation of homology models. In addition, two sets of decoy compounds were generated, one following the DUD methodology [17] and one random subset of the ZINC database [18]. Both of the decoy collections contained 2000 compounds; DUDs were randomly picked to narrow the count of the set.

The homology models of the B2AR were constructed. Nine crystal templates were used for this purpose: serotonin receptors 5-HT1B and 5-HT2B, adenosine receptor A2A, adrenergic receptor beta-1, chemokine receptor CXCR4, dopamine receptor D3, histamine receptor H1, and muscarinic receptors M2 and M3 (Table 2). The sequence alignment was performed manually and only for the transmembrane helices. Loops were not modelled. For each template, 20 models were generated with Modeller 9v13 software [19] and were evaluated by AUROC in the docking of cluster centroids from actives and true inactives sets (Table 3). Additional homology models were also prepared in the same manner for deactivated M2R structure (4MQS).

Table 2

Crystal structures used as templates for homology modeling of beta-2 adrenergic receptor

Template	PDB ID	Resolution [Å]
5-HT1BR	4IAR [23]	2.70
5-HT2BR	4IB4 [24]	2.70
A2AR	3QAK [25]	2.71
Beta1R	2Y00 [26]	2.50
CXCR4R	3OE0 [27]	2.90
D3R	3PBL [28]	2.89
H1R	3RZE [29]	3.10
M2R	3UON [30]	3.00
M3R	4DAJ [31]	3.40

Table 3

Number of compounds used for the preevaluation of homology models

Group of compounds	All	Number of clusters	Number of centroids after Ligprep
Actives	271	81	103
Inactives	324	97	173

The three dimensional structures of the compounds, along with protonation states and atom types, were assigned with LigPrep software [20]. For some compounds, several protonation states were generated what increased the initial number of instances. The docking was performed with GLIDE 5.0, with the number of output poses limited to one for both homology models (Table 4), and a collection of B2AR crystal structures (Table 5).

Table 4

Compound counts for retrospective screening scenarios

Group of compounds	Total number of compounds	Number of compounds after Ligprep
actives	271	550
Inactives	324	601
DUDs	2000	2526
ZINC	2000	2557

Table 5

Crystal structures of beta-2 adrenergic receptor used in the study

PDB ID	Resolution [Å]
2RH1	2.40
3D4S	2.80
3NY8	2.84
3NY9	2.84
3NYA	3.16
3KJ6	3.40
2R4R	3.40
2R4S	3.40
3P0G	3.50
3PDS	3.50

After this initial models evaluation, all compounds from a particular group of molecules (actives, true inactives, DUDs, and ZINC) were docked into the constructed homology models and crystal structures. Ligand-receptor complexes received from the docking procedure were represented by the Structural Interaction Fingerprint [9] which have a type of a binary string that describes the interaction of a ligand with each of the amino acids of the protein; the string is divided into nine-bit chunks that refer to particular amino acid residues. The type of interactions that are taken into account include the presence of any interaction, an interaction with the main chain, an interaction with a side chain, a polar interaction, a hydrophobic interaction, a hydrogen bond acceptor, a hydrogen bond donor, an aromatic interaction, and a charged bond.

For each compound that had at least one pose in a population of receptor conformations, the SIFt profile was calculated. On each position in the string, the values were averaged over all models/crystal structures for the given conformational landscape considered (per ligand SIFt profile). The number of receptor conformations used in the experiments ranged from 3 to 20 (10 for the crystal structures).

The per ligand SIFt profiles were input for machine learning experiments conducted with the use of the WEKA package [21]. The task of the ML algorithm was to distinguish active from inactive or decoy compounds. Support vector machines algorithm [22] was used as a classification method with linear function as a kernel. This model was developed by Vapnik [22] with a core concept of seeking the hyperplane separating the binary-labeled data with the maximum possible margin. This can be written as the following optimization problem:with w being the normal vector to the hyperplane and yi being the class to which the particular example is assigned (in case of binary labeled data, y ∈ {− 1, + 1}). C is the parameter that controls the tradeoff between the correct classification and large margin.

However, in real applications, the data are not usually linearly separable, and the application of the kernel trick is required. In our paper, the linear kernel ( for d-dimensional feature space) was applied.

The original optimization problem is transformed to the dual form with the use of Lagrange’s multipliers α:α represents weights that are assigned to particular example x from the training data, and in the dual form, C constitutes the upper bound of α values.

The optimization of C values was performed (the following C values were checked: 0.01; 0.1; 1; 10; 100; 1 000; 10 000). The experiments were carried out in a 10-fold cross-validation mode.

The effectiveness of machine learning methods was measured with MCC being a balanced measure for such kind of experiments and expressed by the following formula:where:

TP – number of true positives,

TN – number of true negatives,

FP – number of false positives

FN – number of false negatives