DockedAC: Empowering Deep Learning Models With 3D Protein-ligand Data For Activity Cliff Analysis

Thank you for your detailed feedback and insightful suggestion regarding the structural similarity assessment in the activity cliff definition. We will address your comments in the responses below.

[W1. Maximum Common Substructure]

We appreciate the suggestion to incorporate explicit structural comparison methods like Maximum Common Substructure (MCS), as it could provide valuable insights.

While computational complexity is a consideration, we have computed the number of shared atoms between the newly generated AC pairs, as shown in the newly added Appendix Figure 14 for one of the targets. The average proportion of the shared atoms is 86.78%, which shows that the identified AC pairs have a high similarity in terms of common substructures.

We will include this analysis in the revised manuscript. This addition will strengthen the rigour of our generated AC pairs and provide deeper insights into the structural aspects of activity cliffs.

[W2. multiple local minimum]

Thank you for highlighting this important question. We agree that considering multiple local minima would provide a more comprehensive view of potential binding modes.

To further enrich our dataset, we will generate top-N scoring poses (N=5-10) and calculate RMSD between different poses. We will include multiple poses in the dataset and provide pose scoring information.

The reason for that we currently only use the top-1 binding mode in the dataset is because, there are no good deep learning models designed to process multiple input conformations. This would be a very important direction for future work, as the binding affinity is the result of a binding/unbinding dynamic process. We will add this discussion to the revised paper.

[W3. analyses about dataset]

Q2: In Section 5.1, when correlating RMSE with RMSE-cliff, is there a potential for bias introduced by comparing AC data with the entire dataset (which includes both AC and non-AC data)?

You are right that RMSE includes both AC and non-AC data and RMSE-cliff is calculated from only the AC data.

We compare RMSE-cliff with RMSE to show that the prediction with AC data is a more difficult task, which can be seen from the result that most of the targets are lying above the line RMSE = RMSE-cliff (in Figure 4). It means that the average RMSE on the AC samples is higher than the average RMSE on all the samples, i.e.,

RMSE = \frac{\\# AC * RMSE_{cliff} + \\# {nonAC} * RMSE_{nonAC}}{\\#AC+ \\#nonAC} < RMSE_{cliff}.

\\# From the above, we can see that $RMSE_{cliff} > RMSE_{nonAC}$. Please let us know if you have any further questions.

Q3: Conceptually, using 3D information should improve AC prediction. However, deriving this conclusion in section 5.1 from the results comparing IGN[2] and other 2D models seems to lack depth in interpretation. It would be better to compare IGN with a 2D version of IGN to minimize bias.

The inputs of the 2D version and the 3D version are not the same. The 3D model has the pocket and ligand 3D structure as additional inputs. So it is not very straightforward to compare the two cases. In our paper, we do such a comparison in two different aspects. The first is to compare IGN or SS-GNN to other 2D graph neural network models, as shown in section 5.1. We chose IGN and SS-GNN because they are based on similar graph neural network architectures. For example, SS-GNN uses GIN as the backbone and uses the 3D information to enhance the edge features in the GIN. In the second aspect, we use a 3D model as a feature extractor and combine the 3D features with the handcrafted 2D features and employ MLP for affinity prediction, as shown in the Appendix Figure 8 and Table 6. The results from the second aspect also show that the features from 3D data can enhance the AC prediction.

Q4: In Section 5.3, is the sample size too small for calculating the p-value? Additionally, in Table 3, the percentage of AC is at most 0.6%. Why did the authors choose to show the correlation with a maximum AC data percentage of only 0.6%?

Thanks for your careful examination of our paper. 0.6 should be the proportion instead of the percentage. So the percentage is 60%. We have corrected the caption and figure label.

Thank you for your insightful feedback. Your suggestions will help strengthen our benchmark's utility and comprehensiveness.

Issue 1. [Regarding the evaluation of existing 3D binding affinity prediction models]

We appreciate your feedback regarding the evaluation scope of our benchmark. Our initial inclusion of IGN and SS-GNN was motivated by their graph-based architectures, enabling direct comparison with 2D GNN models (GCN, GIN) to demonstrate the impact of 3D structural information.

Following your suggestions, we have now extended the evaluation of PIGNet [1], RTMScore [2], TANKBind [3], DSMBind [4], KarmaDock [5] and EquiScore [6] on our collected dataset and the averaged results across the targets are in the following table (measured by Pearson correlation). We also include the full results of each target in the revised manuscript (Appendix Figure 15).

In general, the performance of these methods is not as good as that of the PDBbind test set. In the targets that have many homologous proteins in the PDBBind training set, the model tends to work better. For example, CHEMBL2147_Ki (DSMBind Pearson=0.688, pdb id 2j2i) has 103 homologous samples in the PDBBind dataset. CHEMBL2971_Ki (DSMBind Pearson=0.671, pdb id 4jia) has 61 homologous samples. While the proteins of CHEMBL219_Ki, CHEMBL228_Ki and CHEMBL233_Ki can not be found in the PDBBind dataset, DSMBind Pearson=-0.021, -0.087 and 0.033 respectively.

In addition, one of the key advantages of 3D approaches is their potential for cross-target applicability. To explore this, we retrain one of the 3D methods (IGN) on multiple targets to assess the cross-target performance. Though we cannot get the full results due to the limited time, from the current results, there are two observations. First, when trained with multiple targets, the performance is comparable to training from scratch on the specific target. Second, when the protein type (Protease) of the testing targets is not in the training targets, the performance drops from 0.9 to 1.4. This is consistent with previous studies [7][8]. When the targets are split by the protein sequence identity of 30%, the RMSE increases from 1.27 (random split) to 1.429 (30% identity split) (See Appendix Table 7 in the revised manuscript).

[1] Moon, Seokhyun, et al. "PIGNet: a physics-informed deep learning model toward generalized drug–target interaction predictions." Chemical Science 13.13 (2022): 3661-3673.

[2] Shen, Chao, et al. "Boosting protein–ligand binding pose prediction and virtual screening based on residue–atom distance likelihood potential and graph transformer." Journal of Medicinal Chemistry 65.15 (2022): 10691-10706.

[3] Lu, Wei, et al. "Tankbind: Trigonometry-aware neural networks for drug-protein binding structure prediction." Advances in neural information processing systems 35 (2022): 7236-7249.

[4] Jin, Wengong, Caroline Uhler, and Nir Hacohen. "SE (3) denoising score matching for unsupervised binding energy prediction and nanobody design." NeurIPS 2023 Generative AI and Biology (GenBio) Workshop.

[5] Zhang, Xujun, et al. "Efficient and accurate large library ligand docking with KarmaDock." Nature Computational Science 3.9 (2023): 789-804.

[6] Cao, Duanhua, et al. "Generic protein–ligand interaction scoring by integrating physical prior knowledge and data augmentation modelling." Nature Machine Intelligence (2024): 1-13.

[7] SS-GNN: A Simple-Structured Graph Neural Network for Affinity Prediction. ACS Omega. 2023, 8(25): 22496–22507.

[8] Pre-Training of Equivariant Graph Matching Networks with Conformation Flexibility for Drug Binding. Advanced Science, 2022, 9(33): 2203796

Thank you for detailed response. I'll increase the score.

Thank you for your careful review and insightful suggestions. We will try to address each point to clarify and improve our manuscript.

Weakness

The authors need to clarify the splitting of their dataset to ensure proper benchmarking

Question: Regarding line 234, could the authors provide more detail about why they split the dataset this way?

Data splitting is an important point for proper benchmarking. In the main paper, we train on one protein at a time, and benchmark the performance for each target. The reason is to compare with the 2D methods, which do not have protein information and can not be trained across proteins.

Besides performance improvement, another advantage of using 3D information is the potential for cross-target applicability. To explore this question and enrich our benchmark results, we combined the targets of the Ki as labels together and trained the IGN model again.

We use the following two settings:

• out-of-distribution (OOD): excluding one type of target for training (Protease)
• in-domain: using all the targets of Ki for training

We show the results of testing on the four Protease targets (Appendix Table 7). It is easy to see that when the protein type (Protease) of the testing targets is not in the training targets, the performance (avg. RMSEcliff) drops from 0.9 to 1.4. The results of testing on all the targets are also reported (Appendix Figure 13). From the current results, when trained with multiple targets, the performance is comparable to training from scratch on a specific target. The current results indicate that there is still a long way to go to fully exploit the multi-target 3D data and make the model generalize to new targets.

We have made it clear about the data splitting and model training, and added the new results to the maniscript.

The authors' chosen 3D GNNs (or 3D/equivariant transformers, which in my view may be most interesting to explore for this problem) are not up-to-date. For example, new models such as the Equiformer v2 architecture [1] could be included as well.

We fully agree that incorporating more new models would significantly enhance the benchmark results.

The main reason for using IGN and SS-GNN for benchmarking is that these two methods are based on graph neural networks and their results are more comparable to 2D GNN models like GCN or GIN.

The Equiformer v2 architecture [1] is for property prediction of a single molecule. It is not straightforward to extend it for affinity prediction from a pocket-ligand complex. But following your advice and suggestions from other reviewers, we evaluate the performance of several other more recent models for binding affinity prediction, including PIGNet [2], RTMScore [3], TANKBind [4], DSMBind [5], KarmaDock [6] and EquiScore [7].

All these methods have model weights trained with the PDBbind dataset and we test them on our collected dataset. The averaged results across the targets are in the following table (measured by Pearson correlation). The full results across 52 targets are presented in the revised manuscript (Appendix Figure 15).

Due to the time limit, we re-train the model EquiScore [7] with our dataset and compare with IGN and SS-GNN. The performance of EquiScore is slightly better than SS-GNN, but worse than IGN.

[1] Liao, Yi-Lun, et al. "EquiformerV2: Improved Equivariant Transformer for Scaling to Higher-Degree Representations." The Twelfth International Conference on Learning Representations.

[2] Moon, Seokhyun, et al. "PIGNet: a physics-informed deep learning model toward generalized drug–target interaction predictions." Chemical Science 13.13 (2022): 3661-3673.

[3] Shen, Chao, et al. "Boosting protein–ligand binding pose prediction and virtual screening based on residue–atom distance likelihood potential and graph transformer." Journal of Medicinal Chemistry 65.15 (2022): 10691-10706.

[4] Lu, Wei, et al. "Tankbind: Trigonometry-aware neural networks for drug-protein binding structure prediction." Advances in neural information processing systems 35 (2022): 7236-7249.

[5] Jin, Wengong, Caroline Uhler, and Nir Hacohen. "SE (3) denoising score matching for unsupervised binding energy prediction and nanobody design." NeurIPS 2023 Generative AI and Biology (GenBio) Workshop.

[6] Zhang, Xujun, et al. "Efficient and accurate large library ligand docking with KarmaDock." Nature Computational Science 3.9 (2023): 789-804.

[7] Cao, Duanhua, et al. "Generic protein–ligand interaction scoring by integrating physical prior knowledge and data augmentation modelling." Nature Machine Intelligence (2024): 1-13.

Thank you for your thoughtful review and valuable suggestions. We will address each of your points below.

For ECFP, what proportion of samples have an ECFP similarity > 0.9, given that similarities > 0.5 are already quite rare?

Thanks for this insightful question. The proportion of pairs with an ECFP similarity > 0.9 is indeed very small, less than 1e-4 out of all the pairs. However, when we annotate the samples appearing in these pairs as AC samples, the proportion of the AC samples represent a significantly larger fraction of all samples, approximately 10%.

Overall, when considering all three similarity metrics, the proportion of the AC samples ranges from 15.7% to 43.2% in our dataset. We will include specific statistics on the proportion of samples with similarity > 0.9. This will help better characterize our dataset's similarity landscape.

In pairs of activity cliffs with docking conformations, what noticeable differences exist in their 3D structures? It has not been verified whether docking conformations can accurately capture the true differences in interaction patterns.

We appreciate your thoughtful comments regarding the reliability of docking conformations in the context of activity cliffs. These are important considerations that merit careful discussion.

To validate our docking protocol, we have compared the docking poses of compounds with known crystal structures in [1]. The observed root-mean-square deviation (RMSD) values (median = 2.55Å) fall well within accepted standards for reliable pose prediction. Building on this foundation, we are implementing template docking strategies for cases where input ligands exhibit high structural similarity to reference ligands, which will further enhance pose accuracy.

In the manuscript, we also have discussed the limitations of docking, including potential inaccuracies in conformational sampling and scoring functions. A possible way to make the docking pose better is to run a short molecular dynamics simulation. This would be an important future work for us to enhance the dataset.

Can generated 3D docking conformations, which may differ from the crystal structures, reliably reflect the actual interaction differences between activity cliffs?

We have examined the structural differences in 3D conformations between some activity cliff pairs in [1]. As we have identified the binding site and use a relatively large exhaustiveness parameter for docking, the docking poses for these samples are quite accurate.

We acknowledge your concern about the docking conformations. We have shown that for 3D deep learning models, using the docked strucutral information is able to improve the performance. And as shown in [2], a deep learning model trained on crystal structures can still be used to make predictions with the docked complex. Its performance is quite close to the crystal structure setting because docking error is relatively low.

[1] Husby, Jarmila, et al. "Structure-based predictions of activity cliffs." Journal of chemical information and modeling 55.5 (2015): 1062-1076.

[2] Jin, Wengong, Caroline Uhler, and Nir Hacohen. "SE (3) denoising score matching for unsupervised binding energy prediction and nanobody design." NeurIPS 2024.