A Quantum-Inspired Method for Three-Dimensional Ligand-Based Virtual Screening
Maritza Hernandez
[
[email protected]
Guo Liang Gan
[
Kirby Linvill
Carl Dukatz
[
Jun Feng
Govinda Bhisetti
[
Abstract
Measuring similarity between molecules is an important part of virtual screening (VS) experiments deployed during the early stages of drug discovery. Most widely used methods for evaluating the similarity of molecules use molecular fingerprints to encode structural information. While similarity methods using fingerprint encodings are efficient, they do not consider all the relevant aspects of molecular structure. In this paper, we describe a quantum-inspired graph-based molecular similarity (GMS) method for ligand-based VS. The GMS method is formulated as a quadratic unconstrained binary optimization problem that can be solved using a quantum annealer, providing the opportunity to take advantage of this nascent and potentially groundbreaking technology.
In this study, we consider various features relevant to ligand-based VS, such as pharmacophore features and three-dimensional atomic coordinates, and include them in the GMS method. We evaluate this approach on various datasets from the DUD_LIB_VS_1.0 library.
Our results show that using three-dimensional atomic coordinates as features for comparison yields higher early enrichment values. In addition, we evaluate the performance of the GMS method against conventional fingerprint approaches. The results demonstrate that the GMS method outperforms fingerprint methods for most of the datasets, presenting a new alternative in ligand-based VS with the potential for future enhancement.
1QBit]1QB Information Technologies (1QBit), 458-550 Burrard Street, Vancouver, BC, V6C 2B5, Canada
1QBit]1QB Information Technologies (1QBit), 458-550 Burrard Street, Vancouver, BC, V6C 2B5, Canada
\alsoaffiliationSchool of Computing Science, Simon Fraser University, 8888 University Drive,
Vancouver, BC, V5A 1S6, Canada
Accenture Labs]Accenture Labs, Accenture PLC, 415 Mission Street, Suite 3300, San Francisco, CA, 94105, USA
Biogen]Biogen, 225 Binney Street, Cambridge, MA, 02142, USA
1 Introduction
The continued need for the development of innovative new medicines faces major challenges due to the increasing costs of drug development and the high failure rate of potential drug candidates. For example, only 19% of new drugs that enter clinical trials eventually get FDA approval 1. Thus, there remains high demand for innovation and opportunities to improve the drug discovery process. In recent years, the VS of small molecules has become a routine and integral part of the drug discovery process.
In silico VS is a critical, early phase drug discovery approach that enables the identification of potential candidate molecules from a large molecular database, in a high-throughput manner. A number of high-throughput screening (HTS) methods have been developed that offer alternative or complementary strategies for VS. They can be divided into two categories: structure based and ligand based. Structure-based approaches, such as docking algorithms, are based on 3D structural information of the protein target, and a scoring function is used to measure how well a ligand binds to the active site 2. In situations where structural information of the binding site is unknown, ligand-based approaches such as similarity searching and pharmacophore mapping are utilized for VS 3. The conventional molecular representation used for ligand-based VS methods is based on 2D fingerprint representations. A fingerprint representation is a binary vector in which each entry indicates the presence or absence of a substructure in a given molecule. In general, they are computationally inexpensive and simple to use; however, they do not consider all of the relevant molecular features. One of the most widely used fingerprint representations is extended-connectivity fingerprints (ECFP) 4. ECFPs are circular fingerprints constructed by encoding an atom’s neighbourhood while expanding iteratively into adjacent neighbourhoods until a given number of iterations is reached.
Alternative molecular representations are intended to characterize the 3D nature of molecules 5. These methods are based either on molecular shape or pharmacophore shape alignment. Molecular shape methods aim to measure similarity based on the maximum degree of overlap of molecular shapes (or volumes). Examples of algorithms that use a shape-based approach are OpenEye Scientific’s rapid overlay of chemical structures (ROCS) 6, 7 and Schrödinger’s Phase 8, 9. Graph matching algorithms have also been considered for the molecular similarity problem. In a molecular graph representation, atoms and bonds are represented by nodes and edges, respectively. In similarity searching applications, graph matching algorithms are commonly based on optimization problems, for instance, the maximum common subgraph (MCS) problem 10 or the optimal assignment problem 11. Three-dimensional approaches take into account conformational properties; hence, they are more computationally expensive than 2D approaches. More recently, there has been a growing interest in applying machine learning methods to ligand-based VS. One example is the molecular graph convolutions 12 machine learning architecture. This method represents molecules as graphs in a deep learning system. Although it does not outperform all fingerprint-based methods, the authors introduce an alternative ligand-based VS method that is still being explored.
Hernandez et al. 13 proposed a graph-based molecular similarity approach that can be implemented with the use of a quantum annealer 14. The algorithm finds the maximum weighted common subgraph (MWCS) of two molecular graphs. The main difference between their algorithm and previous methods is that it finds the MWCS by solving the maximum weighted co-k-plex problem of an induced graph. This method has been reported to yield improvements in accuracy over conventional fingerprint methods when predicting mutagenicity in small molecules. The maximum co-k-plex problem is, in general, NP-hard, which means the time to solve the problem increases exponentially with the number of variables 15. Quantum annealing is a promising approach, with the potential to harness quantum mechanical effects, to solve hard optimization problems: it may be able to address the co-k-plex problem and, consequently, the molecular similarity problem more effectively than classical approaches. A study assessing the performance of a quantum annealer in solving the molecular similarity problem was performed by Hernandez and Aramon 16, and provides useful insights into new techniques for using the quantum annealer and addressing some of its hardware limitations.
The purpose of this study is to extend the applicability of the GMS approach introduced by Hernandez et al. 13 to ligand-based VS experiments. In this paper, we study how various molecular features, including 3D coordinates, affect similarity score and, therefore, the performance of the VS experiments. The GMS methodology described in this paper has been implemented in the 1QBit SDK 17. The performance of VS, using the GMS method, is compared against conventional fingerprint methods. Previous studies comparing various methods for VS have, in general, concluded that ligand-based 2D methods tend to yield better performance than docking methods, and that 2D fingerprint approaches generally outperform 3D-shape-based methods 18, 19. Remarkably, the 3D approach implemented in the GMS method generally outperforms its 2D counterpart for the 13 target proteins used in this study.
2 Similarity Methods
Similarity is a fundamental concept in the field of chemoinformatics. Developing methods to evaluate a measure of similarity is, however, a difficult task due to the nature of the concept of similarity: “like beauty, it is in the eye of the beholder” 20. The criteria used to define similarity can vary accross applications. For instance, organic chemists may be interested in classifying molecules in terms of the core scaffold fragment and its substructures, whereas physical chemists may focus on physicochemical properties, such as excluded volume and electrostatic properties. The similarity method introduced by Hernandez et al. 13 compares molecules according to chemical descriptors. In this work, we have modified the criteria for considering two molecules as similar by including weighted pharmacophore features. Additionally, the new similarity criteria allow the inclusion of partial matches. In this section, we describe how features and their relevance are incorporated in the GMS method.
The Graph-Based Molecular Similarity Method
The GMS method consists of three steps: 1) modelling molecules as graphs; 2) solving the co-k-plex problem; and 3) calculating a coefficient that measures the similarity between the two input molecules. We describe these three steps in an earlier paper 13. In the sections that follow, we present an overview of these processes with a special focus on the variations implemented in this study. A general scheme of this method is shown in Fig. 1.
Step 1: Modelling Molecules as Graphs
The first step of the GMS method is to model molecules as graphs. Individual atoms and ring structures in the molecule are represented by individual vertices in the graph. Here, atoms connected cyclically are referred to as ring structures, without differentiating whether they are aromatic. Bonds connecting a pair of atoms or rings are represented by edges that connect the respective pairs of vertices. If two vertices representing rings share one or more atoms, an edge between these vertices is then added and labelled an artificial bond, emphasizing that it is not intended to represent a natural chemical bond. An illustration of a molecular graph representation is shown in Fig. 2.
Formally, let G=(V,E,LV,LE) be a labelled graph representing a molecule, where V is the set of vertices, E is the set of edges, LV is the set of labels assigned to each vertex, and LE is the set of edge labels. Each label encodes a specific property or feature of an atom, ring, or bond. The set of features considered in this work is summarized in Table 1. It is evident that not every feature carries the same relevance; hence, we use weights to represent the relevance of each feature, which can be determined by an expert in the field. The weighting schemes used in this work are shown in Table 4. All features have been generated using RDKit 21.
2.0.1 Step 2: Solving the Maximum Co-k-Plex Problem
Given two molecular graphs, we are interested in finding the maximum co-k-plex of a third graph. This third graph is called a conflict graph and can be induced from the graphs being compared. Its construction is illustrated in Fig. 3.
Conflict Graph.
Formally, in a conflict graph, vertices represent possible mappings (or matchings), and edges represent conflicts between vertices. The notions of matches and conflicts are based on the criteria used to define similarity. Our algorithm allows the incorporation of different similarity criteria. For example, one criterion can be that two nodes match if they have the same atomic number; another criterion can be that two nodes match if they have the same pharmacophore feature, regardless of whether they are the same element. When comparing two molecular graphs to form a conflict graph, our algorithm compares only atoms to atoms and rings to rings. In addition, we allow both exact and partial matching. The term “partial match” refers to the matching of two nodes, where at least one of their labels matches exactly. An example of such matches is shown in Fig. 3.
In determining partial matches, we may wish to consider only a subset of matching features. To that end, we introduce the concept of criticality for each of the features described in Table 1. One is able to designate a feature as being “critical” or “non-critical” based on screening requirements. When comparing two nodes, if both nodes contain at least one matching feature marked as critical, they are included in the conflict graph. In the event that no feature is marked as critical, but both sets of features are an exact match, they will still be included in the conflict graph.
Maximum co-k-plex solution.
The objective of the maximum weighted co-k-plex problem is to identify the largest weighted set of vertices such that there are at most k−1 edges in the conflict graph. In this work, we set k=1, that is, the maximum weighted co-1-plex problem is equivalent to the maximum weighted independent set. The solution to the maximum weighted co-1-plex problem is a binary vector whose size is equivalent to the number of nodes in the conflict graph. The non-zeros in the solution vector identify a pair of vertices from the molecular graph. With this information, the solution can be mapped to the original molecular graphs to identify the common substructures.
The maximum weighted co-k-plex problem is formulated as a quadratic unconstrained binary optimization (QUBO) problem 13. Having the problem formulated as a QUBO problem allows us to use a quantum annealer; however, it is not our aim in this paper to study the annealer’s performance. Our focus is on investigating the performance of the GMS method in the context of ligand-based VS. In the “Results and Discussion” section we discuss in more detail the use of quantum annealers in ligand-based VS experiments. The solutions to the maximum weighted co-1-plex problems in this paper have been obtained using the classical heuristic parallel tempering Monte Carlo with isoenergetic cluster moves (PTICM) solver, also known as the “borealis” algorithm 22. The parameters used for the solver are given in Section S2 of the Supporting Information document.
2.0.2 Step 3: Similarity Measure
There exists a large number of similarity measures in the literature for graph-based molecular similarity methods. These measures are formulated in terms of subset relations between two graphs being compared and their common substructure. The similarity metric used in this work is based on the convex combination of two existing similarity measures—Bunke and Shearer, and asymmetric 23. In addition, we incorporate the information regarding the weights of each molecular feature. Given two molecules A and B, we denote the set of features of molecule A by FA and the set of features of molecule B by FB, and associate each feature f with a weight wf. After solving the maximum weighted co-1-plex problem, we recover the common features between molecules A and B and denote the set of common features by FAB. Hence, the similarity score is defined by
[TABLE]
with δ∈[0,1],
[TABLE]
3 Experimental Procedures
In order to perform an experimental validation of the GMS method for ligand-based VS, we have adopted standardized experimental procedures. This involves selecting an appropriate dataset and performance metric for a ligand-based VS approach. The experiments reported in this paper and the dataset used are based on the work of Jahn et al. 11, in which they tested the optimal assignment approach on molecular graphs as a ligand-based VS method. To evaluate the performance of VS using the GMS application, we report four metrics: the receiver operating characteristics (ROC) enrichment (ROCE), the area under the curve (AUC) of ROC, and their arithmetic weighted versions: awROCE and awAUC. In the sections that follow, we present further details regarding the dataset and performance evaluation used in our work.
3.1 Dataset
The chemical library used for screening is a modified version of the Directory of Useful Decoys (DUD): Release 2 24, 25, a standard dataset for benchmarking virtual screening, which contains a set of active structures for 40 target proteins. For each active compound, there are 36 inactive structures referred to as “decoys”. Decoys have similar physical properties but a dissimilar topology. We conduct experiments on the 13 target classes reported by Jahn et al. 11. The original DUD dataset is not suitable for ligand-based VS experiments, as it was designed for the evaluation of docking techniques. Good and Oprea 26 suggested and modified the active structures of the dataset using a lead-like filter and a clustering algorithm, aiming to reduce retrieval bias due to the presence of analogous structures.
The set of decoys was modified in a similar way by Jahn et al. 11 to reduce the introduction of artificial enrichment due to their having similar physical property values. As these actives and decoys contain only topological information, the researchers generated 3D atomic coordinates for each structure initially using CORINA 27, then optimized using MacroModel 9.6 28. We use the same modified set of actives and decoys in our experiments.
In order to screen for actives, ligand-based VS experiments require a query—a known biologically active structure that reacts with a target protein. Previous works utilized bound ligands of complexed crystal structure, extracted directly from the RCSB Protein Data Bank (PDB) 11, 29, 30. In situations where ligand-based VS is performed, the conformation adopted by a ligand upon binding to a receptor is often unknown; therefore, it seems advisable to test our ligand-based VS method using ligands with conformations generated by standard conformer models. In this work, we retrieve the 3D conformer of the query ligand from the PubChem website 31. The PubChem repository 32 generates a 3D conformer model that represents all possible biologically relevant conformations for a given molecule. Table 2 presents the 13 target classes; the number of actives, decoys, and clusters for each target class; the PDB code of the complexed crystal structure which contains the ligand; and the PubChem Compound ID number (CID) for the query. For all molecules in our dataset, we use RDKit to generate the molecular information. Note that for the target class FXa, there is a total of 64 actives, but RDKit is unable to generate one of the molecules.
3.2 Fingerprint Methods
To assess whether the GMS method has an advantage over traditional similarity methods in ligand-based VS experiments, we evaluate the performance of molecular fingerprint methods. Generally, there are two types of 2D fingerprints: dictionary-based and hash-based fingerprint methods 33. Dictionary-based methods have a fixed number of bits and each bit represents a certain type of feature of a substructure. Unlike dictionary-based fingerprints, hash-based fingerprints can be used to encode any new types of substructure features.
In this work, we select one of the most widely used fingerprints from each category. The dictionary-based fingerprint selected is the Molecular ACCess System (MACCS), which has 166 bits. The MACCS fingerprint is calculated using RDKit. The hash-based fingerprint is called a circular fingerprint, also known as an ECFP (extended-connectivity fingerprint). It uses a fixed number of iterations to generate identifiers for each atom based on its neighbours. At each iteration, only neighbours connected through a certain number of bonds are considered. After there are no further iterations, all atomic identifiers are collected, and all duplicated identifiers are removed. The resulting set of identifiers is then converted to a bit string with a hash function. In order to consider ECFP fingerprints, we use Morgan fingerprints generated by RDKit. Morgan fingerprints are built by applying the Morgan algorithm, whose default atom invariants use similar connectivity information as the one used for ECFP fingerprints. We also consider feature-based invariants information, similar to that used for functional-class fingerprints (FCFP). Landrum 34 presents details regarding the use of Morgan fingerprints as being equivalent to ECFPs/FCFPs.
3.3 Experimental Setup
As discussed in previous sections, our algorithm operates on a set of molecular features as defined in Table 1, where each feature is assigned a criticality and a weighting value. Details on the use of criticality and weighting values are given in the section “Similarity Methods”.
In order to gain a comprehensive understanding of the effect of the features considered in this work, we generate 12 criticality schemes (CS) and 10 weighting schemes (WS). A particular configuration of CS and WS is referred to as a similarity criteria setting. The detailed similarity criteria settings used in this work are shown in Section S1 of the Supporting Information document.
In both Tables 3 and 4, we present three selected CSs and WSs, respectively.
The CSs presented in Table 3 consider the features “Atomic number” and “Pharmacophore features” critical, that is, two nodes are considered a match if at least one of these features match. The rest of the features are considered non-critical. Three-dimensional coordinates are considered in CS7 and CS9, whereas CS3 does not consider this feature. The only difference between CS7 and CS9 is in how the feature that strictly compares atoms in the rings is set. CS7 imposes a critical constraint on ring comparison, whereas this condition is not imposed in CS9.
We consider a baseline weighting scheme WS (WSB) where each feature is weighted equally; hence, WSB acts a control. Also interested in studying the relevance of assigning a higher weight to rings than to individual atoms, we introduce WSB-5, where rings are given a weight of 5, which is proportional to the average number of atoms in a typical ring structure in our dataset. The rest of the WSs used in this work are selected based on previous knowledge; they have not been optimized.
For example, to set weighting scheme WS4, we consider that there are various types of non-covalent interactions between a ligand and receptor, like a salt bridge interaction, a hydrogen bond, an aromatic interaction, and a hydrophobic interaction. Generally, a salt bridge interaction has the highest interaction energy, followed by a hydrogen bond, an aromatic interaction, and a hydrophobic interaction. Therefore, we assign the highest weight to basic and acidic pharmacophore features; less weight to the hydrogen bond donors, acceptors, and aromatic centres; and the lowest weight to hydrophobic groups. In Table 4, we present these three selected WSs.
The following VS experimental procedure is performed for all similarity criteria settings: for each target class, a query is obtained from the PubChem database and compared against the respective set of active and decoy structures using the GMS method. For each comparison, we obtain the similarity score. The similarity scores between the query and the active and decoy molecules are then sorted to produce a ranking of molecules. Based on the ranking of molecules, we compute the VS performance measures awROCE, ROCE, awAUC, and AUC.
3.4 Performance Evaluation
Enrichment has traditionally been the standard measure used to characterize the performance of VS methods. It can be defined as “the ratio of the observed fraction of active compounds in the top few percent of a virtual screen to that expected by random selection” 35. The enrichment factor is, however, considered a poor performance measure because it depends on an extrinsic variable, that is, the ratio of active to decoy molecules. To address this issue, Jain and Nicholls 35 suggested two alternative measures. One is AUC, which is commonly used in other fields such as machine learning. Formally, AUC is defined as
[TABLE]
where N actives and N decoys is the number of actives and decoys in the dataset, respectively, and Ndecoys seeni is the number of decoy molecules that are ranked higher than the i-th active structure in the ranking list.
AUC values are a global measure that considers the whole dataset and, therefore, does not represent the concept of “early enrichment”. Hence, the other measure (ROCE) incorporates the concept. Specifically, it reports the ratio of true positive rates to false positive rates at different enrichment percentages. For example, “enrichment at 1%” is the ratio of actives observed with the highest 1% of known decoys (multiplied by 100). The ROC enrichment for a false positive rate (FPR) of x% is given by the expression
[TABLE]
where Nactives selectedx% and Ndecoys selectedx% is the number of actives (true positives) and decoys (false positives) retrieved in the range containing a false positive rate of x%, respectively. The enrichment percentages used in this work are 0.5%, 1.0%, 2.0%, and 5.0%.
Another important aspect to consider when evaluating the performance of a VS method is the retrieval of new scaffolds 11. Ligand-based methods could generate artificially higher enrichment results when the dataset contains analogous structures. To reduce this bias induced by structurally similar structures, Clark and Webster-Clark 36 proposed two weighted schemes for the standard ROC and AUC calculations: the harmonic weighted scheme and the arithmetic weighted scheme. In our work, we implement the arithmetic weighting scheme as done by Jahn et al. 11, denoted by awROCE. In order to use this metric, the actives in the dataset used for the VS experiment need to be clustered according to analogous structures. awRoce is given by the expression
[TABLE]
where Nc is the number of clusters, wi=Ni1 is the weight of each structure in the i-th cluster with Ni molecules, and Aix% is the total number of active structures retrieved from the i-th cluster at x% of the FPR.
Likewise, the arithmetic weighting scheme for AUC is denoted by awAUC and given by the equation
[TABLE]
where Ndecoys seenij is the number of decoys retrieved earlier than the i-th active molecule in cluster j, and Nj is the number of molecules in cluster j.
4 Results and Discussion
The performance of the GMS method was evaluated for each combination of 12 CSs and 10 WSs. The extensive set of results for the GMS method is presented in Sections S3 and S4 of the Supporting Information document.
In this section, we present the results for a selected set of similarity criteria settings of the GMS method that are representative of our analysis. Specifically, we include nine similarity criteria settings which consist of CS3, CS7, and CS9, each of them with the WSs WSB, WSB-5, and WS4.
In Table 5, we present the mean values for awAUC, awROCE, AUC, and ROCE over all 13 target classes shown in Table 2. The mean awROCE and mean ROCE are reported with decoy rates of 0.5%, 1.0%, 2.0%, and 5.0%. As mentioned earlier, awAUC and awROCE were introduced in order to reduce the inflated enrichments by considering structurally analogous molecules in the dataset. The similarity criterion CS9WSB-5 yielded the highest value for each of the metrics considered in this study. In particular, we observe that similarity criteria with WSB-5 generally result in higher scores than similarity criteria with the baseline WSB, suggesting it would be advisable to assign higher weights to rings.
4.1 Comparison against Other Methods
We also evaluated the performance of various fingerprint methods with the objective of comparing GMS against the most common method used for VS experiments. Specifically, we evaluated MACCS and Morgan fingerprints. In the case of Morgan fingerprints, we considered two numbers of bits, 1024 and 2048; four different radii, 1, 2, 3, and 4; and we kept the default option for the atom-based invariant feature. In summary, we evaluated 16 variations of Morgan fingerprints. The overall performance for fingerprint methods is presented in Table 11 in the Supporting Information document. Additional sets of results for each target class are detailed in Section S4.
Additionally, we considered the results of the optimal assignment methods presented by Jahn et al. 11. There are two main reasons for selecting these methods. First, optimal assignment methods act on molecular graphs; and second, the ligand-based VS results presented in this paper report awROCE as well as ROCE. The results for the optimal assignment methods were retrieved from their supplementary material, which is publicly available. We should note that their results were calculated with query molecules retrieved from the protein data bank and the authors corrected the bond lengths. The results using the optimal method could vary if the query were to have a different conformation. In Fig. 4, we show the mean awROCE and the mean ROCE at four percentage values for one variation for each of the GMS, Morgan fingerprint, and optimal assignment methods. The selected variation we report for each method is the one with highest mean value of awROCE at 0.5%. The Morgan fingerprint method, with 2048 bits, a radius of 4, and atom-based invariants, was selected from among the various fingerprint methods. Among the optimal assignment methods reported by Jahn et al. 11, a two-step hierarchical assignment approach (2SHA) was selected. From among the different criticality and weighting schemes of the GMS method, CS9WSB-5 was selected. Overall, the GMS method yielded higher enrichment values than the fingerprint and optimal assignment methods.
Table 6 shows the ROCE and awROCE at 0.5% enrichment for each of the 13 targets for the GMS and fingerprint methods.
Src—Tyrosine kinase.
In addition to the overall results already presented, we would like to highlight the results for the target kinase c-Src with the GMS method against fingerprint methods. In Table 7, we present the results for a GMS method and a Morgan fingerprint method. The GMS method yielded higher scores than the Morgan fingerprint method, in particular, awROCE at early enrichment at 0.5%. The awROCE metric was introduced to reduce the bias produced by comparing analogous structures. This high score of 80.22 indicates that the GMS method is able to distinguish the subtle conformational difference of the inhibitor for the c-Src kinase. Although a number of kinase inhibitors have been developed, designing a highly selective kinase inhibitor remains a challenge. We anticipate that our approach can help in the design of such selective inhibitors.
4.2 Impact of Three-Dimensional Coordinates as a Matching Feature
Previous studies comparing the performance of 2D versus 3D molecular similarity methods 19 have shown that 2D methods outperform 3D methods both in terms of achieving early enrichment and in computational time. In the GMS application, we have introduced a simple approach to using the 3D coordinates of atoms. From our VS experimentation, we observed that CSs with the 3D feature set to “ON” yield better results than CSs with the 2D feature set to “ON”. More specifically, the mean awROCE values at 0.5% ranged from 21.9 to 30.02 for CS1 to CS6, and from 31.11 to 43.39 for CS7 to CS12. The same was true for the mean ROCE values at 0.5%, whose values ranged from 30.37 to 37.66 for CS1 to CS6, and from 42.18 to 50.99 for CS7 to CS12.
In particular, let us consider the target class PDGFrb. One possible reason that could explain the poor performance of 3D methods on the PDGFrb dataset is the conformer optimization performed on the molecules in this dataset. In order to understand how different molecular conformations affect the performance of VS experiments, we generated various low-energy conformations of 10 active molecules out of 124, and selected the conformer that had the best overlay score. No modification was performed on the query molecule.
In Fig. 5(a), we show the overlap of the query PDGFrb and ZINC03832219 molecules with the conformation as retrieved from the DUD_LIB_VS_1.0 library. The similarity score obtained using the GMS method is 0.615, whereas for the Morgan fingerprint it is 0.729. Fig. 5(b) shows a conformation for the ZINC03832219 molecule with a higher degree of overlap: in this case, the similarity score using the GMS method increased to 0.887, whereas the Morgan fingerprint remained at 0.729.
Table 8 shows the VS performance of CS9WSB-5 for two PDGFrb datasets. One dataset is the complete original set with the molecules curated by Jahn et al. 11. The second set replaced 10 active molecules and the query molecule with modified conformations; we refer to this set as an “optimized dataset”. The optimized dataset significantly improved results for early enrichment; however, the overall scores for AUC and awAUC remained the same. In the future, it would be advisable to generate a certain number of conformations for all molecules, choosing the conformer that has the best similarity score.
The GMS method outperformed conventional techniques for most of the targets tested in this work by achieving higher early enrichment values, but a question remains regarding its computational cost. Finding the optimum solution for a GMS method comparison quickly grows impractical; however, heuristic methods can be used to find an optimal or near-optimal solution. A challenge that remains is to improve and scale these methods to more accurately and efficiently tackle molecular comparisons using large commercial databases. Quantum annealers have been argued to have the potential to take advantage of quantum phenomena such as quantum tunnelling 37 and entanglement 38 to solve optimization problems. The GMS algorithm has been formulated as a QUBO problem. Although this formulation is hardware agnostic, that is, classical optimization techniques can solve this problem, it also corresponds to the class of objective functions native to the quantum annealing hardware. This approach allows for implementation using the most effective solver (classical or quantum) available in a rapidly changing field.
Quantum annealers are a nascent, growing technology: every new generation shows a reduction in the amount of noise and an increase in the number of qubits (i.e., the basic unit of quantum information). The number of qubits is related to the number of variables an optimization problem can have. The most recent quantum annealer consists of 2048 qubits, and the number of variables (fully connected) that can be mapped to it is 66. The GMS problems generated using the DUD targets varied in size and density (a relation between the number of variables and their connections).
In Fig. 6, we illustrate the relation between the density and the number of variables for each problem in two cases. Fig. 6(a) and Fig. 6(b) correspond to the criticality schemes CS3 and CS9, respectively. Both schemes yielded problem sizes ranging up to 200 variables. The difference lies in their densities. Criticality scheme CS9 will generate denser problems, as it includes the feature of having 3D coordinates, whereas CS3 does not. From these figures, we can observe that it is not yet possible to solve all the generated problems using a quantum annealer. Therefore, a rigorous, exhaustive study of the quantum annealer’s performance is not yet feasible. The objective of this paper has been to evaluate the GMS formulation in terms of success metrics such as early enrichment. The study of the performance of this emergent and potentially groundbreaking technology is a subject for future work.
5 Conclusion
In this paper, we have presented a new methodology for determining molecular similarity, which enables the implementation of ligand-based virtual screening using a quantum computer. Our new graph-based molecular similarity (GMS) method solves the maximum weighted co-k-plex problem of an induced graph to find the maximum weighted common subgraph (MWCS). Solving the MWCS problem tends to be time consuming on classical computers, as it is, in general, NP-hard. Recent advances in quantum computing, and the availability of quantum annealing devices, offer alternatives for solving these classically hard problems 39. The GMS method has been formulated such that it can be implemented on quantum annealers. This formulation has enabled the future use of these devices when they have sufficiently improved (e.g., when their number of qubits and connectivity increase).
The advantage of using graphs is that they are able to encode any molecular information perceived as relevant. In our implementation of the method, we have incorporated both 2D and 3D descriptors into the molecular graphs. The highlight of the method is its flexibility, which allows the user to either include or exclude different molecular features according to their relevance to the problem. We tested this flexibility by implementing several combinations of descriptors. Various similarity criteria and combinations of molecular features used in the GMS method were evaluated on 13 datasets from the DUD_LIB_VS_1.0 library. We identified a particularly successful configuration of features that uses an equal weighting of 1 for all features except the rings; the rings were assigned a higher weight of 5. Our results demonstrate that the GMS method, where similarity criteria include three-dimensional coordinates, has, in general, a higher early enrichment value than the similarity criteria that do not include this feature. We also found the results to be sensitive to the 3D conformation of the molecules used. Further research is needed to find the optimal conformation of query molecules that do not have crystal structures in complex with target proteins in an unbiased manner. Overall, we conclude that the GMS method outperforms conventional fingerprint and optimal assignment methods for most of the 13 DUD targets.
{acknowledgement}
The authors thank Arman Zaribafiyan and Dominic Marchand for their early contributions to the project; Takeshi Yamazaki for valuable discussions, collaboration on molecular figures, and suggestions for revisions to the manuscript; Rudi Plesch, Nick Condé, and SungYeon Kim for their contribution to the development of the GMS application; Jayne Drew for digitizing and editing the GMS method figures; Shawn Wowk for his support; Kausar N. Samli for his comments on an early version of the paper; and Marko Bucyk for his helpful comments and for editing the manuscript. Additionally, the authors thank Giuliano Berellini for insightful discussions at the beginning of the project; Jeff Elton for his expertise and guidance; Simon Sekhon, Brandon Hedrick, and Spencer Herath for their work on the proof-of-concept molecular comparison application; and Charles Rozea, Teresa Tung, and Daniel Garrison for their support.
6 Supporting Information
6.1 S1 Criticality and Weighting Schemes
For each feature presented in Table 1, we have set its value to critical (C), non-critical (NC), or off (OFF). Bond order, formal charge, and degree have been set to NC and pharmacophore features have been set to C. For the remaining set of features, we have set various combinations of values as presented in Table 9. Each combination is called a criticality scheme (CS). In total, we have generated 12 CSs.
In addition to identifying the features as C or NC, we have assigned them a weighting value to reflect its relevance in the virtual screening (VS) experiments. Each combination of weighting values is called a weighting scheme (WS). We set one WS as baseline WSB, with every feature equally weighted to act as a control. In total we have 10 WSs, shown in Table 10.
6.2 S2 Solver Parameters
The algorithm used to solve the maximum weighted co-1-plex problem is the parallel tempering Monte Carlo with isoenergetic cluster moves (PTICM) heuristic solver. The performance of replica-exchange algorithms, such as the PTICM algorithm, depends on the parameters used, especially the temperature schedule selected. In the main paper, the temperature schedule has been based on the geometric schedule for each replica. The low and high temperature values are determined based on the values of the coefficients of each quadratic unconstrained binary optimization problem instance, and the number of replicas has been chosen to be two.
6.3 S3 Overall VS Performance for the GMS and Fingerprint Methods
6.4 S4 VS Performance Across 13 targets in the DUD_LIB_VS_1.0 Library
6.4.1 ACE
6.4.2 AChE
6.4.3 CDK2
6.4.4 COX-2
6.4.5 EGFr
6.4.6 FXa
6.4.7 HIVRT
6.4.8 InhA
6.4.9 P38 MAP
6.4.10 PDE5
6.4.11 PDGFrb
6.4.12 Src
6.4.13 VEGFr-2
Figure 1
Figure 2
Figure 1
Figure 2
Figure 3
Figure 6
Figure 7
Figure 8
Figure 9
Figure 10
Figure 11
Figure 12
Figure 13
Figure 14
Figure 15
Figure 16
Figure 17
Figure 18
Figure 19
Figure 20
Figure 21
Figure 22
Figure 23
Figure 24
Figure 25
Figure 26
Figure 27
Figure 28
Figure 29
Figure 30
Figure 31
Figure 32
Figure 33
Figure 34
Figure 35
Figure 36
Figure 37