Integrating Protein Dynamics into Structure-Based Drug Design via Full-Atom Stochastic Flows

Q3: "Although relevant mathematical models are provided, there is a lack of essential intuitive explanations."

A3: We have indeed provided intuitive explanations for each mathematical component of our work. Flow models naturally model the transition between two distributions. Therefore, our framework intuitively employs flow models to represent the transition of pocket structures from the apo to the holo state, alongside the ligand generation process. The detailed mathematics in Section 3 specifically outlines how we model changes in pocket conformation and ligand during this process. If there are any questions regarding the intuitive explanations of specific mathematical parts, we are more than willing to discuss them further.

Q4: "The analysis and treatment of different states (apo and holo) in the dataset are not sufficiently explored. Moreover, the specific reasoning behind the ratio of the various states (apo to holo) remains unclear."

A4: We specify the dataset curation process in Appendix A, including the data source (see Line 813-830), holo and apo definition and the relevant treatments (see Line 833-842, Line 880-904). We keep the top 10 clusters (or fewer if the total number of clusters are smaller than 10) in every clustered complex MD data as different holo conformations (see Line 840-841). That is to say, each complex has one apo pocket structure predicted by alphafold2 and no more than 10 holo structures extracted from MD simulation. We also present some important figures for dataset visualizations and analysis. Figure 9 shows the number of comformation in dataset, from which we make clustering and define the holo structures. Figure 11 shows the distribution of molecular properties for ligands in holo conformations. To compare the apo and holo pockets in terms of binding affinity with ligands, we present the distribution of vina score and vina min for apo and holo complexes in Figure 12.

Q5: "It would be beneficial to consider comparing a broader range of graph-based generative models to conduct a thorough evaluation of the model's performance."

A5: In our study, we have already included several representative graph-based generative models as baselines. These models span different categories, such as autoregressive models and diffusion models. It's important to note that these baselines are typically limited to rigid pockets. Our work, however, introduces a novel and practical approach by considering protein dynamics in structure-based drug design. This is where our significant contribution lies, as we provide a suitable algorithm and model tailored to this new setting.

Q6: "Reviews of computer science papers typically encourage the inclusion of anonymous code, accompanied by straightforward and easily testable data."

A6: Upon acceptance of this paper, we will open-source both the code and the curated dataset, and offer a user-friendly interface for ease of use. We are also open to discussing further implementation details if needed. Please refer to Q6 & A6 of Reviewer gwDc for specific details on hyperparameters and model architectures.

Q7: "What challenges does this method have in dealing with protein flexibility? How can we learn from the direction of the solution?"

A7: We recognize that data scarcity is a significant barrier to accurately modeling flexible proteins. Although the MISATO dataset marks crucial progress, obtaining stable, long-term MD simulation data for complexes remains challenging due to high computational costs. Furthermore, advocating for the public release of more high-quality complex data is essential to overcoming this challenge.

A current trend is to incorporate more prior knowledge rather than relying solely on data-driven approaches. For instance, we've introduced full-atom representation and interaction loss in protein dynamics modeling, which prevents the model from learning atom-level protein-ligand interactions too implicitly, potentially hindering its learning.

In the future, a promising direction will be how to better incoporate physical rules into protein modeling. For example, how protein force-fields can be integrated in the modeling process to help models understand protein dynamics. A physical-informed architecture may offer a more optimal solution to tackle this challenge.

Thank you for your feedback. Please see below for our responses to the comments.

Q1: "The results pertaining to the main proposed evaluation algorithm appear to align with previous studies. Several algorithms referenced are well established in the literature. The authors should: clearly specify which algorithms are original and unique to this study; explicitly indicate which algorithms are derived from existing works, rather than from the authors' own proofs."

A1: Our work introduces a novel setting in structure-based drug design (SBDD) where protein dynamics are considered, necessitating the creation of new evaluation algorithms.

The evaluation algorithms derived from existing works include:

• Ligand property evaluation: This includes metrics like QED, SA, Vina Score, Lipinski, logP, High Affinity, and Complete Rate. These are widely used in previous SBDD works, such as TargetDiff [1]. Please see Lines 476-484 for details and related references.
• RMSD: This metric is commonly used to measure the structural difference between two structures and is frequently employed in prior research, such as AlphaFold2 [2]. In our work, it measures how closely our generated pocket structures resemble the real holo states.
• Protein-ligand non-covalent interaction profiling: This is also utilized in previous works, e.g., [3,4].

The evaluation algorithms original and unique to this study are:

• Cover ratio: Based on RMSD, this metric demonstrates the diversity and detailed capability of our method to generate holo-like states. Details can be found on Lines 508-512.
• Pocket volume difference: We propose using pocket volume differences to evaluate how generated pocket structures resemble real holo states. While the evaluation perspective is novel, the calculation of pocket volume employs well-established tools like POVME 3 [5].
• Binding Affinity evaluation considering protein dynamics: Considering protein dynamics in SBDD, we extend this consideration to our evaluations. We assess all methods using flexible docking and predicted affinity based on DynamicBind [6], introducing a novel evaluation protocol. For more details, refer to Q7 & A7 for Reviewer FBxC.

References:

[1] Guan, J., Qian, W. W., Peng, X., Su, Y., Peng, J., & Ma, J. (2023). 3d equivariant diffusion for target-aware molecule generation and affinity prediction. ICLR 2023.

[2] Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., ... & Hassabis, D. (2021). Highly accurate protein structure prediction with AlphaFold. Nature, 596(7873), 583-589.

[3] Lee, J., Zhung, W., & Kim, W. Y. (2024). NCIDiff: Non-covalent Interaction-generative Diffusion Model for Improving Reliability of 3D Molecule Generation Inside Protein Pocket. arXiv preprint arXiv:2405.16861.

[4] Zhang, Z., Shen, W. X., Liu, Q., & Zitnik, M. (2024). Efficient generation of protein pockets with PocketGen. Nature Machine Intelligence, 1-14.

[5] Wagner, J. R., Sørensen, J., Hensley, N., Wong, C., Zhu, C., Perison, T., & Amaro, R. E. (2017). POVME 3.0: software for mapping binding pocket flexibility. Journal of chemical theory and computation, 13(9), 4584-4592.

[6] Lu, W., Zhang, J., Huang, W., Zhang, Z., Jia, X., Wang, Z., ... & Zheng, S. (2024). DynamicBind: Predicting ligand-specific protein-ligand complex structure with a deep equivariant generative model. Nature Communications, 15(1), 1071.

Q2: "The paper employs various professional terms and abbreviations; however, the backgrounds and specific definitions for these terms are not adequately clarified (e.g., the term, stochastic full-atom flow)."

A2: Thank you for highlighting this concern. The term "stochastic full-atom flow" combines three aspects: "stochastic," "full-atom," and "flow." Here’s what each term signifies:

• "Flow" indicates that we use a flow model as our generative model.
• "Full-atom" implies that our model explicitly captures atom-level protein-ligand interactions, as opposed to only modeling the protein pockets at the residue level.
• "Stochastic" refers to the stochastic differential equation (SDE) variant of DynamicFlow, where stochasticity is introduced to enhance robustness.

These meanings are intended to be intuitive and straightforward. However, as suggested, we will ensure they are more clearly defined in future versions of our paper. We have introduced biological terminologies in Section 1 (Introduction) and covered mathematical terms in Section 3.1 (Background and Preliminaries), where we provide sufficient explanations. If there are any other professional terms or abbreviations that need clarification, please let us know, and we will gladly offer more detailed explanations.

Thank you for your positive feedback. Please see below for our responses to the comments.

Q1: "This work relies on training datasets constructed from MD simulations, which may lead the model to learn the simulated physics of MD rather than capturing real-world structural distributions. While various 3D generative models utilized simulation datasets (e.g., CrossDocked2020), I think it would be a minor issue. For future work, authors might consider leveraging true experimental datasets PDBbind to enhance reliability."

A1: Thanks for your suggestion. Indeed, simulated data may induce bias. Nevertheless, MD simulation is usually more reliable than docking methods as used in CrossDocked2020, though it is more computationally expensive. And the MISATAO dataset is actually sourced from PDBBind. Besides, we have filtered the MD simulated data to further enhance reliability. Please refer to Appendix A for details of data processing. We would like to follow your suggestion to simultaneously leverage more experimental data to improve our work as a future work.

Q2: "Bias in the predicted pocket structure selection for analysis."

A2: Thanks for your suggestion. Our work aims at exploring more holo-like states given an initial conformation of pocket. Thus, we select pockets yielding the best results based on Vina score, with the intention of identifying more suitable holo-like states for ligand design.

Additionally, we calculated the volume for randomly selected pockets and compared the volume differences between apo states and our generated pockets (via DynamicFlow-ODE and DynamicFlow-SDE) against real holo states for each target. We have reported their mean and median in the table below. More details concerning the volume calculation and a specific example can be found in Appendix H of the revised version (marked in blue).

The results confirm that our methods successfully discover holo-like pockets. Furthermore, we assessed binding affinity via flexible docking on randomly selected and all generated ligands to thoroughly evaluate performance. Our methods outperformed all baselines under these conditions. For more details, please refer to Q7 & A7 for Reviewer FBxC.

Q3: "Questions about applicability in real-world applications."

A3: In real-world applications with only apo structures available, the pocket region can be effectively defined using a geometric center and radius. This method is efficient because the pocket areas in apo and holo states are generally similar, even if their conformations differ. Due to the lack of a clearly defined pocket boundary in receptor proteins, our model can leverage this ambiguity. Both ligand-centric and receptor-centric methods are expected to produce similarly defined pocket regions for SBDD tasks, allowing our model to adapt and apply these definitions to apo protein structures generated by tools like AlphaFold, unconstrained by the rigid boundaries present in holo structures.

In the few cases where protein conformational changes might be significantly large, experts can select or define the pocket in the apo state. This expert input can propose diverse pocket definitions, facilitating the design of potential ligands using our model.

Q4: "The proposed models and baselines are distribution learning-based models. Therefore, QED, SA, Lipinski, logP should be similar to Reference ligands. (No ↑ or ↓.)"

A4: Thanks for pointing this out. We agree that from a distribution learning standpoint, generated ligands with more similar statistics to reference ligands indicate a superior model. In our work, we followed conventions from other studies, like TargetDiff, by using "↑" or "↓" to denote preferences in drug design. We will include a note on this in the revised version. From the distribution learning perspective, our methods outperform others across nearly all metrics, as they most closely approximate the reference molecules in terms of property statistics.

Q5: "What is the generation time scale?"

A5: We benchmark the inference time of baselines and our methods for generating 10 ligand molecules given the same pocket on 1 Tesla V100-SXM2-32GB. The default number of function evaluations (NFE) is 1000 for TargetDiff and IPDiff and 100 for our method.

As the results show, our methods are capable of generating high-quality ligands while simultaneously modeling protein dynamics at a fast speed, demonstrating a significant advantage in computational efficiency.

Thank you for response. Most of my concern is addressed. I'll maintain my score (8).

Q3: "There is no ablation study showing the impact of these architectural changes (i.e., atom-level SE(3)-equivariant geometrical message-passing layers and residue-level Transformer layers) on model performance. Providing an ablation to assess their effects on sample quality (e.g., Vina score) or computational efficiency (e.g., FLOPs) would strengthen the evidence of these improvements."

A3: Based on your suggestion, we conducted ablation studies to evaluate the impact of different architectural components on model performance. We implemented a baseline denoted as "w/o residue-level Transformer", which uses only atom-level SE(3)-equivariant geometrical message-passing layers. In this setup, atom-level output features are aggregated into residue-level features without employing a residue-level Transformer for further extraction, and these aggregated features are used to predict the residue frames’ translation, rotation, and torsion angles.

Additionally, we developed a baseline referred to as "w/o atom-level EGNN", which transforms the atom-level protein-ligand complex graph into a heterogeneous graph, where each node represents either a residue (with C-alpha coordinates, rotation vectors, and torsion angles as input features) or a ligand atom. In this variant, since we do not explicitly reconstruct the full atom representation of the pocket, the atom interaction loss is not applied.

The results are shown in the following table:

The results indicate that our proposed architecture significantly enhances binding affinity and is vital for effectively modeling protein-ligand interactions and protein dynamics.

Both variants ("w/o residue-level Transformer" and "w/o atom-level EGNN") are more computationally efficient due to their reduced model sizes. However, despite using both residue-level and atom-level models, our method maintains acceptable inference speed because our flow model can generate high-quality ligand molecules in fewer steps. (Refer to Q5 & A5 for Reviewer JxBu for a comparison of inference time between the baselines and our methods.)

Q4: "How did the author select the 50 test pockets other than having no overlap with the training set? "

A4: We ensured no overlap by verifying that for each holo pocket in the test set, the PM-score against any holo pocket in the training set is less than 0.95. The PM-score quantifies binding-site similarity using structural descriptors like residue nature and interatomic distances, calculated via PocketMatch [1]. We plan to explore additional similarity measures and data splitting methods in future work.

References:

[1] Nagarajan, D., & Chandra, N. (2013, February). PocketMatch (version 2.0): A parallel algorithm for the detection of structural similarities between protein ligand binding-sites. In 2013 National Conference on Parallel Computing Technologies (PARCOMPTECH) (pp. 1-6). IEEE.

• Methodology:
  • Protein Modeling: FlexSBDD employs a residue-level model for the protein pocket, whereas we utilize both residue-level and atom-level models simultaneously, leveraging atom37 mapping. This approach allows us to more precisely capture protein-ligand interactions at the atomic level, enhancing the accuracy and detail of our modeling.
  • Ligand Modeling: On the ligand side, we construct flow models for atom positions, atom types, and bond types simultaneously in an end-to-end manner. In contrast, FlexSBDD does not incorporate bond modeling within its flow model, instead generating bonds through empirical post-processing rules. Our approach allows for more integrated and cohesive modeling of ligand structures.
  • Discrete Variable Modeling: FlexSBDD uses continuous vectors to represent discrete variables (i.e., atom types) and utilize standard flow matching for continuous variables, employing "norm" for self-normalization to mimic probabilities. This introduces a lack of rigor due to the inference gap created by "norm." Conversely, we apply rigorous discrete flow matching using continuous-time Markov chains (CTMC) to model both atom and bond types, ensuring a more precise and theoretically robust representation. For detailed mathematical insights, see Section 3.1 and Lines 291-309.
  • Torsion Angles: For torsion angles, both FlexSBDD and our approach employ flow matching on the manifold of hypertorus, originally proposed for full-atom peptide design by Li et al. [1]. However, given the amino acid sequence in SBDD, we can explicitly address cases where certain residues have side-chain torsion angles with -rotation symmetry (e.g., of ASP). This is a more natural choice than FlexSBDD's method, which overlooks symmetry-induced angle period differences. For more details, see Lines 270-288 and Appendix B.
  • SDE Variants: Both DynamicFlow-ODE (ours) and FlexSBDD use ODEs to model transitions between apo and holo states and the ligand generation process. However, we also introduce an SDE variant to enhance robustness, with experimental results demonstrating that the DynamicFlow-SDE variant outperforms the DynamicFlow-ODE. For more details, refer to Section 3.3.
  • Interaction Loss: FlexSBDD models predict the vector field directly, while our approach predicts "clean" samples and reparameterizes them into vector fields. This allows us to introduce an interaction loss focused on atom distances, enhancing the learning of protein-ligand interactions from ground-truth data. Our experiments show that this interaction loss improves the model's understanding of these interactions and enhances the binding affinity of generated ligands.
• Evaluation: FlexSBDD assesses generated small molecule ligands based on QED, SA, Binding Affinity (measured by Vina), and profiles of protein-ligand interaction. We evaluate baselines and our methods from these perspectives, and also add an evaluation of how similar the generated pocket structures are to actual holo states by comparing pocket volume and RMSD. For details, see Lines 508-514, Figure 5, and Figure 6.

These differences underline our unique approach to incorporating protein dynamics in SBDD.

References:

[1] Li, Jiahan, et al. "Full-Atom Peptide Design based on Multi-modal Flow Matching." ICML 2024.

Thank you for your feedback. Please see below for our responses to the comments.

Q1: "Why not start with CrossDock2020 and BindingMOAD before moving into new dataset?"

A1: CrossDocked2020 enhances PDBBind with docking and filtering, but MD simulations offer greater accuracy and provide multiple valid holo states rather than just one. BindingMOAD resembles PDBBind as it comprises protein-ligand crystal structures. The MISATO dataset, however, includes molecular dynamics simulations for approximately 20,000 experimental protein-ligand complexes and is slightly larger than BindingMOAD. Since part of our work is focused on discovering holo-like states in the context of structure-based drug design (SBDD), a dataset enriched with diverse holo states is preferred. Thus, starting with a dataset like MISATO aligns well with our objectives.

Q2: Comparison with FlexSBDD published on NeurIPS 2024.

A2: Thank you for bringing FlexSBDD to our attention. According to ICLR Review Guidelines (https://iclr.cc/Conferences/2025/ReviewerGuide), contemporaneous papers—those published within four months of our submission—need not be compared. FlexSBDD became publicly available on September 29, 2024, just days before our deadline on October 1, 2024, showing we are independent works. We acknowledge FlexSBDD and our work both integrate protein flexibility or dynamics into SBDD but present differences in various aspects. Although their code has not been open-sourced, we will cite and discuss their work in future revisions.

Key distinctions between our work and FlexSBDD include:

• Motivation: FlexSBDD primarily seeks to incorporate protein flexibility into SBDD for optimizing complex structures and ligands. However, it overlooks the role of thermodynamic fluctuations that govern protein flexibility and conformational shifts, leading to diverse conformations with different ligands. In contrast, our work delves into the physics underlying these dynamics. We illustrate this by examining the DFG-in and DFG-out states of Abl kinase, emphasizing our motivation to integrate comprehensive protein dynamics into SBDD, beyond merely addressing flexibility.

• Data: FlexSBDD derives most of its apo data by augmenting holo data through relaxation/sidechain repacking. In contrast, we use AlphaFold2 to predict our apo data, potentially resulting in greater conformational changes. Additionally, our holo states are diverse, providing multiple states for each protein-ligand pair through molecular dynamics simulations, enabling a more thorough exploration of pocket conformational changes. This aligns with our motivation.

(continued on the next session)

Q9: "The validity of the generated pockets is not assessed. For example, the AF PLDDT may be used."

A9: Assessing the validity of the generated pocket structures poses a significant challenge. We have evaluated the pocket volume distribution, and our findings indicate that the generated pocket structures exhibit volumes more akin to MD-simulated holo pockets compared to apo pockets. The AF pLDDT metric is unsuitable for our assessment for two primary reasons: (i) AF pLDDT provides a measure of per-residue local confidence specifically for AF-predicted structures and cannot be computed for structures not predicted by AF; and (ii) it does not account for the presence of ligands. We acknowledge the importance of validating generated pocket structures and propose developing a related benchmark as future work.

Q10: FREED and FREED++ are not cited.

A10: FREED [1] and FREED++ [2] are ligand-based drug design methods, differing from structure-based drug design approaches. They utilize fragment-based molecule generation models combined with reinforcement learning algorithms, leveraging desired properties as rewards for designing molecules. Importantly, the 3D structures of pockets are not inputs to these models. In our revision, we will cite these methods and discuss the distinctions between their approach and ours.

References:

[1] Yang, S., Hwang, D., Lee, S., Ryu, S., & Hwang, S. J. (2021). Hit and lead discovery with explorative rl and fragment-based molecule generation. Advances in Neural Information Processing Systems, 34, 7924-7936.

[2] Telepov, A., Tsypin, A., Khrabrov, K., Yakukhnov, S., Strashnov, P., Zhilyaev, P., ... & Kadurin, A. FREED++: Improving RL Agents for Fragment-Based Molecule Generation by Thorough Reproduction. Transactions on Machine Learning Research.

Q11: Typos: "Line 214: string repetition" and "line 848: mistake (are can)".

A11: Thanks for pointing these out. We have fixed the typos in the revision.

Q7: A valid way to compute different settings for baselines and DynamicFlow is to use MD trajectories-based methods, such as MMPBSA.

A7: We agree that incorporating molecular dynamics in evaluation is important. However, relying on MD trajectory-based methods such as MMGBSA or MMPBSA can be highly cumbersome. These approaches employ different solvation models and require extensive computational resources, with MD simulations often taking months to complete for thousands of systems. Therefore, as an alternative, we opted for a deep-learning-based approach for flexible docking and scoring to achieve reliable and scalable evaluation.

Specifically, for each generated ligand designed by the baselines and our methods, we employ DynamicBind [1], a geometric deep generative model tailored for "dynamic docking", to generate 10 protein-ligand complex structures. DynamicBind also includes a model that predicts an "affinity" score, which estimates the negative logarithm of the binding affinity in concentration units. We then calculate the weighted average of these predicted binding affinities to derive the final "affinity" score, where a higher "affinity" score indicates better binding potential.

For each target, we assess the affinity of a randomly selected generated ligand (Single), the highest affinity among 10 generated ligands (Best over 10), and the best affinity across all 100 generated ligands (Best over all). We report the mean, standard deviation, and median of these affinities across 50 targets. The results are summarized as follows:

The results demonstrate that our methods outperform all baseline models across all evaluation settings. This highlights the strength of our approach in designing ligands with high binding affinity.

References:

[1] Lu, W., Zhang, J., Huang, W., Zhang, Z., Jia, X., Wang, Z., ... & Zheng, S. (2024). DynamicBind: Predicting ligand-specific protein-ligand complex structure with a deep equivariant generative model. Nature Communications, 15(1), 1071.

Q8: "Please add the DynamicsFlow models to Table 2 for easier comparison." "Table 2 shows that the diffusion methods work better for the ligand prediction given the generated pocket than the proposed model itself."

A8: Table 2 presents the performance of baseline models using either apo states or the pocket structures generated by DynamicFlow as inputs. Specifically, the entry labeled "Pocket2Mol" represents Pocket2Mol with apo states as input, while "Pocket2Mol + Our Pocket" represents Pocket2Mol utilizing pocket structures generated by DynamicFlow. The results highlight that DynamicFlow effectively discovers more appropriate holo states, thereby facilitating the design of ligands with high binding affinity.

The results show that some baselines outperform our proposed model in ligand prediction when using the generated pockets, which aligns with our expectations. This is largely because the baselines have the advantage of more extensive training data, whereas our method is specifically designed to discover pocket structures that closely resemble actual holo states. Consequently, with an accurate holo-like state, baselines might surpass our model in performance. However, this does not undermine the utility of our model; rather, it underscores its capability to discover pocket structures that closely align with real-world scenarios.

Q5: "The paper lacks the baselines for joint pocket and ligand generation or the motivation for their absence."

A5: Thanks for highlighting the related works [1,2,3,4]. We have cited and discussed these in Appendix L in the newly-updated revision. Although these works pertain to protein-ligand complex modeling, they address distinct tasks.

Our focus is on structure-based drug design (SBDD) considering protein dynamics, where we start with the apo state (initial pocket structure) and aim to generate the holo state and binding ligands. In our case, detailed ligand information, including both topology graphs and 3D structures, is not provided.

[1] concentrates on pocket design where the topology graph and initial 3D structure of the ligand are provided and the goal is to generate a compatible pocket for binding.

[2] focuses on protein-ligand complex structure generation where the protein sequence and the topology graph (i.e., 2D graph) of the ligand molecule are provided and only their 3D structures need to be generated.

[3] focuses on pocket representation learning via pretraining on pseudo-ligand-pocket complexes instead of SBDD.

[4] represents a standard SBDD method with rigid-pocket input. Although molecular dynamics were mentioned in this work, they refer to dynamics induced by the forward process of the diffusion model. In our work, we have compared our methods with various similar baselines [5,6]. They are all diffusion-based SBDD methods with rigid-pocket input, with slight differences in models or algorithms.

References:

[1] Zhang, Z., Lu, Z., Zhongkai, H., Zitnik, M., & Liu, Q. (2023). Full-atom protein pocket design via iterative refinement. Advances in Neural Information Processing Systems, 36, 16816-16836.

[2] Gao, B., Jia, Y., Mo, Y., Ni, Y., Ma, W. Y., Ma, Z. M., & Lan, Y. Self-supervised Pocket Pretraining via Protein Fragment-Surroundings Alignment. In The Twelfth International Conference on Learning Representations.

[3] Nakata, S., Mori, Y., & Tanaka, S. (2023). End-to-end protein–ligand complex structure generation with diffusion-based generative models. BMC bioinformatics, 24(1), 233.

[4] Huang, L., Xu, T., Yu, Y., Zhao, P., Chen, X., Han, J., ... & Zhang, H. (2024). A dual diffusion model enables 3D molecule generation and lead optimization based on target pockets. Nature Communications, 15(1), 2657.

[5] Guan, J., Qian, W. W., Peng, X., Su, Y., Peng, J., & Ma, J. (2023). 3d equivariant diffusion for target-aware molecule generation and affinity prediction. ICLR 2023.

[6] Huang, Z., Yang, L., Zhou, X., Zhang, Z., Zhang, W., Zheng, X., ... & Yang, W. (2024). Protein-ligand interaction prior for binding-aware 3d molecule diffusion models. ICLR 2024.

Q6: "The evaluation pipeline seems unfair and may be misleading: baseline models are designed to work with holo-state pockets, whereas, in the experiment, they are utilized to generate ligands for apo-state pockets."

A6: Our work represents a pioneering effort to integrate protein dynamics into structure-based drug design (SBDD), introducing a novel experimental setting. Given its unique and innovative nature, achieving a completely fair comparison with existing baseline models is inherently challenging. Our experimental setup is intentionally crafted to simulate scenarios where complete holo-state structures may not be readily available. This highlights the necessity of developing solutions that can effectively generate ligands using apo-state pockets, thus addressing an often overlooked yet critical aspect of the drug design process.

Thank you for your detailed feedback. Please see below for our responses to the comments.

Q1: "How did you transition from 19437 proteins to 16972 complexes? How are ligands selected, and why was the original number of proteins reduced?"

The MISATO dataset is curated using 19,443 protein-ligand complexes data from PDBbind (release 2022). According to the author, structures from PDBbind were excluded whenever non-standard ligand atoms or inconsistencies in the protein starting structures were encountered, resulting in 16972 complexes for MD simulation. We further exclude complexes with oligopeptide ligands, resulting in 12,695 complexes for further processing (see L833). Detailed data processing procedures are provided in Appendix A.

Q2: "What does "100-frame" mean? Every 100th frame of the MD trajectory? Or 100 frames from each MD simulation? If the second is true, how were the 100 frames selected from the MD trajectory?"

A2: The MISATO dataset collects 100 snapshots for each protein-ligand complex from the 8 ns MD trajectory with systematic sampling. We then cluster the 100 snapshots for each complex using an RMSD threshold of 1 Å (see Line 893).

Q3: Do you filter the protein-ligand complexes depending on the average RMSD throughout the trajectory? Please explain Figure 10 b.

A3: Yes, we filter the clustered holo structures based on the average RMSD_Ligand of the MD trajectory, applying a 3 Å threshold (see Line 966-967). A large RMSD_Ligand suggests potential unreliability in the MD trajectory of the protein-ligand complex, rendering the data questionable. The RMSD_Ligand measures the root-mean-square deviation of the ligand after aligning the protein with its native structure.

Figure 10 b shows the change in the number of complexes along our data processing procedures. Stage A represents the original MISATO dataset with 16,972 complexes. At Stage B, We filter out complexes where ligands are peptides, resulting in 12,695 complexes (see Line 886-887). At Stage C, we align proteins in our data and AlphaFold Database by sequence and filter out the unsuccessful cases, resulting in 7,528 complexes (see Line 945-958 for details). At stage D, we remove the data with RMSD_Ligand smaller than 3 Å, resulting in 5,692 complexes (see Line 964-967 for details).

Q4: "The dataset is simulational, which makes it less representative than the PDB."

A4: The MISATO dataset originates from PDBbind, where the protein-ligand structures are experiment-based, as indicated in Line 877. While these structures provide a solid experimental foundation, it's important to recognize that protein-ligand complexes are inherently dynamic, and their holo states are not singular. Incorporating molecular dynamics (MD) simulations enhances the dataset by introducing additional dynamic conformational information, which aids the model in exploring a broader range of valid holo states. Furthermore, our processed dataset can be utilized for other significant tasks, such as conformational sampling of protein-ligand complexes, thereby making a valuable contribution to the research community.

Q6: Hyperparameter and other architectural details.

A6: The hyperparameters about the training loss are provided in Q5 & A5 and directly in Equation 18. We use AdamW [1] as the optimizer with learning rate 0.0002, beta1 0.95, and beta2 0.999. Gamma $\gamma$ controls the stochasticity of the stochastic flow (see Equations 19, 20, and 21). We use 2.0, 0.005, 1.0, 2.0 as the values of gamma $\gamma$ for residue frames' translation, rotation, torsion angles, and ligand atom positions.

Our model consists of an atom-level SE(3)-equivariant graph neural network and a residue-level Transformers. The number of total parameters is 15.9 M. The total estimated model parameter size is 63.401 MB.

We include more details about the model architecture as follows:

Q7: Details of each $\phi_i$ in Section 3.4 and how this parameterization relates to the final flow matching objective.

A7: Each $\phi$ in Section 3.4 is an SE(3)-equivariant graph neural network. As introduced in Section 3.4 (especially Line 407-413) and Figure 3, these SE(3)-equivariant graph neural networks are atom-level and their outputs are used to predict ligand atom types and positions (which are further used to compute the flow matching losses for ligand atom/bond types and ligand atom positions) and serve as inputs of the residue-level Transformer. The outputs of residue-level Transformer are used to compute flow matching losses for protein residue frames' translation, rotation, and torsion angles. We hope our explanation helps you understand how these EGNNs relate to the final flow matching objectives.

Q8: "How were the hidden states used to predict atom positions and atom/bond types? Similarly in residue-level transformer, L412, final updated frames were used as predictions - for what? It was again unclear how torsion angles were predicted based on final residue level hidden states."

A8: The final SE(3)-equivariant features of ligand atom are directly used as predicted ligand atom positions without any further transformation. The final SE(3)-invariant features of ligand atom/bond are used to predict atom/bond types by a linear layer. The related details of the residue-level Transformer can also be found in FrameDiff [2]. Specifically, the SE(3)-invariant hidden states, translation and rotation of residue frames are updated after each Transformer block. The translation and rotation of the final frame are viewed as the final prediction without any further transformation or post process. The hidden states are used to predict the torsion angles by a linear layer, where each torsion angle of each residue is a scalar value. These operations are very simple and common.

References:

[1] Loshchilov, I., 2017. Decoupled weight decay regularization. arXiv preprint arXiv:1711.05101.

[2] Yim, J., Trippe, B.L., De Bortoli, V., Mathieu, E., Doucet, A., Barzilay, R. and Jaakkola, T., 2023. SE (3) diffusion model with application to protein backbone generation. ICML 2023.

Thank you for your feedback. Please see below for our responses to the comments.

Q1: About reproducibility.

A1: To help researchers better understand our framework, we will open-source both the code and curated dataset upon acceptance of this paper and provide a user-friendly interface. Additionally, we will provide further details about the model architecture and hyperparameters, as discussed in Q6 & A6, to enhance understanding.

Q2: "In Table 2, where are the results for DynamicFlow?"

A2: . Table 2 shows the performance of the rigid-pocket SBDD methods with our refined pocket conformation (i.e., holo pocket structures generated by DynamicFlow). More specifically, the entry "TargetDiff" corresponds to TargetDiff with apo pockets as input, and the entry "TargetDiff + Our pocket" corresponds to TargetDiff with holo states generated by DynamicFlow as input. This experiment shows that the holo states discovered by DynamicFlow might serve as better inputs for the rigid-pocket SBDD methods and improve their performance when real holo pockets are not available. We will include more descriptions about this in the revision to enhance clarity.

Q3: "In Figure 1, what are protein and ligand embeddings, where are they computed in the proposed workflow and how are they being used in complex graph and ligand graph respectively?"

A3: Figure 1 shows different holo conformations of Abl kinase with corresponding binding ligands as an example to illustrate the motivation of our work. In Figure 3, protein and ligand embeddings are derived from the encodings of protein atom features and ligand atom and bond features, respectively, through an embedding layer (i.e., learnable linear transformation). The protein atom feature contains its atom37 representation and residue type. (Atom37 is an all-atom representation of proteins where each heavy atom corresponds to a given position in a 37-dimensional array. This mapping is non amino acid specific, but each slot corresponds to an atom of a given name. Note that atom37 is widely used in protein modeling [1]) We concatenate the one-hot encodings of these two features (whose dimensions are 37 and 20, respectively) to derive the protein atom encodings (whose dimension is 57). We use the one-hot encoding of the atom type as the ligand atom encoding. We only consider explicitly modeling "C, N, O, F, P, S, Cl, Br" in ligands, so the dimension is 8. For ligand bond types, we consider "non-bond, single, double, triple, aromatic", so the dimension of ligand bond encoding is 5. The protein and ligand atom features are used as the initial node features in the complex graph. And the ligand atom and bond features are used as the initial node and edge features, respectively, in the ligand graph. The above encodings are common in modeling proteins and small molecules. We will include more details about this in the revision to improve clarity.

References:

[1] Jumper, J., Evans, R., Pritzel, A., Green, T., Figurnov, M., Ronneberger, O., Tunyasuvunakool, K., Bates, R., Žídek, A., Potapenko, A. and Bridgland, A., 2021. Highly accurate protein structure prediction with AlphaFold. Nature, 596(7873), pp.583-589.

Q4: "Given the number of components in the workflow, I would suggest including an aggregated workflow figure with end to end pipeline starting from the apo state input to the molecule and holo state output, complete with the final flow matching loss."

A4: Thanks for your suggestion. We will include a more comprehensive illustration of the overall workflow to promote understanding in the revision.

Q5: "The overall loss for this work is unclear, while individual losses for structural features for protein and ligand are provided, how are they aggregated is not mentioned."

A5: There are 7 individual losses: 4 continuous flow matching losses for residue frames' translation (Equation 5), rotation (Equation 7), torsion angles (Equation 8), and ligand atom position (same as Equation 5), 2 discrete flow matching losses for ligand atom and bond types (Equation 14), and interaction loss (Equation 18). They are first averaged across all residues or atoms in a training sample and then simply weighted summed with weights: 2.0, 1.0, 1.0, 4.0, 1.0, 1.0, 0.5.

Thank you for taking the time and effort to evaluate our submission!

We have just updated the revision again to include a newly added comprehensive illustration of the overall workflow, as you suggested. Please refer to Appendix D and Figure 14 for further details.