ELIXIR-A: An Interactive Visualization Tool for Multi-Target Pharmacophore Refinement.

Pharmacophore modeling is an important step in computer-aided drug design for identifying interaction points between the receptor and ligand complex. Pharmacophore-based models can be used for de novo drug design, lead identification, and optimization in virtual screening as well as for multi-target drug design. There is a need to develop a user-friendly interface to filter the pharmacophore points resulting from multiple ligand conformations. Here, we present ELIXIR-A, a Python-based pharmacophore refinement tool, to help refine the pharmacophores between multiple ligands from multiple receptors. Furthermore, the output can be easily used in virtual pharmacophore-based screening platforms, thereby contributing to the development of drug discovery.

Introduction

Drug discovery and development is a complex, time-consuming, and expensive process. Computer-aided drug design approaches have the potential to accelerate this process cost-effectively when compared to the laborious traditional compound screening methods. Existing computational techniques include quantitative structure–activity relationship (QSAR), molecular docking-based high-throughput, and pharmacophore-based virtual compound screening. In the usual virtual drug discovery process, molecular docking, pharmacophore models, and 3D QSAR models are often used in combination.[1−3] Existing QSAR or docking approaches do not have the capability of pharmacophore refinement, and thus, there is a need for a technique for refining pharmacophores to identify the best set of pharmacophores for the modern drug discovery process.

The pharmacophore concept was introduced by Ehrlich in early 1909.[4] A pharmacophore is an ensemble of steric and electronic features that is necessary to ensure optimal supramolecular interactions with a specific biological target structure and to trigger (or to block) its biological response.[5] Pharmacophores describe specific ligand–receptor interactions as a generalized pattern. The first 3D pharmacophore screening software was developed by Gund in 1977.[6] Pharmacophore modeling approaches can be broadly classified as ligand-based or receptor-based.[7] Software programs such as LigandScout,[8] DISCO,[9] GASP,[10] GALAHAD,[11] HipHop,[12] HypoGen,[13] MOE (Chemical Computing Group, https://www.chemcomp.com), PharmaGist,[14] MolAlign,[15] and PHASE[16] have been developed to construct ligand-based pharmacophore models. Their performance differs based on the efficiency of the algorithm to handle ligand flexibility and alignment and requires a set of pharmacologically active ligands.[7] On the contrary, receptor-based methods, such as GBPM,[17] LigandScout,[8] Pharmit,[18−20] PyRod,[21] and ZINCPharmer,[22] analyze the receptor–ligand complex structures to isolate essential pharmacophoric features.[23] A 2D pharmacophore fingerprint is a form of a binary code that contains pharmacophore properties.[24−27] These pharmacophore fingerprints containing molecular fragments have been applied with multiple artificial intelligence-related models such as PTML,[28,29] Pharmacoprint,[27] and Pharm-IF.[30]

In situations where ligands are not known for the target receptor, methods such as CavityPlus,[31] GRID,[32] HS-Pharm,[33] Pocket version 4.0,[19,20] Shaper2,[34] GRAIL,[35] and SuperStar[36] have been developed to identify hotspots (highly probable ligand binding sites) on the receptor. Druggability simulations (molecular dynamics simulations conducted in the presence of drug-like organic molecules) assess ligand hotspots while maintaining receptor flexibility. Tools like SILCS-Pharm from the Mackerell lab[37,38] and Pharmmaker from the Bahar lab[39] have been developed to extract pharmacophore features from druggability simulations. However, there is a need for a systematic tool that analyzes and compares multiple pharmacophore models irrespective of their method of construction.

Here, we present the Enhanced Ligand Exploration and Interaction Recognition Algorithm (ELIXIR-A), an open-source, user-friendly application that serves the purpose of both pharmacophore modeling and pharmacophore mapping. ELIXIR-A is a Python-based application that can import pharmacophore models created in Visual Molecular Dynamics (VMD),[40] as well as manual coordinate input from any other platform. In addition, the output files from ELIXIR-A can be easily visualized in VMD[41] and can be exported to pharmacophore-based virtual screening platforms like Pharmit.[18]

Methodology

ELIXIR-A was developed to unify and simplify the interaction data from multiple pharmacophore models. This tool can accept two sets of pharmacophore models created directly by a pharmacophore generating platform. For example, the platform can accept ligand-based pharmacophore models created directly by Pharmit as the input. ELIXIR-A is developed on the modern 3D data processing library Open3D[42] to include the algorithms of fast global registration[43] and colored point cloud registration (colored ICP)[44,45] (Figure ). ELIXIR-A initially prepares each pharmacophore point represented as a point cloud with its radius, and the pharmacophores are color coded into different types. Then, ELIXIR-A aligns two structures and calculates the initial transformation matrix based on their geometric properties only. Finally, the two pharmacophore clouds are superpositioned with their geometric and pharmacophore type information, and the overlaid pharmacophore point is refined with matched points.

Figure 1

Flow diagram of the computational procedures of ELIXIR-A.

Pharmacophore Point Clouds

As a shape similarity approach, ELIXIR-A takes pharmacophore points and applies each of them to a three-dimensional volume with a point cloud. Each pharmacophore type consists of 1000 uniformly distributed points in a sphere. The radius of the pharmacophore cloud is defined in the occupancy category of the pdb file. The pharmacophore type is marked with three characters in the residue name category (Table S1). Throughout the complete analysis, the Biopandas framework is used to handle the point cloud data.[46]

Global Registration with RANSAC Iteration

The first alignment is the global alignment with features. Two pharmacophore point clouds are treated with the same process to calculate a vector of descriptors called the Fast Point Feature Histogram (FPFH).[47] FPFH is a 33-dimensional vector that describes the geometric characteristics and principles for a point. The FPFH vector can search points with similar features.[47] As ELIXIR-A attempts to align pharmacophore sites from multiple ligands or receptors, there may be distinct differences between some areas of the clouds considered as “noise”. It is known from experience that substantial noise levels can affect the alignment of point clouds with FPFH matching. To enhance the search of the FPFH matching algorithm, the random sample consensus (RANSAC) algorithm, proposed by Fischler and Bolles,[48] is included in the registration process. RANSAC can estimate parameters of a mathematical model from a set of observed data with “noise”. The global registration process calculates a preliminary rigid rotation and transformation matrix with their geometric characteristics. A fitness score will be given to evaluate this initial transformation.

Colored ICP

The second alignment is the local ICP alignment with pharmacophore features. The pharmacophore alignment utilizing colored features is inspired by the alignment of the red, green, blue model (RGB) image and its corresponding depth image by Park et al.[44] In general, the colored ICP algorithm extends the color vector as the extra dimension of the point coordinate data. The standard ICP algorithm will repeat the transformation of matrices to find the minimum square distance between two clouds. In comparison, the colored ICP algorithm will find the best fit of the extended paired matrices with the pharmacophore “color” information for each point. ELIXIR-A uses a slightly different number of iterations, fitness values, and root-mean-square deviation (RMSD) values from the default to optimize pharmacophore alignment. ELIXIR-A also supports the user adjustment of these critical parameters to obtain a more appropriate fit for the transformation.

Pharmacophore Refinement

After applying preliminary global transformation and local colored ICP transformation from the point cloud superposition to the first pharmacophore input, two pharmacophore points are partially overlapped. It is necessary to remove the nonoverlapped pharmacophore using the refinement algorithm. The algorithm first creates two 3′ N matrices A and B for two groups of pharmacophore datasets. A is defined as the transforming source group, and B is defined as the fixed target. The Euclidean distance for each point in group A to find the corresponding point in group B is calculated. Then, the Euclidean distance for each point in group B to find the corresponding point in group A is evaluated. In each group, if there is no corresponding point within the threshold distance, this point will be considered irrelevant and removed from the group. According to the definition of the pharmacophore, it is possible for different types of pharmacophores to exist in superposition; therefore, this refinement only calculates the geometric properties and does not consider the color parameters of the pharmacophore. The alignment of the two refined pharmacophore points will be visualized as the van der Waals surface (vdW) in VMD. These files will be saved in the same path as the two input pharmacophore files

Fitness Score Function

ELIXIR-A uses a fitness value to evaluate the effectiveness of the transformation for both alignment algorithms. Fitness is calculated using the formula below to find the volume ratio of overlap. The higher the fitness value, the better the alignment performance

Benchmark Compound Validation

To validate the efficiency of the pharmacophore refinement algorithm, a molecular dataset consisting of active inhibitors and inactive decoys targeting specific protein receptors was designed to reduce the testing bias. All the pharmacophore models were screened on the Pharmit platform, and the Directory of Useful Decoys (DUD-e) dataset was used as the small molecule library.[49] The benchmark receptors were human immunodeficiency virus type 1 protease (HIVPR), acetylcholinesterase (ACES), and cyclin-dependent kinase 2 (CDK2). These proteins are in different protein families (protease, esterase, and kinase). The feature “shared pharmacophores” in the LigandScout 4.4 Demo version was used as the comparison pharmacophore refinement algorithm in this study. Schrödinger Phase (release 2020–1) was used as another comparison software. The E-Pharmacophore plugin prepared the pharmacophore model under the receptor–ligand complex mode. The hypothesis alignment plugin was used for pharmacophore alignment with 2.0 Å feature matching tolerance and three minimum number of matching features. The overlapped pharmacophore was selected by manual inspection.

The molecular weight (MW), octanol–water partition coefficient (logP), and the total polar surface area (TPSA) were used to compare physicochemical properties between the screened active and decoy molecules. These physicochemical properties were calculated by open babel v2.4.0.[50] The properties were compared using the two-tailed Student’s t test.[51] The equality of two variances was given by the F-test.[52] The significance level was set to 0.05 (α = 0.05).

Enrichment Factor

Enrichment factors (EF) of the pharmacophore-based virtual screening represent the ability of the pharmacophore model to find true positive active inhibitors in the database in contrast to random selection. EF can be calculated using the ratio of the active inhibitors in the screened subgroup over the ratio of active inhibitors in the whole database. The higher the EF value, the better the pharmacology model performance

Results

Configuration

ELIXIR-A is designed as a GUI plugin based on VMD. The script itself is pre-configured. The first step is to install VMD and then edit the startup files in the VMD folder. The physical system requirements for ELXIR-A will not exceed the requirements of the VMD software itself. Once the setup files have been updated, the ELXIR-A folder must be placed in the VMD TCL plugin directory. The ELIXIR-A option can be found under the top menu. A Python environment is required to run the script for ELIXIR-A. The package was tested in cPython 3.8.8 with the following packages (Numpy==1.20.1; Open3D==0.13.0; Biopandas==0.2.9.dev0; and Matplotlib==3.3.4).[42,46,53,54]

Execution of Pharmacophore Alignment Jobs

ELIXIR-A provides two ways to introduce the pharmacophore information to the GUI (Figure ). The user can choose to pre-edit the text file with the “.pdb” extension via “Import pdb file:”, or directly use the Pharmit saved session file “.json”, which is the same as the pre-edited pdb file. The guidelines for the text files are described in Table S2.

Figure 2

ELIXIR-A interface. The target receptor can be loaded to visualize the pharmacophore–receptor interactions. According to different data sets, the default parameters can be optimized to improve the alignment performance.

Visualization of Results

The first pharmacophore cluster is recommended as the first primary site. The second cluster is interactively rotated and transformed to retrieve the conformation that best fits the primary site. Once all the data are prepared in the GUI, the “submit” button can be used to analyze the pharmacophore points’ alignment. This demonstration used ELIXIR-A to find refined structure-based pharmacophores from HIVPR (Figure ). These models were prepared in Pharmit and visualized in Maestro.[55] The receptor HIVPR is an essential enzyme of HIV replication. The active site of HIVPR was at the core of the dimerization interface (Figure ). The pharmacophores of the ligand(s) were extracted from Pharmit Engine (Figure ).[56] In the VMD OpenGL Display window, two groups of vdW spheres are represented as pharmacophores after the ELXIR-A alignment. The transformed pharmacophores were presented with solid spheres. The fixed pharmacophores were presented as transparent spheres. The ligand regarding the fixed pharmacophores is also shown in Figure left. The initial transformation has 15.39% of point cloud overlap, and the colored ICP improved the overlay by up to 65.66%. After refinement, two separate pharmacophore spheres in each group did not overlap and were therefore removed. A total of eight pharmacophores shared the superposition space, and this result can be used for further drug discovery. In addition to the HIVPR, the use of similarity refinement for the analysis of two other types of proteins is also applied in Table .

Figure 3

Pharmacological description of the two HIVPR in the crystal structure in complex with inhibitors.

Figure 4

Structure-based pharmacophore refinement between HIVPR-inhibitor complexes utilizing ELIXIR-A. The transformed pharmacophores are presented with solid spheres. The fixed pharmacophores are presented using transparent spheres.

Table 1

ELIXIR-A Aligned Pharmacophore Cluster Selection

						best fitness score
PDB IDa	target	ligand IDb	initial points	refined points	RMSD (Å)	global registration	colored ICP	ref
1YT9	HIVPR	OIS	11	8	0.267	15.39%	65.66%	(57)
2FDE		385	9	8				(58)
1Q84	ACES	TZ4	14	5	0.256	14.65%	42.15%	(59)
2CEK		N8T	7	3				(60)
6INL	CDK2	AJR	7	4	0.287	30.612%	31.82%	(61)
4EOR		4SP	6	4				(62)

The protein complexes were downloaded from the Protein Data Bank (https://www.rcsb.org).

These ligands are the local substrates corresponding to the protein complexes. Detailed information on the ligands is available at https://www.rcsb.org/ligand/(ligand ID).

Benchmark Compound Validation

A benchmark compound study was used to validate the effectiveness of the module in generating enriched screened molecules using refined pharmacophores. For comparison, the same analysis was carried out for the benchmark compounds using the two software packages that possess pharmacophore alignment modules, LigandScout and Phase, which were accessible to the researchers. We used three targets for benchmarking after considering 11 data sets that are freely available for validation in the Zinc15 database, which also would give positive results for all three screening platforms, i.e., ELIXIR-A, Phase, and LigandScout. The usage of three targets is in line with other works for testing screening algorithms.[38,63−65]

ELIXIR-A refined pharmacophores were extracted from the overlaid points (Figure , right). Five features with three different types were built from 20 initial features. The HIVPR pharmacophore model prepared directly by Pharmit could not generate any screened results because the initial model contained too many pharmacophore features, confirming the need for a tool to refine pharmacophores for effective virtual screening. To compare the performance of ELIXIR-A in comparison to other software capable of pharmacophore screening, we ran the screening algorithms using ∼25,000 starting molecules for three distinct receptors.

For HIVPR, the ELIXIR-A algorithm was able to output 19 molecules, with three being active (3/536) with an enrichment factor of 10.69 (Table ). LigandScout was able to retrieve ∼47% active compounds (253/536), while the enrichment factor was only 4.91. Phase found three aromatic ring features based on the alignment of two structure-based complexes. The enrichment factor of 0.45 suggested that this model did not match the true binding modes. HIVPR has four key hydrophobic features located in the drug pocket, which are important for binding.[38,66] ELIXIR-A and LigandScout both found three of the four hydrophobic features for virtual screening (Figure ). Thus, high enrichment factors were obtained using ELIXIR-A and LigandScout for the HIVPR target.

Table 2

Enrichment Analysis for Refined Pharmacophore Models

		active and inactive compounds in the dataset		virtual screening resultsb
approach	target	active	decoy	active	decoy	EF
ELIXIR-A	HIVPR	536	35,750	3	16	10.69
LigandScout				253	3236	4.91
Phase				36	5409	0.45
ELIXIR-A	ACES	408	26,250	201	3349	3.70
LigandScout				190	2913	4.00
Phase				167	8308	1.29
ELIXIR-A	CDK2	335	27,850	189	1597	8.90
LigandScouta
Phase				33	345	7.34

LigandScout lacked sufficient independent matching of pharmacophore features within tolerances to generate alignment.

Some molecules in the dataset have multiple tautomer confirmations. The output will only include a unique active/decoy molecule with minimum RMSD from pharmacophore matching.

Figure 5

Refined pharmacophore models of HIVPR generated by ELIXIR-A (A), LigandScout (B), and Phase (C).

For ACES, ELIXIR-A was able to retrieve ∼50% (201/408) active compounds with an enrichment factor of 3.70 (Table ). LigandScout was able to retrieve ∼47% active molecules with an enrichment factor of 4.00. Phase retrieved ∼40% active molecules with an enrichment factor of 1.29 due to the high number of screened decoys (8038/26250). The ACES active binding pocket was long, narrow, and hydrophobic.[60] Conserved aromatic residues were found at the peripheral site.[67] All three algorithms in this study modeled one hydrophobic feature at the active site and two hydrophobic or aromatic features at the peripheral site (Figure ). Comparatively, ELIXIR-A was able to match the highest number of active compounds in the dataset.

Figure 6

Refined pharmacophore models of ACES generated by ELIXIR-A (A), LigandScout (B), and Phase (C).

For CDK2, ELIXIR-A was able to retrieve ∼55% (189/335) active compounds with an enrichment factor of 8.90 (Table ). Phase retrieved ∼10% (33/335) active molecules with an enrichment factor of 7.34. Multiple CDK2 pharmacophores were reported with different pharmacophore features.[68] Meanwhile, most of them include one hydrogen bond acceptor, one hydrogen bond donor, and two hydrophobic features. The ELIXIR-A model included all these common features and therefore matched half of the active compounds in the dataset (Figure ). The Phase model recognizes the nonpolar features as the aromatic rings, which matched 10% of the active molecules. Here, ELIXIR-A was able to retrieve a much higher number of active compounds as compared to Phase, while LigandScout failed to retrieve any.

Figure 7

Refined pharmacophore models of CDK2 generated by ELIXIR-A (A) and Phase (B).

This benchmark compound validation indicated that ELIXIR-A could be a valuable tool to refine pharmacophores for better enrichment during virtual screening.

The physicochemical properties of ligands screened by various pharmacophore refinement models were compared (Table S3). It is noted that at least one of the physicochemical properties (MW, TPSA, or logP) differed significantly between the active and decoy molecules in the seven models. In the ACES-ELIXIR-A model, MW, TPSA, and logP showed significant differences. In the CDK2-ELIXIR-A model, MW and logP showed a significant difference. These results suggested that ELIXIR-A was able to satisfactorily differentiate between active and decoy molecules in terms of key physicochemical aspects. It is noted that further physicochemical studies such as QSAR could be utilized to further enhance the ELIXIR-A pharmacophore refinement platform.

Discussion

Pharmacophores are a set of steric and electronic features that recognize optimal supramolecular interactions. Structure-based modeling using protein–ligand interactions and ligand-based modeling using common chemical features from a set of active/inactive ligands are the two common approaches for building 3D pharmacophore models.[69,70] ELIXIR-A is not typical pharmacophore modeling software; rather, it is a pharmacophore refinement algorithm that uses sets of pharmacophores as input and aligns these pharmacophores in 3D space to identify any overlap. ELIXIR-A uses a computer vision-inspired ICP variant algorithm to align multiple pharmacophore models with similar geometric and physicochemical properties as point clouds into a refined point cloud. Although researchers have applied this ICP algorithm for drug discovery to match three-dimensional protein structures[71] as well as the alignment of protein binding cavities,[72,73] it has not been used for pharmacophore refinement. ELIXIR-A fills this need.

ELIXIR-A uses a geometry-based approach to find pharmacophore similarity between two protein–ligand binding pockets before calculating the binding energy between the receptor and small ligand molecules using colored ICP. ICP requires a good initial transformation to ensure that the point cloud converges to a minimum acceptable value.[74] Also, the presence of outliers (nonuniform points) can affect the alignment in ICP. ELIXIR-A utilizes FPFH matching and RANSAC iterations to solve the global fitting and outlier problems. Another algorithm that can perform similar calculations is the Kabsch algorithm. The Kabsch algorithm is widely used in bioinformatics and can calculate the RMSD between two 3D protein or pharmacophore structures via rotation and translation.[75,76] When the point correspondences are known, the Kabsch algorithm can be applied. Since point correspondences are not known, the Kabsch algorithm cannot be applied here and thus ICP was used. When using colored ICP, the whole cloud is paired to figure out the corresponding point, and the points in one group can be mapped with multiple corresponding points in the other group. It is also possible to find some important unpaired pharmacophore points when comparing the output results. Another algorithm that is commonly used for ligand-based pharmacophore generation is linear assignment for molecular dataset alignment (LAMDA), which finds globally optimal pairing between objects.[70] LAMDA was not considered for ELIXIR-A due to its limitation in ligand-based pharmacophore generation. Genetic algorithms (GA) were used to calculate the initial transformation position before ICP transformation.[77−80] For example, GA is used for pharmacophore generation via molecular docking-based data sets such as AutoDock Vina.[81] Meanwhile, ELIXIR-A used FPFH matching with RANSAC iteration for global alignment. This is mainly because of the good compatibility between Global Registration and color ICP supported by Open3D. Although ELIXIR-A currently cannot use GA, when the molecular structures of the compared proteins are very different, such a powerful algorithm would be helpful to optimize refinement results.

For the initial pharmacophore analysis, ELIXIR-A and LigandScout used the same method to find the rigid superposition of two pharmacophore clusters, while ELIXIR-A used the ICP variant algorithm and LigandScout used the Hungarian algorithm. In general, the time complexity of the ICP algorithm in ELIXIR-A was O(N2), and the time complexity of the Hungarian algorithm ranged from O(N4) to O(N3). Thus, ELIXIR-A provided better time complexity with the potential to efficiently solve large pharmacophore alignment problems.

Conclusions

ELIXIR-A is a python-based VMD plugin used to help refine pharmacophore models in situations where a large number of pharmacophores are present from multiple models due to various ligand–receptor interactions and many conformations of a ligand on a receptor binding site. The ELIXIR-A GUI can help refine the pharmacophores generated from multiple ligands from multiple receptors using the ICP variant algorithm. The output from ELIXIR-A can be used in virtual pharmacophore-based platforms for compound screening.

Requirements and Availability

The ELIXIR-A package was developed in python3 for the data operations and tcl/tk for the user interface. The package was tested on the Ubuntu 20.04, macOS Big Sur, and Windows 10 systems. The python versions were 3.7.4 or later. The latest version of ELXIR-A is available for download on GitHub (https://github.com/sfernando-BAEN/ELIXIR-A/releases). VMD was downloaded from their official websites. A 64-bit version of VMD is recommended for this package. ELIXIR-A follows the Apache-2.0 license and is open-source on GitHub.