From Static to Dynamic Structures: Improving Binding Affinity Prediction with Graph-Based Deep Learning

TL;DR

Dynaformer, a graph‐based deep learning model, is developed to predict the binding affinities of protein‐ligand binding affinities by learning the geometric characteristics of the protein‐ligand interactions from the MD trajectories, and demonstrated that the approach offer a promising avenue for accelerating the early drug discovery process.

Abstract

Accurate prediction of protein-ligand binding affinities is an essential challenge in structure-based drug design. Despite recent advances in data-driven methods for affinity prediction, their accuracy is still limited, partially because they only take advantage of static crystal structures while the actual binding affinities are generally determined by the thermodynamic ensembles between proteins and ligands. One effective way to approximate such a thermodynamic ensemble is to use molecular dynamics (MD) simulation. Here, an MD dataset containing 3,218 different protein-ligand complexes is curated, and Dynaformer, a graph-based deep learning model is further developed to predict the binding affinities by learning the geometric characteristics of the protein-ligand interactions from the MD trajectories. In silico experiments demonstrated that the model exhibits state-of-the-art scoring and ranking power on the CASF-2016 benchmark dataset, outperforming the methods hitherto reported. Moreover, in a virtual screening on heat shock protein 90 (HSP90) using Dynaformer, 20 candidates are identified and their binding affinities are further experimentally validated. Dynaformer displayed promising results in virtual drug screening, revealing 12 hit compounds (two are in the submicromolar range), including several novel scaffolds. Overall, these results demonstrated that the approach offer a promising avenue for accelerating the early drug discovery process.

Figures (5)

• Figure 1: An illustrative overview of the thermodynamic ensemble and the MD trajectory dataset.A) The conformations in a thermodynamic ensemble generally follow a specific distribution. In principle, the binding affinity is determined by the properties of the thermodynamic ensemble rather than a single conformation. Based on this fact, we built a dataset containing various thermodynamic ensembles sampled by MD simulations. More details can be found in the main text. B) Relationship between the mean ligand RMSD and the binding affinity in the MD trajectory dataset. The mean ligand RMSD is calculated as the average distance deviation of the ligand from its original position in the crystal structure during the simulation. The heatmap and the mean RMSD range in the inset indicate that the ligand dynamics are related to binding affinities, i.e., the higher the mean ligand RMSD, the lower the binding affinity. C) Representative examples of ligand RMSD and their binding affinities, including unstable (PDB ID: 4y3j), intermediate (PDB ID: 3udh), and stable (PDB ID: 2yge) MD trajectories. In the unstable case, the ligand RMSD increased sharply at the 60th snapshot, indicating that the ligand left the binding site. The ligand with intermediate stability remained relatively flexible but stayed near the binding site. In the stable case, the ligand stuck tightly to the binding pocket, yielding low ligand RMSDs.
• Figure 2: Illustration of our virtual screening pipeline with Dynaformer.A) Hit discovery through our deep learning model Dynaformer based on an MD trajectory dataset. First, crystal structures of protein-ligand complexes are collected to form an MD trajectory dataset. Then, the simulated trajectories are utilized for training the binding affinity prediction model. Next, the trained model is evaluated against a benchmark dataset. The model with the best performance is then used as a scoring function for the virtual screen of the compound candidates. Finally, the wet-lab experiment validates the hit molecules for further hit-to-lead optimization. B) An example snapshot of protein-ligand from the MD trajectory dataset and its graph representation. In order to create such a graph representation to feed into the model, atoms from both the ligand and the protein are chosen based on their proximity to the ligand within a specified distance cutoff. The nodes and the edges represent the atoms and their covalent or non-covalent interactions, respectively. C) The architecture of Dynaformer. The node features (i.e., atomic features) are first fed into a multi-head attention module. The structural encodings, including distance, angle, and edge features, are encoded and added as the attention bias. After the final layer, the feature representations of the graph are globally pooled and fused with the pre-calculated fingerprints, which are knowledge-based features capturing structural and chemical properties. The fused representations are then used for the final prediction of the binding affinity. Abbreviations in the figure: Linear, linear layer; Feed-forward, feed-forward neural network; MatMul, matrix multiplication.
• Figure 3: Performance evaluation of binding affinity prediction.A) Scoring power evaluation on the 285 CASF-2016 complexes. The scoring power illustrates the inter-target performance of how well the predictions correlate with the experimentally measured binding affinities. Pearson coefficient (Pearson $\bm{r}$), root mean squared error (RMSE), and standard deviation (SD) are used as evaluation metrics. B) Ranking power evaluation on the 57 classes of target proteins from CASF-2016, in which each target bound to five different ligands. The goal of evaluating the ranking power is to assess the intra-target performances of scoring functions in terms of how accurately they predict the relative order of molecules. Evaluation metrics are measured in terms of the average values of Spearman's coefficient (Spearman $\rho$), Kendall's coefficient (Kendall $\tau$), and predictive index (PI) on each target. C) A detailed scatter plot of the Dynaformer predicted $\hbox{p}K_i$ and experimentally measured $\hbox{p}K_i$ values of the 285 protein-ligand complexes from CASF-2016, demonstrating the effectiveness of Dynaformer in predicting binding affinity with high correlation and low prediction bias. D) Ablation studies of Dynaformer. Here, -F indicates a modified version of Dynaformer without any pre-calculated fingerprints; -3D stands for a modified version of Dynaformer without the structural encoding module; -M stands for a modified version of Dynaformer that is trained without MD trajectories, i.e., only crystal structures are used for training; The results showed consistent improvements in both scoring and ranking power by incorporating MD data.
• Figure 4: Illustrative examples of Dynaformer trained with MD trajectories improves binding affinity prediction.A) An example (PDB ID:2v7a) of improved affinity prediction through incorporating additional entropic information. The head of the ligand is tightly bound in the pocket, resulting in a significant favorable enthalpy change. The random movement of the solvent-exposed ligand tail contributes to a favorable entropy change during ligand binding. B) An example (PDB ID: 3udh) of improved affinity prediction through better protein-ligand interaction modeling. With a rigid structure and no rotatable bonds, the binding affinity relies primarily on the enthalpy changes from binding site interactions. Despite a few favorable interactions, subpockets around the ligand remain unoccupied. Therefore, the ligand moves away from its initial position during the simulation, leading to low binding affinity. The ability of Dynaformer to model these interactions results in better prediction. C) An example of ranking ligands with subtle interaction differences. The three ligands (PDB IDs: 2qbr, 2qbq and 2qbp) have the same scaffold (the shaded area) but different tail structures. Dynaformer better addresses activity cliffs by differentiating between ligands sharing the same scaffold but with varying binding affinities. D) Conformational stability of the three ligands shown in C. The 2qbp ligand exhibited high stability during the simulation, while the 2qbq ligand partially left the original binding site after 1 ns and the 2qbr ligand left the binding site after 9 ns. The varying binding affinities of these ligands can be attributed to the hydrophobic interactions with the second subpocket on the right side shown in C. E) The interaction fraction analysis of the three ligands shown in C. The percentages of simulation snapshots with different protein-ligand interactions are shown. Higher binding affinity corresponded to more stable interactions (e.g., $\pi-\pi$ stacking, hydrogen bond and hydrophilic interaction) observed during the simulation. With such interaction features and better modeling capability, Dynaformer can provide more accurate ligand rankings. ECIF and SIGN are selected as baselines for comparison in A-D.
• Figure 5: Virtual screening against the target HSP90 and the corresponding experimental validation.A) Structures, experimentally measured $K_i$ values, and the SPR sensorgrams of the top 12 molecules with measurable binding affinities from the virtual screening. B) Detailed $K_i$ values of the selected 12 molecules and their ranking statistics by Dynaformer and Autodock Vina. To the best of our knowledge, compounds 1, 3, 4, 6, and 12 show promise as novel candidates for hit-to-lead optimization. C) The correlations of the predictions by Dynaformer and Autodock Vina predictions versus the experimental $K_i$ values of the 12 compounds. D) The possible binding modes of the top three compounds, in which the interacting residues are consistent with the previous studies. Overall, the wet-lab validations demonstrated that Dynaformer can effectively identify hit compounds with favorable binding affinities and novel scaffolds.