FlexSBDD: Structure-Based Drug Design with Flexible Protein Modeling

Authors: zaixi  zhang
mengdi  wang
qi  liu
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-02 03:12:25

Abstract: Structure-based drug design (SBDD), which aims to generate 3D ligand molecules binding to target proteins, is a fundamental task in drug discovery. Existing SBDD methods typically treat protein as rigid and neglect protein structural change when binding with ligand molecules, leading to a big gap with real-world scenarios and inferior generation qualities (e.g., many steric clashes). To bridge the gap, we propose FlexSBDD, a deep generative model capable of accurately modeling the flexible protein-ligand complex structure for ligand molecule generation. FlexSBDD adopts an efficient flow matching framework and leverages E(3)-equivariant network with scalar-vector dual representation to model dynamic structural changes. Moreover, novel data augmentation schemes based on structure relaxation/sidechain repacking are adopted to boost performance. Extensive experiments demonstrate that FlexSBDD achieves state-of-the-art performance in generating high-affinity molecules and effectively modeling the protein's conformation change to increase favorable protein-ligand interactions (e.g., Hydrogen bonds) and decrease steric clashes.

Main paper: 
 introduction: 
 Deep generative models are profoundly impacting drug discovery, particularly within the challenging subfield of structure-based drug design (SBDD) [1, 54, 36] . SBDD focuses on generating drug-like ligand molecules conditioned on target-binding proteins, necessitating precise modeling of complex geometric structures and detailed protein-ligand interactions. Some early attempts adopt autoregressive models to generate 3D ligand molecules atom-by-atom [53, 57] or fragment-by-fragment [85] . To overcome the limitations of autoregressive methods (e.g., error accumulation), recent works [28, 29] leverage non-autoregressive diffusion-based models [31] to predict the distribution of ligand atom types and positions via denoising and have achieved the state-of-the-art performance.Despite the remarkable success, most existing SBDD models treat target proteins as rigid and neglect the conformation change. However, according to the "induced fit" theory in biochemistry [43] , proteins are flexible structures that undergo structural changes upon ligand binding, leading to enhanced interactions and binding affinity. In structural biology, the protein structure that has a bound small molecule is referred to as ligand-bound or holo conformation, and the protein structure without a bound small molecule is called ligand-free or apo conformation [16] . For example, Figure . 1 (a)&(b) shows the aligned apo and holo-structures of two proteins, and the extent of structural change is influenced by the specific properties and structures of the proteins and ligands. The neglect of protein flexibility in SBDD leads to several significant drawbacks: (1) The generated protein-ligand complexes are prone to have sub-optimal protein structures and steric clashes as the protein cannot adaptively adjust structures according to different generated molecules [30] . (2) The chemical search space is overly restricted by existing holo data, limiting the exploration of diverse high-quality drugs. As a result, SBDD models tend to produce molecules that are similar to those fitting the predefined pocket space [53, 62] . (3) A large gap between real-world physical binding process and computational simulations, which may lead to high false positive rates in real-world drug discovery applications i.e., most computationally designed drugs do not have real therapeutic effects [2] . However, several challenges exist for flexible protein modeling in existing SBDD models. Firstly, proteins are macromolecules with thousands of atoms, which brings formidable high degrees of freedom for flexible structural modeling [35] . Secondly, there lack apo-holo structure pair datasets for learning structural changes [16] . Finally, the widely used diffusion-based models are time-consuming for exploring the huge space of flexible protein-ligand complexes [28] .To address the aforementioned challenges, we propose a new method FlexSBDD capable of modeling the flexibility of target protein while generating de novo 3D ligand molecules. To reduce the computational complexity, we focus on the key degrees of freedom in protein structure (i.e., C α coordinate, the orientation of the backbone frame, and the sidechain dihedral angels to determine the full atom structure) inspired by previous works [67, 77, 9] . As for the dataset, we employ the Apobind dataset [3] with additional Apo data generated by OpenMM relaxation [21] and Rosetta repacking [18] as data augmentation. For efficient and stable generation, we adopt a flow-matching framework [49, 5] that defines multimodal conditional flows for different components in protein-ligand complex and use E(3)-equivariant network with scalar-vector dual representation for learning the chemical and geometric information. The input to FlexSBDD is the initialized protein structure (e.g., apo structure). FlexSBDD learns to iteratively update protein-ligand structures and ligand atom types from time = 0 to 1 and finally outputs both the generated 3D ligand molecule and the updated protein structure (holo). Extensive evaluations on benchmarks and case studies show the advantage of FlexSBDD in generating structurally valid protein-ligand complexes with high affinity, more favorable non-covalent interactions, and fewer steric clashes. The code of the paper is provided at https://github.com/zaixizhang/FlexSBDD. We highlight our main contributions as follows:• We propose a flow-matching-based generative model FlexSBDD, capable of modeling protein flexibility while generating de novo 3D ligand molecules.• FlexSBDD not only achieves state-of-the-art performance on benchmark datasets (e.g., -9.12 Avg. Vina Dock score), but also learns to adjust the protein structure to increase favorable interactions (e.g., 1.96 Avg. Hydrogen bond acceptors) and decrease steric clashes.• With a concrete case study on KRAS G12C , a promising target of solid tumor, we demonstrate FlexSBDD's potential to discover cryptic pockets for drug discovery.2 Related Works 
structure-based drug design: 
 Structure-based drug design (SBDD) [87] aims to directly generate 3D ligand molecules inside target protein pockets. LiGAN [62] first uses 3D CNN to encode the protein-ligand structures and generate ligands by atom fitting and bond inference from the predicted atom densities. Several follow-up works have adopted autoregressive models for atom-wise [53, 50, 57, 79, 78] or fragment-wise [27, 60, 85, 78] generation of 3D molecules. For example, Pocket2Mol [57] adopts the geometric vector perceptrons [40] as the context encoder and autoregressively predicts the atom types, atom coordinates, and bond types until the generated molecule is completed. FLAG [85] and DrugGPS [82] predict the next molecular fragment and add it to the partially generated molecule in each round. Recently, powerful diffusion models have started to make a significant impact in SBDD, demonstrating promising results with non-autoregressive sampling [65, 28, 29, 46] . They usually represent the protein-ligand complex as 3D atom point sets, and define diffusion and denoising processes for both continuous atom coordinates and discrete atom types. For instance, TargetDiff [28] proposes a target-aware molecular diffusion process with a SE(3)-equivariant GNN denoiser. While there has been significant progress, existing approaches often neglect the critical aspect of protein flexibility. To address this gap, we explicitly model the conformational change of protein in FlexSBDD. Proteins are intrinsically dynamic entities and flexibility is a critical factor affecting protein's function and behavior in biological systems [72, 37] . The dynamic properties are traditionally modeled with physics-based methods such as Molecular Dynamics simulation (MD) [33] . To overcome the computationally demanding drawback of MD, some deep learning-based methods have been recently proposed, e.g., DynamicBind [51] , FlexPose [20] , SBAlign [68] , NeuralPlexer [61] , and DiffDock-Pocket [59] consider protein flexibility in protein-ligand docking. However, these methods can hardly extended to the challenging de novo ligand generation, leaving it an unsolved problem. 
preliminaries: 
 Notations and Problem Formulation: We model protein-ligand complex as C = {P, G} [84, 86, 76, 83] . The objective of FlexSBDD is to learn a conditional generative model p({P ′ , G}|P) that generates ligand molecule conditioned on the target protein and meanwhile updates the structure of the target protein (P ′ ). Proteins are composed of a sequence of residues (amino acids), each containing 4 backbone atoms (i.e., C α , N , C, O) and a sidechain that identifies the residue type. In this paper, the residue types are assumed known and the full atom structure of an amino acid can be represented by its C α coordinates x i ∈ R 3 , the frame orientation O (i) ∈ SO(3), and maximally 4 sidechain dihedral angles 4 , where i ∈ {1, • • • N p } and N p is the number of residues in a protein. The backbone structure of the residue can be determined according to their ideal local coordinates relative to the C α position x (i) and the orientation O (i) [23] . The sidechain conformations can be derived with the dihedral angles χ (i) as the bond length/angles are largely fixed [81] . With the above notations, a protein structure with N p residues can be compactly represented asχ (i) = {χ i1 , χ i2 , χ i3 , χ i4 } ∈ [0, 2π)P = {x (i) , O (i) , χ (i) } Np i=1 (see the illustration in Figure. 1(c)). Following previous works [74] , we treat each amino acid as a node and integrate backbone orientation and sidechain dihedral angles as node features. Such method enjoys the advantage of fewer nodes and less computational cost. The generated ligand molecule can be represented as a set of atoms:G = {a (i) , x (i) } N l i=1, where a (i) ∈ R na indicate the atom type (n a is the total number of atom types we consider) and x (i) ∈ R 3 denotes the atom coordinate. In the constructed protein-ligand complex 3D graph, we use x (i) to denote both C α coordinates and the ligand atom coordinates for conciseness.Riemanian Flow Matching: Flow Matching (FM) [49, 5, 4] , a simulation-free method for learning continuous normalizing flows (CNFs) [15] , has shown better performance and efficiency than diffusion-based models on a series of biomolecular tasks [9, 77, 69] . To model the complicated protein-ligand complex structures, we need to apply the general flow matching on the Riemannian manifolds [14] . Let M be the manifold space with metric g. q is the probability distribution of data x ∈ M, and p be the prior distribution. The time-dependant probability path on M is defined as p t∈[0,1] : M → R >0 satisfying p 0 = p and p 1 = q. u t (x) ∈ T x M is the corresponding gradient vector of the path on x at time t. Flow Matching aims to learn a neural network v θ (x, t) to approximate the target vector field u t :L F M (θ) = E t,pt(x) ∥v θ (x, t) -u t (x)∥ 2g . However, u t is intractable in practice and an alternative is to use a conditional density path p t (x|x 1 ) with a conditional gradient EQUATIONbecause L F M and L CF M have the same gradients according to previous works [14, 49] . t is sampled from the uniform distribution between 0 and 1. Once the gradient field v θ is learned, we can integrate ordinary differential equations (ODE): d dt ϕ t (x) = v θ (ϕ t (x), t) with ϕ 0 (x) = x, ϕ t (x) = x t and ODE solvers [10] to push data from prior distribution p 0 to the data distribution p 1 . Specifically, a data point from the prior distribution p 0 is C 0 = {P 0 , G 0 }, where P 0 is the initialized protein structure (e.g., the apo conformation) and G 0 is the initialized ligand molecule. The number of ligand atoms is sampled from the reference dataset distribution; the atom types are initialized with uniform distributions; the atom coordinates are initialized with Gaussian distributions inside the protein pocket following [28] . The target data distribution p 1 is the holo protein-ligand complex C 1 from datasets.Question: For each theoretical result, does the paper provide the full set of assumptions and a complete (and correct) proof? Guidelines:• The answer NA means that the authors have not reviewed the NeurIPS Code of Ethics.• If the authors answer No, they should explain the special circumstances that require a deviation from the Code of Ethics. • The authors should make sure to preserve anonymity (e.g., if there is a special consideration due to laws or regulations in their jurisdiction). 
flexsbdd: 
 Figure . 2 shows the overview of FlexSBDD. Given the initialized target protein structure, the goal of FlexSBDD is to generate the binding ligand molecule as well as the adjusted protein structure. In other words, the output is the generated protein-ligand complex. Given the complexity of the protein-ligand system, we first introduce flow matching on the protein backbone (Sec. 4.1), side chain (Sec. 4.2), and ligand atom type (Sec. 4.3) respectively. Then we show E(3)-equivariant network in Sec. 4.4. Finally, we discuss the training and generation procedure of FlexSBDD (Sec. 4.5).Following previous works [9, 41, 77] , The backbone atom positions of each residue in a protein can be modeled as a rigid frame T = (x, O) ∈ SE(3), consisting of C α coordinate x ∈ R 3 (we use x to denote coordinates in the rest of this paper) and the frame orientation matrix O ∈ SO(3). For simplicity, the following deduction focuses on a single residue/frame and can be generalized to all the residues in the protein. The conditional flow of frame T t is defined to be along the geodesic path connecting T 0 (apo frame) and T 1 (holo frame): T t = exp T0 (t log T0 (T 1 )), where exp T represents the exponential map and log T denotes the logarithmic map at T [77, 9] . Specifically, the conditional flow for the Euclidean coordinate vector x t and the orientation matrix O t are defined as:EQUATIONBoth R 3 and SO(3) are simple manifolds and their closed-form geodesics can be derived. The exponential map exp O0 can be computed using Rodrigues' formula and the logarithmic map log O0 is similarly easy to compute with its Lie algebra so(3) [77] . In protein-ligand complex, the ligand atom coordinates have the same data modality and probability path as C α . The loss function of FlexSBDD for protein backbone and ligand coordinates is the summation of the following two terms:EQUATIONEQUATIONwhere the superscript (i) in x (i) and O (i) is the index of the residues/ligand atoms. For conciseness, we use Equ. 4 to represent the coordinate loss with respect to N p α-Carbon and N l ligand atoms.In FlexSBDD, we define sidechain torsion angles on the torus space to model the protein sidechain structures. There are maximally four sidechain dihedral angels for each residue i.e., 4 and there are totally 4N p torsion angles [81] . Since each torsion angle lies in [0, 2π), the 4N p torsion angles of sidechains define a hypertorus T 4Np . The manifold of the hypertorus is parameterized as the quotient space R 4Np /2πZ 4Np , leading to the equivalence relations χ = (χ (1) , . . . , χ (4Np) ) ∼ (χ (1) + 2π, . . . , χ (4Np) ) ∼ (χ (1) , . . . , χ (4Np) + 2π) [39, 81] . We use the linear interpolation paths and the conditional flow for χ is defined as:χ (i) = {χ i1 , χ i2 , χ i3 , χ i4 } ∈ [0, 2π)EQUATIONwhere reg(•) means regularizing the torsion angles by reg(χ) = (χ + π) mod (2π) -π. This leads to the closed-from expression of the loss to train the conditional Torus Flow Matching:EQUATIONThe ligand atom types are denoted as a = {a (i) } Np i=1 where a (i) is the i-th atom probability vector with n a dimensions: a (i) ∈ R na . To build a path, we define a 0 as a uniform distribution over all atom types and a 1 as the one-hot vector indicating the ground truth atom type. The probability path is define as a t = ta 1 + (1 -t)a 0 , and u t (a|a 0 , a 1 ) = a 1 -a 0 . a t is a probability vector because its summation over all types equals 1. Following [47, 70, 12] , we use Cross-Entropy loss CE(•, •) to directly measure the difference between the ground truth type and the predicted one:EQUATIONWe also note the recent progress of sequential flow matching methods [71, 12] , which can be seamlessly integrated into FlexSBDD and are left for future works.FlexSBDD is parameterized with an E(3)-Equivariant Neural Network with scalar-vector dual feature representation to effectively capture the 3D geometric attributes [40, 57] . The scalar features contain basic biochemical knowledge (e.g., residue/atom types), and the vector features contain geometric knowledge of the structure (e.g., direction to the geometric center). The basic building blocks include geometric vector linear (GVL) and geometric vector perceptron (GVP). We also incorporate the geometric vector normalization (GVNorm) and the geometric vector gate (GVGate) for the model's stability and better performance. There are mainly two modules: an encoder that is responsible for encoding the protein-ligand complex 3D graph (see details in Sec. D.2) and a decoder that updates both the coordinates, frame orientation, atom types, and side-chain torsion angles (see details in Sec. D.3). Similar to previous works [57] , the update process satisfies the E(3)-equivariance. More model details are in the Appendix. D.5.Training with Data Augementation: For protein-ligand complexes from the training set, we associate them with apo conformations from Apobind [3] to create apo-holo pairs. Additionally, we create synthetic apo conformations as data augmentation. This is done by first removing the ligands from the holo proteins, followed by applying OpenMM [21] relaxation and Rosetta repacking [18] to these proteins. For each holo-structure, we generate a total of 9 additional structures: 3 only with sidechain repacking, 3 with both structure relaxation and repacking, and 3 with additional random perturbations with up to 30 degrees to the sidechain angles. In each training iteration, we randomly sample from the corresponding pool of apo structure C 0 of the holo-structure C 1 and interpolate to obtain C t . The overall loss function is the weighted summation of the above four loss functions:EQUATIONwhere w atom , w coord , w ori , and w sc are the loss weights and are set to 2.0, 1.0, 1.0, and 1.0 in the default setting. We adopt Adam [42] optimizer for the optimization and finish training on a Tesla A100 GPU.Generation: Starting with the apo structure and an initialized ligand molecule, denoted as C 0 , the generation process of FlexSBDD is the integration of the ODE dCt dt = v θ (C t , t) from t = 0 to t = 1 with an Euler solver [10] . Specifically, for each component in C t , we have:EQUATIONEQUATION)where ∆t is the time step. norm(•) means normalizing the vector to a probability vector such that its summation is 1, and reg(•) means regularizing the angles by reg(τ ) = (τ + π) mod (2π) -π. 
experimental setup: 
 Datasets: Following previous works [65, 57] , we use two popular benchmark datasets for experimental evaluations: CrossDocked and Binding MOAD. Binding MOAD dataset [34] contains around 41k experimentally determined protein-ligand complexes. We further filter and split the Binding MOAD dataset based on the proteins' enzyme commission number [7] , resulting in 40k protein-ligand pairs for training, 100 pairs for validation, and 100 pairs for testing following previous work [65] . CrossDocked dataset [25] contains 22.5 million protein-molecule pairs generated through cross-docking. We use the same data preprocessing and splitting as [52] , only high-quality docking poses (RMSD between the docked pose and the ground truth < 1Å) are kept and 30% sequence identity dataset split is adopted. This produces 100, 000 protein-ligand pairs for training and 100 proteins for testing. We regard the protein-ligand structures in the datasets as holo-structures. The corresponding apo structures are obtained from Apobind and the generated apo structures as described in Sec. 4.5. We note that it is fair to compare FlexSBDD with other baseline methods as the additional apo structures contain no ligand molecules and cannot be used by baselines for training.Baseline Methods: We compare FlexSBDD with five representative methods for SBDD. LiGAN [62] is a conditional VAE model that represents protein-ligand complex as an atomic density grid. AR [52] and Pocket2Mol [57] are autoregressive schemes that generate 3D ligand molecules atomby-atom. TargetDiff [28] and DecompDiff [29] are state-of-the-art diffusion-based models.Evaluation: We comprehensively evaluate the generated molecules from three perspectives: binding affinity and molecular properties, molecular structures, and protein-ligand interactions: (1) Following previous work [28, 29] , we use AutoDock Vina [22] to calculate and report the mean and median of binding affinity-related metrics, including Vina Score, Vina Min, Vina Dock, and High Affinity. Vina Score directly measures the binding affinity based on the generated 3D molecules; Vina Min performs a local structure relaxation before calculation; Vina Dock includes an extra step of re-docking, which serves to reveal the optimal binding affinity achievable; High affinity measures the percentage of generated molecules with higher binding affinity than the reference molecule in the (3) We adopt PoseCheck [30] to evaluate whether methods can establish favorable interaction between protein and ligand. The interactions include Hydrogen bonds and Hydrophobic interactions. We also conduct Steric clashes analysis to examine unphysical structures/interactions.In Table . 1 and 4, we show the binding affinity and the drug-related properties of the generated molecules on two benchmarks. We can observe that our FlexSBDD outperforms baselines by a large margin in affinity-related metrics. For example, FlexSBDD surpasses the strongest baseline DecompDiff by 0.73 and 0.82 in Avg. and Med. Vina Dock on CrossDokecd and 0.92 and 1.03 respectively on Binding MOAD. These gains indicate the strong capability of FlexSBDD to explore high-affinity drug molecules and adjust protein structures for tight binding. As for molecular properties QED and SA, FlexSBDD also achieves competitive performance with baselines. As discussed in [29] , these properties are usually employed as preliminary screening criteria in real drug discovery scenarios as long as they fall into a reasonable range. Finally, the high diversity indicates that FlexSBDD can explore larger chemical space with flexible protein modeling, which is important for early drug discovery. Generation efficiency is also a key factor to consider when sampling a large batch of molecules for screening. A major drawback of widely-used diffusion models is their inference speed, which may require 1000 time steps to produce high-quality samples. In contrast, flow matching methods remove stochasticity from the sampling path and can achieve stable and high-quality generation with much fewer steps (e.g., 20 steps in FlexSBDD). In Figure . 3, we observe that FlexSBDD can generate molecules much more efficiently than autoregressive-based methods such as AR [52] and diffusion-based methods such as TargetDiff [28] and DecompDiff [29] .We further consider steric clashes, hydrogen bonds, and hydrophobic interactions. Steric Clashes happens when two neutral atoms come into closer proximity than the combined extent of their van der Waals radii [63] , indicating energetically unfavorable and physically unrealistic structures. Hydrogen bonds (HBs) [58] and Hydrophobic interactions are polar interactions that significantly contribute to the binding affinity between proteins and ligands (More details in Appendix A.1). In Figure . 3, we show the average number of steric clashes, hydrogen bond donors, acceptors, and hydrophobic interactions in the generated ligands (without redocking). We observe that FlexSBDD can generate ligands introducing fewer clashes and more favorable interactions. For example, the average steric clashes for DecompDiff and FlexSBDD are 6.43 and 1.39 respectively. The average number of HB Acceptors for DecompDiff and FlexSBDD are 1.18 and 1.96 respectively. This could be attributed to the flexible protein adjustment capability of FlexSBDD, which could adaptively adjust protein and ligand conformations to reduce clashes and increase favorable protein-ligand interactions.In Figure . 4, we show examples of the generated ligand molecules for target proteins. Especially, We colored the original holo protein structure green and the updated structure with FlexSBDD cyan for comparison. Firstly, we observe FlexSBDD can generate ligand molecules with higher affinity and comparable QED and SA compared with reference complexes from datasets and molecules generated by DecompDiff. Moreover, the protein structures of FlexSBDD are adjusted to accommodate the generated ligand molecules. Consistent with the prior knowledge in biology [8] , we generally observe that the loop regions in protein structures exhibit greater flexibility, whereas the alpha-helix and beta-sheet regions display more rigidity. To further evaluate the validity of the updated protein structure, we employ self-consistency Template Modeling (scTM) following [48] (more details in Appendix A.2). scTM score ranges from 0 to 1 and a larger scTM score indicates better structural validity. On average, the updated protein structures by FlexSBDD have a scTM score of 0.964, comparable to the score of the original structures from the datasets (0.975). Moreover, to evaluate the validity of sidechain structure, we compute the Mean Absolute Error (MAE) of sidechain angles following previous works [81] . In Table 2 : Jensen-Shannon divergence between bond distance distributions of reference and generated ligands (the lower, the better). "-", "=", and ":" denote single, double, and aromatic bonds.We further conduct sub-structure analysis to evaluate whether FlexSBDD can generate valid molecular conformations. In Table . 2 and Table.5 in Appendix B, we compute different bond distance and bond angle distributions of the generated molecules and compare them against the corresponding reference empirical distributions following [28, 29] . We can observe that our model has a comparable or better performance on all the bond distances and angles, demonstrating the strong capability of FlexSBDD to generate realistic 3D molecules directly. The dynamic nature of proteins frequently results in the formation of cryptic pockets, which can reveal novel druggable sites not found in static structures and make previously "undruggable" proteins into potential drug targets [55] . To study the capability of FlexSBDD to explore cryptic pockets, we take KRAS G12C for a case study, which is a promising target in the treatment of solid tumors, and over 3 decades of efforts have been devoted to discovering its inhibitors (drug molecules) [17, 32] . The binding mode of ARS-1620 (green, PDB id 5V9U) represents the typical binding pocket exploited by previous research, which limits the exploration of high-affinity inhibitors. Here, we take the protein structure of ARS-1620 as the apo structure and generate ligand molecules. By comparing and filtering the generated molecules according to recent literature [45] , we managed to rediscover the cryptic pockets with FlexSBDD. In Figure . 5, the updated structure is colored cyan and the generated ligand molecule is colored orange. We observe that the side chain rotation of residue Histidine-95 (marked with the red rectangle) is consistent with the report in [45] (PDB id 6P8Z), which forms a new subpocket and contributes a lot to binding affinity. This case study demonstrates FlexSBDD's capability to accurately model flexible protein structure, update sidechains to reduce steric clashes, and explore cryptic pockets for drug discovery.We conduct a series of ablation experiments to study the effect of different modules on the generation capability of FlexSBDD: (1) Exp0: In model training, we remove the data augmentation mentioned in Sec. 5.1, (2) Exp1: we replace the geometric vector modules with EGNN [64] adopted in [65] , which only has scalar features without vector features, (3) Exp2: we do not update the backbone structure of the protein (i.e., x, O) (4) Exp3: we do not update the sidechain dihedral angles of the protein (i.e., χ), (5) Exp4: we fix the whole protein structure in ligand molecule generation. We retrain all the FlexSBDD variant models for comparison. The results are present in By comparing results from Exp0 and FlexSBDD, we can find that data augmentation indeed helps boost performance by introducing more diverse apo structures. In comparing Exp1 with FlexSBDD, it is obvious that scalar-vector dual feature representation can benefit ligand molecule generation by well capturing geometrical features. When comparing EXP2, 3, and 4 with FlexSBDD, we observe that the modeling of flexible protein structures including backbone and sidechains is important for FlexSBDD. Specifically, we find modeling the flexibility of sidechain angles is more important than backbone structure as the backbone is more rigid. For example, the average Vina Dock drops to -8.86 for Exp3 (FlexSBDD w/o flexible sidechain) while only drops to -9.05 for Exp2 (FlexSBDD w/o flexible backbone). According to [43] , the sidechains are critical to "induced fit", where they adjust positions to accommodate the ligand and enhance binding affinity. Overall, the FlexSBDD variants still demonstrate competitive performance, showing the advantage of flow-matching architecture.We investigate the influence of two important hyperparameters on the performance of FlexSBDD, the hidden dimension size and the total number of iteration steps T in flow matching. In Figure . 6, we observe the trend of generating higher-quality molecules with larger hidden dimension sizes and more iteration steps. In the default setting, we set the node scaler feature size to 256 and total iteration steps to 20 to achieve a balance between the computational complexity and the generation quality. 
conclusion: 
 In this paper, we propose FlexSBDD, a deep generative model capable of modeling the flexible protein structure for ligand molecule generation. FlexSBDD adopts a flow matching framework for efficient ligand generation and leverages E(3)-equivariant network with scalar-vector dual feature representation to effectively model dynamic structural changes. Extensive experiments show its state-of-the-art performance in generating high-affinity molecules with less steric clashes and more favorable interactions. Potential future works include leveraging FlexSBDD to discover more cryptic pockets and modeling other functional proteins such as antibodies, peptides, and enzymes. 
broader impacts: 
 Question: Does the paper discuss both potential positive societal impacts and negative societal impacts of the work performed?Answer: [Yes]Justification: The paper discussed potential societal impacts.Guidelines:• The answer NA means that there is no societal impact of the work performed.• If the authors answer NA or No, they should explain why their work has no societal impact or why the paper does not address societal impact. • Examples of negative societal impacts include potential malicious or unintended uses (e.g., disinformation, generating fake profiles, surveillance), fairness considerations (e.g., deployment of technologies that could make decisions that unfairly impact specific groups), privacy considerations, and security considerations.• The conference expects that many papers will be foundational research and not tied to particular applications, let alone deployments. However, if there is a direct path to any negative applications, the authors should point it out. For example, it is legitimate to point out that an improvement in the quality of generative models could be used to generate deepfakes for disinformation. On the other hand, it is not needed to point out that a generic algorithm for optimizing neural networks could enable people to train models that generate Deepfakes faster. • The authors should consider possible harms that could arise when the technology is being used as intended and functioning correctly, harms that could arise when the technology is being used as intended but gives incorrect results, and harms following from (intentional or unintentional) misuse of the technology. • If there are negative societal impacts, the authors could also discuss possible mitigation strategies (e.g., gated release of models, providing defenses in addition to attacks, mechanisms for monitoring misuse, mechanisms to monitor how a system learns from feedback over time, improving the efficiency and accessibility of ML). 
safeguards: 
 Question: Does the paper describe safeguards that have been put in place for responsible release of data or models that have a high risk for misuse (e.g., pretrained language models, image generators, or scraped datasets)?Answer: [NA]Guidelines:• The answer NA means that the paper poses no such risks.• Released models that have a high risk for misuse or dual-use should be released with necessary safeguards to allow for controlled use of the model, for example by requiring that users adhere to usage guidelines or restrictions to access the model or implementing safety filters. • Datasets that have been scraped from the Internet could pose safety risks. The authors should describe how they avoided releasing unsafe images. • We recognize that providing effective safeguards is challenging, and many papers do not require this, but we encourage authors to take this into account and make a best faith effort.12. Licenses for existing assets Question: Are the creators or original owners of assets (e.g., code, data, models), used in the paper, properly credited and are the license and terms of use explicitly mentioned and properly respected?Answer: [Yes] Justification: The existing assets are properly cited and credited.Guidelines:• The answer NA means that the paper does not use existing assets.• The authors should cite the original paper that produced the code package or dataset.• The authors should state which version of the asset is used and, if possible, include a URL. • The name of the license (e.g., CC-BY 4.0) should be included for each asset.• For scraped data from a particular source (e.g., website), the copyright and terms of service of that source should be provided. • If assets are released, the license, copyright information, and terms of use in the package should be provided. For popular datasets, paperswithcode.com/datasets has curated licenses for some datasets. Their licensing guide can help determine the license of a dataset. • For existing datasets that are re-packaged, both the original license and the license of the derived asset (if it has changed) should be provided. • If this information is not available online, the authors are encouraged to reach out to the asset's creators.13. New Assets Question: Are new assets introduced in the paper well documented and is the documentation provided alongside the assets? Answer: [NA] Guidelines:• The answer NA means that the paper does not release new assets.• Researchers should communicate the details of the dataset/code/model as part of their submissions via structured templates. This includes details about training, license, limitations, etc. • The paper should discuss whether and how consent was obtained from people whose asset is used. • We recognize that the procedures for this may vary significantly between institutions and locations, and we expect authors to adhere to the NeurIPS Code of Ethics and the guidelines for their institution. • For initial submissions, do not include any information that would break anonymity (if applicable), such as the institution conducting the review. 
 
 Appendix: We consider steric clashes, hydrogen bonds, and hydrophobic interactions in protein-ligand interaction analysis with PoseCheck [30] . Steric Clashes happens when two neutral atoms come into closer proximity than the combined extent of their van der Waals radii [63] , indicating energetically unfavorable and physically unrealistic structures. In PoseCheck, a clash is counted when the pairwise distance between a protein and ligand atom falls below the sum of their van der Waals radii, allowing a clash tolerance of 0.5 Å. Hydrogen bonds (HBs) represent a form of molecular interaction where a hydrogen atom, covalently linked to an element of high electronegativity like nitrogen, oxygen, or fluorine, engages with another electronegative atom [58] . These bonds are crucial in numerous protein-ligand interactions [13] and necessitate precise geometric alignments to form [11] . HBs are directional, bestowing distinct roles on the atoms involved: the hydrogen covalently bonded to the electronegative atom acts as a "donor", while the atom that receives the HB is known as an "acceptor". Hydrophobic interactions are a type of non-covalent bonding that occurs among hydrophobic molecules or moieties within an aqueous setting. Driven by water's propensity to hydrogen bond with itself, these interactions result in the segregation of non-polar entities, compelling them to cluster together away from the water-rich environment [56] . This behavior significantly contributes to the binding affinity between proteins and ligands.Following [48] , scTM takes a generated structure and feeds it into ProteinMPNN [19] , a state-ofthe-art structure-conditioned sequence generation method. With a sampling temperature of 0.1, we generate eight sequences per input structure and then use OmegaFold [75] to predict the structure of each putative sequence. We follow [48] that substitutes AlphaFold2 with OmegaFold for better sequence prediction performance. Finally, scTM is measured by computing the TM-score [80] , a metric of structural congruence of the OmegaFold-predicted structure and the original updated structure by our FlexSBDD. scTM scores range from 0 to 1, with higher numbers corresponding to the increased likelihood that an input structure is designable (higher structural validity). In FlexSBDD, the sidechain torsion angles are predicted and the sidechain conformations are derived based on the dihedral angles and the ideal bond length/angles. To evaluate the validity of sidechain structure, we compute the Mean Absolute Error (MAE) of sidechain angles (degrees) following previous works [81] in Table . 6. We also compare the results with NeuralPlexer [61] , the state-of-theart protein-ligand complex structure prediction. In the Table 7 : Effect of different modules on the generation performance of FlexSBDD. We show the average number of interactions here. We highlight the best two results with bold text and underlined text, respectively.Besides Vina Scores, we also try to leverage other docking methods to evaluate the binding affinity [24] . In Table . 8, we further leverage Glide [26] to evaluate the generated ligand molecules, which demonstrates superior ability in filtering active compounds. Specifically, we calculate the min-inplace GlideSP score following [24] , where the ligand structure undergoes force-field-based energy minimization within the receptor's field before scoring. In the table below, we observe that FlexSBDD can also achieve the best score on Glide, demonstrating its strong performance in generating proteinbinding molecules. To construct the protein-ligand KNN graph, we set k as 8 (each node is connected to its nearest 8 neighbors). In the default setting, we use a hidden size of 256, 128, 128, and 64 for the scalar features of nodes, scalar features of edges, vector features of nodes, and vector features of edges, respectively. The Encoder and Decoder have 6 layers respectively with the number of attention heads set as 4.The number of integration steps in flow matching is 20 for FlexSBDD. The hyperparameters for the loss function: w atom , w coord , w ori , and w sc are selected based on grid search ({0.5, 1.0, 2.0, 3.0}). w atom , w coord , w ori , and w sc are set to 2.0, 1.0, 1.0, and 1.0 in the default setting. To train FlexSBDD, we use the Adam [42] as our optimizer with a learning rate of 0.001, betas = (0.95, 0.999), and batch size 4 for 500k iterations. It takes around 36 hours on n one NVIDIA GeForce GTX A100 GPU to complete the training.We show the pseudo codes of FlexSBDD training and generation in Algorithm 1 and 2.A series of previous works on biochemistry tasks [40, 57, 79, 20] have shown that representing the nodes and edges in 3D graphs with scalar-vector dual features can greatly improve performance. In Sample C 0 from the corresponding apo structure pool (Apobind and data augmentation)x t = (1 -t)x 0 + tx 1 ; O t = exp O0 (t log O0 (O 1 )) χ t = (1 -t)χ 0 + t • reg(χ 1 -χ 0 ); a t = ta 1 + (1 -t)a 0 ; // Interpolation L ← w atom L atom + w coord L coord + w ori L ori + w sc L sc // calculate loss according to Equ. 9; θ ← Update(θ, ∇ θ L F M ); return v θAlgorithm 2: Generation algorithm of FlexSBDD Input: Total number of integration steps T, and trained model v θ steps ← 0, t ← 0, ∆t ← 1/T ; Initialize C 0 with the apo structure and sampled ligand atoms;while steps ≤ T -1 do x (i) t+∆t = x (i) t + v θ (x (i) t , t)∆t; O (i) t+∆t = O (i) t exp v θ (O (i) t , t)∆t ; χ (i) t+∆t = reg(χ (i) t + v θ (χ (i) t , t)∆t); a (i) t+∆t = norm a (i) t + v θ (a (i) t , t)∆t ; t ← t + ∆t ; return C 1FlexSBDD, all nodes and edges in the target protein P and the generated molecules G are assigned with both scalar and vector features to better capture the 3D geometric information. The scalar features contain basic biochemical knowledge (e.g., residue/atom types), and the vector features contain geometric knowledge of the structure (e.g., direction to the geometric center). In the rest of this paper, we use "•" and "→" overheads to explicitly indicate scalar features and vector features (e.g., v and ⃗ v).We adopt the geometric vector linear (GVL) and the geometric vector perceptron (GVP) as the main building blocks to enhance the information flows between the scalar features and the vector features and achieve E(3)-equivariance [40] . The details of GVP and GVL are shown in Appendix D.5. Briefly, they propagate the vector features into the scalar features by row-wise norm and propagate the scalar features to the vector features through gating. The GVP further applies extra non-linear transformations to both the scalar and vector features, following the output of GVL:EQUATIONwhere ( v, ⃗ v) could be any pair of scalar-vector features. Incorporating vector features is essential for our model's ability to directly and precisely update the positions of atoms and model the conformation change of the flexible protein. In our model, we also incorporate the geometric vector normalization (GVNorm) and the geometric vector gate (GVGate) for model's stability and better performance. Specifically, GVNorm combines the layer normalization [6] with the vector normalization [40] ; GVGate performs skip connection and fuses features from different blocks.In FlexSBDD, we represent the protein pocket-ligand complex as a k-nearest neighbor (KNN) graph in which nodes represent protein residues or ligand atoms and each node is connected to its k-nearest neighbors. The input scalar features of residues are onehot embeddings of residue types and the input scalar features of ligand atoms are initialized with a uniform distribution over all atom types. The input vector features of residues are computed based on the coordinates of backbone and sidechain atoms, while the vector features of ligand atoms are the Euclidean vectors pointing to the geometric center of the ligand molecule (see Appendix D.4). The scalar edge features are the radial basis function (RBF) distance encodings [66] and the vector edge features are the relative coordinates. In the flow matching model, to further incorporate time step information, we embed time with sinusoidal embedding [73] and concatenate it with the input scalar node features following [28] .Generally, the encoder of FlexSBDD follows the message-passing paradigm to update the features. Denote the feature of node i at the l-th layer as ( v(l) i , ⃗ v (l) i ) and the edge feature between node i and j as ( ė(l) ij , ⃗ e (l) ij ) (we use subscript here for the index of nodes instead of time step in the main paper). Each layer consists of a message-passing module M l and an update module U l :EQUATIONwhere we useM (l) v = ( ṁ(l) i , ⃗ m (l) i ) and M (l) e = ( ṁ(l) ij , ⃗ m (l) ij )to denote the calculated messages for node i and the edge between node i and j. The message-passing module is based on an attention mechanism [73] . The query ( qi , ⃗ q i ), key ( kj , ⃗ k j ), value ( uj , ⃗ u j ), and edge bias ( ḃij , ⃗ b ij ) are first calculated with GVLs (the layer superscripts are omitted here for simplicity):EQUATIONThen the attention weights for the scalar are computed as:EQUATIONSimilarly for the vector features:EQUATIONwhere ⊙ denotes the Hadamard product, 1 is the vector with all entries as 1, h s and h v denote the hidden dimension size of the scalar and vector features respectively. softmax j means perform softmax over the j index. âij and âij are the attention weights (scalar and vector channel) between node i and j. The message is obtained as:EQUATIONwhere the message features with respect to node i are obtained by applying GVL to the weighted summation of the neighboring value features. Finally, the update module U l is formulated as:EQUATIONwhere the GVGate fuses the information from the (l -1)-th layer and the calculated messages from the l-th layer. GVNorm is appended to normalize the features.The node/edge features of the decoder are initialized from the output of the encoder and are updated the same as the encoder. The decoder of FlexSBDD further updates the protein C α and ligand atom coordinates as follows:EQUATIONNeurIPS Paper Checklist Guidelines:• The answer NA means that the paper has no limitation while the answer No means that the paper has limitations, but those are not discussed in the paper. • The authors are encouraged to create a separate "Limitations" section in their paper.• The paper should point out any strong assumptions and how robust the results are to violations of these assumptions (e.g., independence assumptions, noiseless settings, model well-specification, asymptotic approximations only holding locally). The authors should reflect on how these assumptions might be violated in practice and what the implications would be. • The authors should reflect on the scope of the claims made, e.g., if the approach was only tested on a few datasets or with a few runs. In general, empirical results often depend on implicit assumptions, which should be articulated. • The authors should reflect on the factors that influence the performance of the approach.For example, a facial recognition algorithm may perform poorly when image resolution is low or images are taken in low lighting. Or a speech-to-text system might not be used reliably to provide closed captions for online lectures because it fails to handle technical jargon. • The authors should discuss the computational efficiency of the proposed algorithms and how they scale with dataset size. • If applicable, the authors should discuss possible limitations of their approach to address problems of privacy and fairness. • While the authors might fear that complete honesty about limitations might be used by reviewers as grounds for rejection, a worse outcome might be that reviewers discover limitations that aren't acknowledged in the paper. The authors should use their best judgment and recognize that individual actions in favor of transparency play an important role in developing norms that preserve the integrity of the community. Reviewers will be specifically instructed to not penalize honesty concerning limitations. 


Decision: This work explores the use of dynamic information in molecular design, introducing a method called FlexSBDD based on flow matching and utilizing E(3)-equivariant neural networks. Experimental results demonstrate that the model achieves state-of-the-art performance across multiple benchmarks. After a thorough rebuttal, the reviewers found this work to be universally innovative and timely, given the growing capabilities in exploring protein dynamics for drug discovery. Despite some minor concerns, this is a good paper that I recommend for acceptance to NeurIPS.

Reviews: review: paper_summary: paper_summary: In this paper, the authors focus on the flexible protein setting in the structure-based drug design task. They propose a method named FlexSBDD, which is based on flow matching and utilizes E(3)-equivariant neural networks. The experiments show the advantages of the proposed FlexSBDD.
strength_weakness: strength_weakness: 1. The method does not explicitly model the interaction between the ligand and the protein, especially the pocket. The authors might consider building an external interaction graph between the residues in the pocket and the atoms of the ligand.
2. I would like to see more focused results on the binding interface.
questions: questions: See the weaknesses. Additionally, I noticed a work [1] that is closely related to this paper, also employing flow matching. However, the authors have not cited or discussed the differences between their work and this one.


**Minor Concern:**

Line 240 "Table.": Since "Table" is written in full without abbreviation, there is no need to add a period after "Table". This issue occurs in multiple places.



**Reference:**

[1] Schneuing, Arne, et al. "Towards Structure-based Drug Design with Protein Flexibility." ICLR 2024 Workshop on Generative and Experimental Perspectives for Biomolecular Design.
limitations: limitations: The authors have discussed the limitations of FlexSBDD.
paper_summary: paper_summary: The author identified an important missing factor in current SBDD modeling, i.e. protein structural change upon binding, and proposed an E(3)-equivariant flow matching framework named FlexSBDD that jointly models protein flexibility and molecule generation. This paper augmented apo structures for each holo protein in the training set via structure relaxation, Rosetta repacking and random perturbations. FlexSBDD achieves SOTA binding affinities, QED and diversity benchmarked on CrossDocked2020 and Binding MOAD dataset, together with fewer clashes and more HB donors and acceptors.
strength_weakness: strength_weakness: - Typo: (Line 90-91) However, these methods can hardly [be] extended to the challenging de novo ligand generation, leaving it an unsolved problem.
- Typo: (Line 591) Protein Sturcture Analysis => Protein Structure Analysis
- Ablation studies suggest that the biggest performance gain come from data augmentation. However, I feel that augmentation shouldn't matter that much, since the whole pipeline only implicitly utilizes the protein structural change and the explicit outcome is the generated ligand itself. Could the authors explain the role of data augmentation in your method?
- Since apo (unbound) and holo (bound) state proteins are of great focus in this paper, it seems to me that more reasonable metrics for protein analysis would be based on some molecular dynamics, instead of scTM or something that are not guaranteed to output holo proteins for given ligand molecules.
questions: questions: - How would FlexSBDD behave if applied to an inpainting scenario (with ground truth holo-protein)?
- For a fair comparison, I would recommend the authors to try some flexible docking tools that also take the protein flexibility into account, and see how the performances of FlexSBDD and other baselines change.
- The authors could elaborate a bit more on the evaluation of protein structures. For example, why is FlexSBDD superior to SOTA protein-ligand complex structure prediction method? Under what setting are they being evaluated and compared?
limitations: limitations: N/A
paper_summary: paper_summary: In this research paper, the authors introduce FlexSBDD, a novel model that employs flow matching for the generation of flexible protein-based molecules. Initially, the model sample a noisy ligand based on an empirical distribution. Subsequently, it conducts flow matching, performing on both geometric characteristics and atomic types. Through comprehensive experimentation, FlexSBDD has demonstrated superior performance compared to previous methodologies, including TargetDiff, Pocket2Mol, and DecomDiff.
strength_weakness: strength_weakness: 1. As discussed in the Dynamic-Bind, many flexible aspects of protein residues involve changes in backbone atoms, such as transitions from DFG-in to DFG-out conformations. The current implementation of FlexSBDD does not demonstrate its capability for conformational search in more rigorous settings, such as when the apo pocket exhibits substantial dissimilarity from the holo pocket.
2. Although FlexSBDD represents a commendable effort, it fails to address a fundamental issue in flexible-pocket generation: evolving the appropriate atomic number through the generation process. FlexSBDD initializes ligands from an empirical distribution without considering pocket structures, resulting in a somewhat stochastic generation process: sampling a large number of ligands may artificially expand the pocket, leading to an increase in biased docking scores. The illustration provided with the 1a2g example supports this observation.
3. The model architecture of FlexSBDD, which includes flow matching on various geometries, has been previously implemented in other works such as PPFlow, diminishing the architectural novelty of FlexSBDD.
4. The illustrations in Figure 4, including the structures from 4yhj and 1fmc, display unusual bond topologies in molecules generated by FlexSBDD, such as a ring containing two double bonds. This issue may stem from an oversight in bond modeling within the architecture.
questions: questions: 1. In the evaluation using Vina dock metrics, it remains unclear which protein structures were utilized for benchmarking. Could you specify whether the benchmarks employed the initial holo structures or those updated by FlexSBDD?
2. Could you elaborate on the statement made in lines 216-217? “We note that it is fair to compare FlexSBDD with other baseline methods as the additional 217 apo structures contain no ligand molecules and cannot be used by baselines for training.” I do not fully understand the points here. 
3. The code provided with the submission appears to lack both training and inference components, which undermines the credibility of the reported results. Could you address this omission?
4. Regarding the prediction of side-chain conformations, it appears that the analysis is limited to the mean squared error (MSE) of chi angles without considering the orientation within the residue frame prediction. Could you discuss the rationale behind this methodological choice?
limitations: limitations: The limitations and broader impacts are well discussed in Appendix E of the main paper.
paper_summary: paper_summary: This paper presents a new deep generative model, FlexSBDD, which advances the field of structure-based drug design (SBDD) by accounting for the flexibility of proteins when generating 3D ligand molecules. This approach addresses the shortcomings of traditional SBDD methods that assume proteins are rigid, leading to less effective drug interactions. To obtain apo-holo data pairs, this paper utilizes Apobind to generate apo strutures from known holo structures. This paper adopts advanced flow matching to learn the apo-holo dynamics of protein and generate 3D molecules in surprising 20 steps. Experiments on CrossDocked and Bingding MOAD show that FlexSBDD can generate drug-like molecules with highest affinity compared to auto-regressive and diffusion-based baselines.
strength_weakness: strength_weakness: This paper focuses on protein flexibility for SBDD. However, the evaluation (except for a case study in 5.4) mainly follows normal SBDD, while the readers may concern more about: 1. Quality of Apobind-generated apo structures and model-generated holo structures, which is not fully discussed and evaluated, 2. Why protein flexibility can help improve affinity for fixed holo structures? Thus, the reviewer concerns the thoroughness and presentation of the experiments. And the reviewer encourages the authors to elaborate on further evaluation and more discussion specific to protein-flexible SBDD.
I will raise my score if more flexible SBDD evaluations and discussions are presented.
questions: questions: 1. The reviewer believes FlexSBDD may generate different holo structures than ground-truth (GT). 1. The generated molecules should have higher affinities on generated holos than GT holos, right? 2. How different are the generated and GT holos? 3. If the GT ligand is provided, can the model generate accurate GT holos?
2. In Table 2, FlexSBDD shows much better bond length distribution than other baselines. How is bond generated? And why is it so good? And why is bond angles (Table 5) not much better?
3. In Table 3, why data augmentation boosts so much (and the most) on performance? It's hard to understand the connection between apo structures and affinity in fixed holos.
4. Affinity and other metrics have strong correlation with atom numbers. Please include atom numbers in Table 1 for fairer comparison.
I will raise my score if the questions are properly answered.
limitations: limitations: see weakness above

Comments: title: Thanks for your response!
Thanks for your support!
Further Response to Reviewer 7REV
Further Response to Reviewer Y45d
Further Response to Reviewer jwrm
comment: Dear Reviewer,

Thanks for your valuable suggestions and support! We will include the above discussions in our revised paper. Thanks!

Bests,  
Authors
Thanks for your response. I will keep my original score.
Dear Reviewer,

Thanks for your support! We are glad that our response addressed all your concerns.

Bests,  
Authors
Thanks for the authors' detailed response. It addressed all my concerns. I have raised my score to 7 in hopes that this paper gets accepted.
**Comment 5**: In Table 3, why data augmentation boosts so much (and the most) on performance? It's hard to understand the connection between apo structures and affinity in fixed holos.

**Response 5**: Thanks for the insightful comment! Data augmentation plays an important role in FlexSBDD. As indicated in lines 186-195 of the submitted paper, we take apo-holo structure pairs for training. In each training iteration, we sample the apo structures $\mathcal{C}_0$ and holo-structures $\mathcal{C}_1$ and interpolate to obtain $\mathcal{C}_t$, i.e., FlexSBDD is supervised to learn the protein structural changes and the apo data play critical roles. However, existing apo-holo data pairs from Apobind are quite limiting (~10K data points) and directly training FlexSBDD on it leads to severe overfitting according to our experiments. The low-quality of the generated protein structure would directly lead to inferior Vina scores because of steric clashes and few protein-ligand interactions. **Therefore, we propose to use data augmentation to increase the training dataset size, cover more diverse apo-holo transition paths, and boost FlexSBDD’s generalization capability.** The ablation studies show that the data augmentation contributes a lot to FlexSBDD’s performance. 

We will include more detailed discussions and analysis of the role of data augmentation in our revised paper. 


**Comment 6**: Affinity and other metrics have strong correlation with atom numbers. Please include atom numbers in Table 1 for fairer comparison.

**Response 6**: Thanks for the constructive suggestion! In FlexSBDD, the number of ligand atoms is sampled from the reference dataset distribution. We report the average number of ligand atoms below. Generally, the average num of atoms of FlexSBDD is comparable to the reference and other baseline methods. DecompDiff has the most number of atoms. We will include these statistics into the final version.

| Methods    | Reference | LiGAN | AR | Pocket2Mol | TargetDiff | DecompDiff | FlexSBDD|
|------------|--------------------|------|------------------------|------|---------------|------|------|
| Avg. Num of Atoms  | 22.8     | 19.9| 17.7   | 24.2| 29.4     | 23.0 |
**Comment 5**: In the evaluation using Vina dock metrics, it remains unclear which protein structures were utilized for benchmarking. Could you specify whether the benchmarks employed the initial holo structures or those updated by FlexSBDD?

**Response 5**: For the Vina scores of FlexSBDD, we use the updated protein structure. As for the other baselines, we follow previous works to use the target structure from the test set for evaluation.

**Comment 6**: Could you elaborate on the statement made in lines 216-217? “We note that it is fair to compare FlexSBDD with other baseline methods as the additional 217 apo structures contain no ligand molecules and cannot be used by baselines for training.” I do not fully understand the points here.

**Response 6**: Thanks for the detailed question! In FlexSBDD, we associate holo structures from training dataset (CrossDocked and Binding MOAD) with apo conformations from Apobind [3] to create apo-holo pairs for training. We want to note that the additional Apobind dataset does not contain protein-ligand structures (i.e., only protein structures) and cannot be used by baseline methods for training. The Apobind dataset employed by FlexSBDD will not bring data leakage or additional advantage. Therefore, it is fair to compare FlexSBDD with baseline methods. We will make the statement clearer in our revised paper. 

**Comment 7**: The code provided with the submission appears to lack both training and inference components, which undermines the credibility of the reported results. Could you address this omission?

**Response 7**: Thanks for the valuable comment! We have uploaded the training and inference codes. We will open-source all the codes upon paper acceptance. 

**Comment 8**: Regarding the prediction of side-chain conformations, it appears that the analysis is limited to the mean squared error (MSE) of chi angles without considering the orientation within the residue frame prediction. Could you discuss the rationale behind this methodological choice?

**Response 8**: Thanks for the detailed comment! We use the mean squared error (MSE) of chi angles to evaluate the prediction of sidechain conformations following previous works [9-11]. Generally, lower MSE indicate more precise sidechain structure prediction. In table 6 of the paper, we can observe that FlexSBDD achieves better performance in generating valid sidechain structures. 

During rebuttal, we perform additional analysis of the side chain prediction. For example, we follow DynamicBind to conduct a comprehensive analysis six distinct conformational changes across the picosecond level to millisecond level (molecular dynamics), each exemplified by a case from PDBbind. In the following table, we report Δpocket RMSD (including side chain and backbone) of DynamicBind and FlexSBDD, which measures the relative decrease in pocket RMSD (crystal structure as reference) compared with the AlphaFold structures. **A negative Δpocket RMSD indicates that the predicted aligns more closely with the crystal structure compared with the AlphaFold prediction.** We observe that FlexSBDD achieves competitive performance that improves the AlphaFold prediction to have lower pocket RMSD, even though it is not specifically designed for dynamic docking.

| Methods    | 6QGF | 6PGO | 6N8X | 6UWV | 6ROT | 6S9X | 
|------------|--------------------|------|------------------------|------|---------------|------|
| DynamicBind  | -0.669     | -1.140| -2.297    | -0.465| -2.327      | -5.245 |
| FlexSBDD| -0.680    | -0.976| -1.159  | -0.504| -1.148   | -3.083 | 

[9] Zhang Y, Zhang Z, Zhong B, et al. Diffpack: A torsional diffusion model for autoregressive protein side-chain packing. Advances in Neural Information Processing Systems, 2023.  
[10] McPartlon M, Xu J. An end-to-end deep learning method for protein side-chain packing and inverse folding[J]. Proceedings of the National Academy of Sciences, 2023, 120(23): e2216438120.  
[11] Dong T, Yang Z, Zhou J, et al. Equivariant flexible modeling of the protein–ligand binding pose with geometric deep learning[J]. Journal of Chemical Theory and Computation, 2023, 19(22): 8446-8459.
**Comment 4**: For a fair comparison, I would recommend the authors to try some flexible docking tools that also take the protein flexibility into account, and see how the performances of FlexSBDD and other baselines change.

**Response 4**: Thanks for the valuable suggestion! In the following table, **we combine the top-2 SBDD baselines with DynamicBind [51], the flexible docking to compare with FlexSBDD.** Specifically, after generating the ligand molecules with TargetDiff/DecompDiff, we further dock the ligand to the target protein with DynamicBind. For fair comparison, we show the results of Vina Dock, QED, SA, and Diversity below. We observe that applying flexible docking as post-processing indeed improve the Vina Dock of baselines. FlexSBDD still achieves the best results as an end-to-end generative model for de novo ligand generation while adjusting protein structures. We will include these new results and discussions in our revised paper. 

| Methods    | Vina Dock (↓) Avg. | Med. | High Affinity (↑) Avg. | Med. | QED (↑) Avg. | Med. | SA (↑) Avg. | Med. | Diversity (↑) Avg. | Med. |
|------------|--------------------|------|------------------------|------|---------------|------|-------------|------|--------------------|------|
| TargetDiff  | -8.17              | -8.25| 62.3%                  | 63.0%| 0.48          | 0.48 | 0.58        | 0.58 | 0.72               | *0.71*  |
| DecompDiff  | *-8.89*            | *-8.97*| *69.1%*             | *74.5%*| 0.45         | 0.43 | 0.61        | 0.60 | 0.68               | 0.68  |
| **FlexSBDD**| **-9.12**          | **-9.25**| **78.5%**         | **84.2%**| **0.58**    | **0.59** | *0.69*    | *0.73* | **0.76**          | **0.75**  |

- **Bold**: Best results
- *Italic*: Second best 

[51] Wei Lu, Ji-Xian Zhang, Weifeng Huang, Ziqiao Zhang, Xiangyu Jia, Zhenyu Wang, Leilei Shi, Chengtao Li, Peter Wolynes, and Shuangjia Zheng. Dynamicbind: Predicting ligand-specific protein-ligand complex structure with a deep equivariant generative model. 2023.

**Comment 5**: The authors could elaborate a bit more on the evaluation of protein structures. For example, why is FlexSBDD superior to SOTA protein-ligand complex structure prediction method? Under what setting are they being evaluated and compared?

**Response 5**: Thanks for the valuable suggestion! In Appendix B.3, we evaluate the sidechain prediction of FlexSBDD and compare it with SOTA protein-ligand complex structure prediction method NeuralPlexer. **We use the inpaint setting for experiments: the NeuralPLexer model is asked to jointly predict the structure for a cropped spherical region within 6.0 Å of any ligand atom by inpainting all the amino acid and ligand atomic coordinates from scratch; FlexSBDD generate the ligand structure and updates the apo structure to holo.** To evaluate the validity of sidechain structure, we compute the Mean Absolute Error (MAE) of sidechain angles (degrees) of the inpainting region. In Table 6, we observed that FlexSBDD achieves better results on MAE. This could be attributed to the advance flow matching framework, scalar-vector dual representation architecture, and data augmentation strategy. We will include more details of the experimental settings and discussions in the revised version. 

**Comment 6**: Typos: (Line 90-91) However, these methods can hardly [be] extended to the challenging de novo ligand generation, leaving it an unsolved problem. (Line 591) Protein Sturcture Analysis => Protein Structure Analysis

**Response 6**: Thanks for the detailed comments! We have corrected our typos and will submit the updated paper in the final version.