In this section, we first conduct the experiments on a benchmark to demonstrate ParetoDrug’s remarkable ability to generate molecules with multiple desired properties including the binding affinity and drug-like properties when compared with various baselines. Meanwhile, we also give the statistical analysis of the generated molecules of ParetoDrug. Then we use ParetoDrug to perform the case studies for the multi-objective target-aware drug discovery task, multi-target drug discovery task, and multi-target multi-objective drug discovery task respectively. In these case studies, ParetoDrug is able to generate the Pareto Dominate molecules over the known drug ligands in terms of the specified molecule property objectives, which exhibits the promising molecule discovery potential of ParetoDrug.Benchmark experimentsIn the benchmark experiments, we follow the settings as Qian et al.15 where there are 100 protein targets sampled from the public database of protein-ligand pairs BindingDB31 as the test set. For each test protein target, we generate 10 candidate molecules for evaluation. All 1000 candidate molecules are evaluated by a set of molecule property metrics, and the scores are averaged for an overall comparison. Please refer to Supplementary Information A and B for a detailed experimental and hyperparameter setup. We use several important metrics to evaluate the generated molecules, including docking score, uniqueness, LogP, QED, SA score, and NP-likeness described as follows.
Docking score. Binding energy is regarded as a general indicator to describe the binding affinity between molecule ligands and target proteins. Specifically, we utilize a free and widely used tool called smina32 to compute the binding affinity. We use the negative value of the output by smina as the docking score. The higher the docking score is, the better the molecule is docked into the target protein.
Uniqueness. Drug design models should be able to generate different molecules conditioning on different target proteins. The higher the uniqueness value is, the more sensitive the model is to the specified target protein. This metric is computed as follows:$${{{{\rm{Uniqueness}}}}}( \% )=\frac{\#({{{{\rm{Set}}}}}({\cup }_{{S}_{{{{{\rm{p}}}}}}\in {{\mathbb{S}}}_{{{{{\rm{p}}}}}}}{{{{\rm{Set}}}}}({M}_{{s}_{{{{{\rm{p}}}}}}})))}{\#({\cup }_{{S}_{{{{{\rm{p}}}}}}\in {{\mathbb{S}}}_{{{{{\rm{p}}}}}}}\,{{{{\rm{Set}}}}}({M}_{{S}_{{{{{\rm{p}}}}}}}))}\times 100 \% ,$$
(1)
where \({{\mathbb{S}}}_{{{{{\rm{p}}}}}}\) indicates the set of test proteins, \({M}_{{S}_{{{{{\rm{p}}}}}}}\) denotes the collection of generated molecules for the target protein \({S}_{{{{{\rm{p}}}}}}\in {{\mathbb{S}}}_{{{{{\rm{p}}}}}}\), # counts the number of molecules, and Set is an operator to remove the repeated molecules in the given set.
LogP. A large LogP value indicates the substance is lipophilic, while a small LogP value means it is easy to dissolve in water. According to Ghose filter33, the LogP value of a druggable molecule should range from −0.4 to +5.6.
QED. This score measures the drug-likeness and ranges from 0 to 1. A higher QED score indicates that a molecule is more likely to be a potential drug-like compound, with the desired molecular properties such as hydrogen bond acceptor, hydrogen bond donor, and polar molecular surface area34.
SA score. The synthetic accessibility (SA) score indicates how difficult one molecule is to synthesize, which is calculated based on a combination of fragment contributions and a complexity penalty35. The range of the estimated SA metric is from 1 (easy to make) to 10 (very difficult to make).
NP-likeness. Natural products play an important role in the history of drug discovery. Many drugs are natural products and their derivatives. The higher the score is, the more likely the molecule is to be a natural product. The calculated NP-likeness is typically in the range from -5 to 536.
The reported results of “Known ligands”, SBMolGen, LiGANN, SBDD-3D, and BeamLmser are from AlphaDrug15. The “Known ligands” indicates the original molecules binding to protein targets in the database. The results of LiGANN10 were collected on the web-based application provided in the original paper. SBMolGen37 is developed from ChemTS38 for target-specific molecular generation. The results of SBDD-3D18 were based on the released codes and trained model published by the authors. BeamLmser applies the beam search on the pretrained Lmser Transformer15. The beam size of BeamLmser is set at 10 to collect 10 molecules for each test protein target. Besides the above representative baselines, we also test three recent advanced methods. The first is Pocket2Mol20, which uses the equivariant generative network and autoregressive sampling scheme to generate three-dimensional molecules. For Pocket2Mol, we utilize the official codes and trained model for sampling molecules. The second is TargetDiff12, which develops a three-dimensional equivariant diffusion model to sample molecules. For TargetDiff, we also use the officially released trained model and codes for sampling. We keep the sample numbers of Pocket2Mol and TargetDiff at 100 for each test protein, which is the default configuration to ensure the quality of generated molecules. To make a fair comparison with other methods, for each test protein target, we randomly select 10 molecules from the generated 100 molecules of Pocket2Mol and TargetDiff for the evaluation. The third is CProMG28, which proposes a multi-constraint autoregressive model to generate small molecules with controllable properties. We use the official codes and default configurations of CProMG to generate 10 molecules for each test protein with the pretrained CProMG-VQSLT model, which is trained to control multiple property metrics including the docking score, LogP, QED, and SA score that are evaluated here.Besides the above basic generative models, there also emerges another kind of fundamental approach that integrates the powerful MCTS-based searching technique to better control the molecule generation procedure of the pretrained autoregressive generative models with the simulation feedback, and AlphaDrug and the proposed ParetoDrug fall into this kind. For AlphaDrug15 which utilizes MCTS with the pretrained Lmser Transformer model to generate molecules based on given protein targets, we run the official codes and set iteration times (IT) at 150 when selecting the next atom symbol in MCTS. For ParetoDrug which conducts Pareto MCTS with the same pretrained Lmser Transformer model, we also set IT at 150 and let it optimize all objectives (docking score, LogP, QED, SA score, and NP-likeness) synchronously except the unoptimizable Uniqueness, which is a statistic metric for all generated molecules. In addition, we set the metric value of LogP as 1 if the molecule’s LogP value is in the range of [ − 0.4, 5.6], and 0 otherwise. After each Pareto MCTS, ParetoDrug obtains a global pool of Pareto optimal molecules. We choose the molecule with the largest reward vector summation value from the pool, which means this molecule has top rankings in each property metric. When testing, we collect 10 generated molecules for each test protein target from AlphaDrug and ParetoDrug.Additionally, we compare a multi-objective drug discovery algorithm REINVENT 439 while its generation model is not conditioned on the protein information. It uses a reinforcement learning algorithm to generate optimized molecules compliant with a user-defined property profile defined as a multi-component score. We let REINVENT 4 optimize the docking score, LogP, QED, SA, and NP while setting their weights in the property profile all at 0.2. For each test protein target, we collected 10 molecules with the highest multi-component scores during the training process of REINVENT 4.The results are shown in Table 1 and the direction of the arrow in the table means a better property score. The 95% confidence intervals for property scores of RL/MCTS are included. As we see, in terms of the docking score, ParetoDrug demonstrates superiority over all baselines except AlphaDrug. However, AlphaDrug is a single-objective target-aware drug discovery method that only optimizes the binding affinity. As AlphaDrug and ParetoDrug have the same iteration budgets (IT=150) for each atom symbol in sequence but ParetoDrug needs to optimize multiple objectives including the binding affinity, it is expected that ParetoDrug has a lower docking score than AlphaDrug. Meanwhile, although the docking score of ParetoDrug decreases slightly, other metrics including QED, SA score, and NP-likeness are improved significantly compared with AlphaDrug. Notably, QED changes from 0.4 to 0.6 (50% improvement) while NP-likeness changes from -0.9 to -0.4 (55.6% improvement). For the special LogP metric, although the average LogP value of AlphaDrug falls into the druggable molecule range, only 52.7% generated molecules of AlphaDrug satisfy the LogP range constraint if tested individually. On the contrary, 96.5% (83.1% improvement over AlphaDrug) generated molecules of ParetoDrug satisfy the LogP range constraint. These impressive results demonstrate that ParetoDrug is able to address the multi-objective target-aware drug discovery task by discovering novel compounds that possess multiple satisfactory properties including the binding affinity. On the other hand, we observe that the pretrained autoregressive Lmser Transformer with beam search (BeamLmser) cannot generate molecules with higher docking scores than the most recent TargetDiff. But with MCTS replacing beam search, AlphaDrug greatly boosts Lmser Transformer’s performance to find molecules with stronger binding affinity than BeamLmser even with the same docking time budgets15. Furthermore, ParetoDrug proposes the multi-objective Pareto MCTS to replace the MCTS used in AlphaDrug. With the same iteration times, ParetoDrug significantly improves multiple molecule properties compared with AlphaDrug while maintaining the docking score at the same level. Additionally, when compared with the multi-constraint conditional generation method CProMG, ParetoDrug has advantages in docking score, Uniqueness, SA score, and NP-likeness. In addition, the Uniqueness of CProMG is only 26.9% as it generates the same molecules for different protein targets, which is undesirable in de novo target-aware drug discovery tasks. Lastly, for REINVENT 4 which does not belong to the kind of target-aware drug discovery methods, we could see although it achieves superior performance in some metrics such as QED and NP, its docking score is much lower than ParetoDrug as it does not encode the protein-ligand prior to its generation model. This also indicates the importance of incorporating the protein target information into the molecule generation process as in the generative target-aware drug discovery methods.Table 1 Average metric scores of generated molecules of each method (n = 1000 molecules) on the sampled 100 test proteinsNext, we conduct the statistical analysis with kernel density estimate40, which is analogous to a histogram but endowed with benefits such as smoothness and continuity. The property distributions of molecules generated by TargetDiff, AlphaDrug, and ParetoDrug are shown in Fig. 1. For TargetDiff, here we use 10 molecules with the highest docking scores among the generated 100 molecules for each test protein to make an aligned comparison. We can see that although the docking score distributions of the three methods are similar while AlphaDrug is slightly better, other property distributions present differently. For LogP, ParetoDrug satisfies the range constraint of [ − 0.4, 5.6] while TargetDiff tends to generate more molecules with LogP values below the lower bound and AlphaDrug tends to generate more molecules with LogP values above the upper bound. Besides, ParetoDrug is able to generate more molecules with high QED values than the other two methods especially when QED is larger than 0.8 that molecules are very likely to be potential initiators of a drug candidate. Meanwhile, TargetDiff’s molecules are with significantly higher SA values than ParetoDrug, which means that TargetDiff’s molecules are much harder to synthesize. These statistical findings demonstrate that ParetoDrug has better molecule distributions than AlphaDrug and TargetDiff when taking multiple properties into account. More comparisons of computational efficacy, score distributions of a specific target, and the diversity of generated molecules between ParetoDrug and other methods could be referred to Supplementary Information C (and Supplementary Table 1), D (and Supplementary Fig. 1), and E (and Supplementary Table 2).Fig. 1: The molecular property distributions of generated molecules (n = 1000 molecules) by ParetoDrug, AlphaDrug, and TargetDiff respectively.A The docking score (kcal ⋅ mol−1) distributions of each method. B The LogP value distributions of each method. C The QED value distributions of each method. D The SA score distributions of each method. E The NP-likeness distributions of each method.Case studies for multi-objective target-aware drug discoveryHere we use two case studies of disease protein targets to show the molecule discovery ability of ParetoDrug for the multi-objective target-aware drug discovery tasks. The molecule objectives optimized by ParetoDrug are the docking score, LogP, QED, SA score, and NP-likeness. Additionally, the binding affinity is further validated by MM-GBSA41,42, which is a more accurate metric than docking scores but computationally expensive. For the analysis of the protein-ligand interactions, we use PLIP43 and detailed can be referred to Supplementary Information F.Case 1: targeting FXRNon-alcoholic fatty liver disease (NAFLD) is defined as the excessive and abnormal intracellular accumulation of lipids in the liver, primarily in the form of triglycerides44,45. Currently, NAFLD has been the most common cause of chronic liver disease, especially in Western countries, and the estimated prevalence of NAFLD is approximately 30% in the general population46,47. One of the best-known drugs for NAFLD is Tropifexor, which acts as an agonist of the farnesoid X receptor (FXR). The structural basis of Tropifexor as a potent and selective agonist of FXR is shown in Fig. 2A (PDB ID: 7D42)48. In this case study, we use ParetoDrug to discover potential drug molecules with desired computational properties for FXR. Using ParetoDrug, we collect 10 molecules and find four Pareto Dominate molecules compared with Tropifexor. The chemical structures of Tropifexor and the discovered ligands by ParetoDrug for FXR are shown in Fig. 2B. Table 2 shows the property metrics of different ligands. Our ParetoDrug model discovers multiple ligands that outperform Tropifexor on all the optimized properties. Especially, the SA scores of the new ligands are much lower than Tropifexor, which means that they are easier to synthesize. We also run AlphaDrug and TargetDiff to collect molecules for FXR, however, no Pareto Dominate molecule over Tropifexor is found for the two methods. For example, the best molecule from AlphaDrug (with the most number of better properties than Tropifexor) has the Docking score at 12.9, LogP at 4.1, QED at 0.3, SA at 2.6, and NP at -1.59. Compared with it, Compound 4 generated by ParetoDrug has 3 better properties (Docking score, QED, and NP) and 1 worse property (SA). Meanwhile, the best molecule from TargetDiff (with the most number of better properties than Tropifexor) has the Docking score at 12.8, LogP at 4.9, QED at 0.48, SA at 7.7, and NP at 1.01. Compared with it, Compound 4 generated by ParetoDrug has 3 better properties (Docking score, QED, and SA) and 1 worse property (NP). As shown in Fig. 2C, the docked poses and interactions of these four discovered compounds are quite different compared with Tropifexor. More specifically, one hydrogen bond forms between Tropifexor and the amino-acid residue MET265. At the same time, Compounds 1, 3, and 4 with new scaffolds form new hydrogen bonds with other residues (Compound 1 with THR288 and TYR369, Compound 3 with HIS294, and Compound 4 with THR288) while no hydrogen bond forms between Compound 2 and FXR.Fig. 2: Static structural analysis of ligands binding to FXR (PDB ID: 7D42).Tropifexor is an agonist of FXR and Compounds 1–4 are found by ParetoDrug. Hydrogen bonds are displayed in yellow dashed lines and π-π interactions are in red. A Solvent-accessible surfaces of the binding pocket of FXR for Tropifexor and Compounds 1 to 4. B Chemical structures of Tropifexor and Compounds 1–4. C The binding poses of Tropifexor and Compounds 1–4.Table 2 Metrics of generated molecules for protein target FXRBesides the docking score, MM-GBSA rescoring based on molecular dynamics simulations is used to further computationally validate the discovered compounds41. MM-GBSA uses molecular mechanics with generalized Born surface area to determine highly potential inhibitors for targets. The detailed settings of MM-GBSA are provided in Supplementary Information G. In this case, we use MM-GBSA to further validate the generated Pareto optimal molecules with promising docking scores. As shown in Table 2, the MM-GBSA scores indicate that the generated molecules have the same level of binding free energies as the known drug Tropifexor. Surprisingly, although two hydrogen bonds formed between Compound 1 and FXR, the MM-GBSA scores show no significant difference. Possibly because hydrogen bonding interaction is ignored in the MM-GBSA calculation.Case 2: targeting PI3K-γ
Follicular lymphoma (FL) is a systemic neoplasm of the lymphoid tissue displaying germinal center B-cell differentiation, which belongs to a cancer that involves certain types of white blood cells known as lymphocytes. FL represents 5% of all hematological neoplasms and about 20-25% of all new non-Hodgkin lymphoma diagnoses in Western countries49. One of the best-known drugs for FL is Copanlisib, which has been shown to affect the survival and spread of cancerous B-cells. The structural basis of the PI3K-γ related to FL in complex with Copanlisib is shown in Fig. 3A (PDB ID: 5G2N)50. Here we use ParetoDrug to discover potential drug molecules with desired computational properties for PI3K-γ. We collect 10 molecules and find one Pareto Dominate molecule (Compound 5) compared with Copanlisib (Fig. 3B). As shown in Fig. 3C, Compound 5 found by ParetoDrug has three hydrogen bonds with surrounding residues (VAL882, ASP836, and LYS833) and two π-π stackings with surrounding residues (TRP812 and TYR867). Notably, the hydrogen bonds to VAL882 and LYS833 as well as π-π stacking to TYR867 also appear in Copanlisib’s docking interactions. Meanwhile, Table 3 shows the computational metric values of Copanlisib and Compound 5. Compound 5 is better than Copanlisib in terms of the optimized molecule metric objectives. However, the MM-GBSA scores indicate that the binding strength of Compound 5 decreases compared with Copanlisib (-46.48 kcal ⋅ mol−1 vs. -55.51 kcal ⋅ mol−1). The possible reason is that hydrogen bonding and π-π interactions are not considered in the energy terms of MM-GBSA. We also run AlphaDrug and TargetDiff to collect molecules for PI3K-γ and no Pareto Dominate molecule over Copanlisib is found by the two methods. For example, the best molecule from AlphaDrug (with the most number of better properties than Copanlisib) has the Docking score at 12.4, LogP at 4.35, QED at 0.41, SA at 5.0, and NP at 1.29. Meanwhile, the best molecule from TargetDiff (with the most number of better properties than Copanlisib) has the Docking score at 11.6, LogP at 3.2, QED at 0.56, SA at 3.9, and NP at 0.18. Although with some good properties, the two top molecules from AlphaDrug and TargetDiff cannot dominate the drug Copanlisib with worse SA.Fig. 3: Static structural analysis of ligands binding to PI3K-γ (PDB ID: 5G2N).Copanlisib is a drug binding to PI3K-γ and Compound 5 is found by ParetoDrug. Hydrogen bonds are displayed in yellow dashed lines and π-π interactions are in red. A Solvent-accessible surfaces of the binding pocket of PI3K-γ for Copanlisib and Compound 5. B Chemical structures of Copanlisib and Compound 5. C The binding poses of Copanlisib and Compound 5.Table 3 Metrics of the generated molecule for protein target PI3K-γWhile the in silico computational metrics of molecules discovered by ParetoDrug show promise in comparison to existing drugs, it is crucial to acknowledge that these molecules are still far from being drugs. Drug discovery is an extremely complicated process, and the current metrics for molecules cannot perfectly reflect the physicochemical properties required for a compound to be a drug. Nevertheless, we clearly see ParetoDrug’s promising potential in addressing multi-objective target-aware drug discovery tasks.Case study for multi-target drug discoveryMulti-target drug discovery can be considered a special case of multi-objective drug discovery where each protein target is going to be regarded as an objective to optimize. Until now, the study of multi-target target-aware drug discovery remains underexplored as it is challenging to consider the information of multiple protein targets at the same time to derive one ligand that could bind to all these given targets. Meanwhile, previous generative target-aware drug discovery works mainly focus on the single-target situation as there lack data sets to train the multi-target conditioned generative models. In this case study, we use ParetoDrug to perform a multi-target target-aware drug discovery task to design dual-functional inhibitors for both the HIV protease (HIV-PR) and HIV reverse transcriptase (HIV-RT). ParetoDrug is slightly modified to be compatible with this kind of task, and details are given in the Method section.The crystal structures of HIV-PR and HIV-RT used here are 3A2O and 4G1Q51. Both structures are complexes with potent inhibitors solved at high resolution. We compare with LigBuilder V351, the first de novo multi-target drug design program, and the variants of Pocket2Mol and TargetDiff extended by combining with screening. There are three different strategies in LigBuilder V3, including multi-target de novo design, multi-target growing, and multi-target linking. The best molecules for each strategy from the original paper are reported here. For the variants of Pocket2Mol and TargetDiff, we use each method to generate 100 molecules for each target. Then we use smina to screen the generated 200 molecules of each method to find the best molecule that has the best docking scores for both targets. We call the variants as Pocket2Mol-screen and TargetDiff-screen, and report the best molecules of each variant.Results of both the docking scores and MM-GBSA scores for each method’s best molecule are shown in Table 4. For ParetoDrug, we additionally report another two top molecules (Compound 7 and Compound 8) which are almost the same good as Compound 6 in terms of the docking scores. As shown in Table 4, the promising docking scores and MM-GBSA scores of Compounds 6–8 demonstrate that they are potential strong dual inhibitors to the given two protein targets in this task. Additionally, Fig. 4 shows Compounds 6–8 and their docking poses and interactions with HIV-PR (PDB ID: 3A2O) and HIV-RT (PDB ID: 4G1Q). Interestingly, the similar structures of Compounds 6 to 8 in Fig. 4B indicate that ParetoDrug found a chemical subspace of strong inhibitors for both the HIV-PR and HIV-RT targets.Table 4 Docking and MM-GBSA scores of the generated molecules by baselines and ParetoDrug for protein targets HIV-PR and HIV-RTFig. 4: Static structural analysis of ligands binding to both the HIV-PR (PDB ID: 3A2O) and HIV-RT (PDB ID: 4G1Q).Compounds 6–8 are found by ParetoDrug. Hydrogen bonds are displayed in yellow dashed lines. A Solvent-accessible surfaces of the binding pockets of HIV-PR (left) and HIV-RT (right) for Compounds 6–8. B Chemical structures of Compounds 6–8. C The binding poses of Compounds 6 to 8 with HIV-PR. D The binding poses of Compounds 6–8 with HIV-RT.Case study for multi-target multi-objective drug discoveryWe have shown ParetoDrug’s promising ability for both the multi-objective target-aware drug discovery and multi-target drug discovery tasks separately. Naturally, a more attractive and challenging task is the multi-target multi-objective drug discovery task where the generated molecules need to bind to a given set of protein targets while manifesting other desired computational molecule properties. To the best of our knowledge, there are no published works yet to specifically address this kind of task. Here we use Lapatinib as a case study to evaluate ParetoDrug’s ability for the multi-target multi-objective drug discovery task. Lapatinib is a dual tyrosine kinase inhibitor that interrupts the EGFR pathway and inhibits HER4/ErbB4 Kinase for the treatment of breast cancer52,53. 1XKK and 3BBT are respectively the PDB IDs of crystal structures for Lapatinib binding to EGFR54 and the HER4/ErbB4 kinase55. In this task, we configure ParetoDrug to bind to both protein targets while optimizing LogP, QED, SA, and NP metrics synchronously with IT at 150. We compare the generated molecules in the global Pareto pool with the known dual-inhibitor drug Lapatinib. Impressively, plenty of Pareto Dominate molecules over Lapatinib are discovered by ParetoDrug, as shown in Fig. 5. Furthermore, the QED values of Compound 12 and Compound 14 are greater than 0.8, which indicates the two molecules are potential initiators for a drug.Fig. 5: Property metric values of Lapatinib and its Pareto Dominate molecules (Compounds 9 to 15) found by ParetoDrug.EGFR means the ligand’s docking score to EGFR (PDB ID: 1XKK) while HER4/ErbB4 means the ligand’s docking score to HER4/ErbB4 kinase (PDB ID: 3BBT), and the unit for the two docking scores is kcal ⋅ mol−1.