DeltaDock: A Unified Framework for Accurate, Efficient, and Physically Reliable Molecular Docking

The initial step of this module involves using RDKit [27] to generate a 3D conformer of the input ligand, as depicted in Fig. 1. Binding site prediction models including P2Rank and DSDP are adopted to extract druggable binding sites, and the binding sites predicted by these different methods are combined to form a set of candidate binding sites, denoted as $S=\{\varsigma_{1},\varsigma_{2},...\}$, where $\varsigma_{i}$ represents the geometric center of $i$-th binding site. For CPLA, the protein pocket $\rho_{i}$ is defined as the residues within 15.0 Å  to $\varsigma_{i}$.

To map the ligand and pockets into the embedding space, the ligand encoder Attentive-FP (AFP) [28] and protein encoder Geometric Vector Perceptron (GVP) [29] are employed. These encoders first extract informative ligand node and protein node representations, and the feature extraction process can be formally expressed as:

$\displaystyle H^{\mathcal{L}}=AFP(\mathcal{G}^{\mathcal{L}}),\;H^{\mathcal{P}*}=GVP(\mathcal{G}^{\mathcal{P}*}),$

(1)

where $H^{\mathcal{L}}$ is the ligand embedding matrix of shape $|\mathcal{V}^{\mathcal{L}}|\times d$ and $H^{\mathcal{P}*}$ is the protein residue embedding matrix of shape $|\mathcal{V}^{\mathcal{P}*}|\times d$. The ligand representations $m^{\mathcal{L}}$ and pocket representations $m^{\rho}_{i}$ are then obtained by pooling ligand nodes embedding and pocket nodes embedding:

$\displaystyle m^{\mathcal{L}}=Sum(H^{\mathcal{L}},\mathcal{V}^{\mathcal{L}}),\;m^{\rho}_{i}=Sum(H^{\mathcal{P}*},\mathcal{V}^{\rho}_{i}),$

(2)

where $\mathcal{V}^{\rho}_{i}$ is the protein node set of $i$-th pocket $\rho$, and the pooling operation is sum pooling. For the pocket encoder, we input the entire protein residue graph $\mathcal{G}^{\mathcal{P}*}$ rather than just the protein pocket residue graph, to incorporate global protein information into the pocket representation.

With ligand representation $m^{\mathcal{L}}$ and pocket representation $m^{\rho}_{i}$ in hand, we calculate the cosine similarity score:

$$s_{i}=\frac{m^{\mathcal{L}}\cdot m^{\rho}_{i}}{\left\|m^{\mathcal{L}}\right\|_{2}\cdot\left\|m^{\rho}_{i}\right\|_{2}}.$$

(3)

For the candidate pockets $S=\{\varsigma_{1},\varsigma_{2},...\}$, the similarity score $s_{+}$ between the target pocket and the ligand is expected to be higher than others. Thus, we propose the contrastive learning objective:

$$L=-\frac{1}{N}\cdot log\frac{exp(s_{+}/\tau)}{\sum_{i}exp(s_{i}/\tau)},$$

(4)

where $\tau$ is the temperature paramter. For blind docking, the pocket with the highest similarity score with the ligand is selected for the next docking step.

For an iterative refinement module, an initial structure is needed as a starting point. Previous work on molecular 3D conformer generation [17] demonstrates the importance of a good initial structure. Therefore, Bi-EGMN adopts a rapid GPU-accelerated sampling method proposed by Huang et al. [23] to sample a high-quality initial $\mathcal{X}^{\mathcal{L}}$. In this work, the search steps number and the search copy number are set to 40 and 384, respectively. Details about the search box setting can be found in the Appendix.B.3.1.

With input initial structure $\mathcal{X}^{\mathcal{L}}$, we iteratively update it to improve its accuracy. As discussed in Sec.1, the modeling of an entire binding pocket structure is crucial for the success of the process. Current methods either ignore the atom-level structure or model the full-atom pocket structure directly. The latter approach can significantly elevate the computational resource demand, particularly when dealing with large pockets. To overcome these challenges and maintain high docking accuracy and efficiency, we propose a bi-level strategy in this work. In the following sections, we first present the details of the bi-level strategy. Subsequently, we discuss the Bi-EGMN layer, which is used to perform refinement, as depicted in Fig. 1.

Bi-level strategy. The first refinement level is the residue level, where the protein residues within a 40.0 Å  cubic region centered at the geometric centers of ligands are considered as pocket $\rho$. Previous work demonstrates such a range is large enough to cover the binding pocket [30]. In this context, as the full-atom structure of proteins is not considered, the pocket residue graph $\mathcal{G}^{\mathcal{\rho}*}$ is adapted. The second level is the atomic level, where we set the ligand structures refined through $T$ rounds of residue level refinement as the reference structure. In this level, protein atoms within a $6.0$ Å  radius of the ligand atoms are considered to construct pocket atomic graph $\mathcal{G}^{\mathcal{\rho}}$ for modeling the fine-grained interaction. The ligand coordinates $X^{a,\mathcal{L}}$ output by the last layer of atomic level refinement correspond to the final predicted structure $\hat{\mathcal{X}}^{\mathcal{L}}$.

Bi-EGMN Layer. The bi-level E(3)-equivariant graph matching network (Bi-EGMN) layer is the model designed to calculate the protein-ligand interaction and refine the structures. More specifically, this layer adheres to the message-passing paradigm [31] and consists of four functions: intra-message function, inter-message function, aggregate function, and update function.

The intra-message function works to extract messages $m_{i,j}$ and $\hat{m}_{i,j}$ between a node $i$ and its neighbor nodes $j$ from the same molecule graph. $m_{i,j}$ is later used for the updating of node features and $\hat{m}_{i,j}$ for the updating node coordinates. $\forall(i,j)\in\mathcal{E}_{\mathcal{P}}\cup\mathcal{E}_{\mathcal{L}}$, this function can be formally written as :

$\displaystyle d_{i,j}^{(l)}=||x_{i}^{(l)}-x_{j}^{(l)}||,\;m_{i,j}=\varphi_{m}(h_{i}^{(l)},h_{j}^{(l)},d_{i,j}^{(l)},),\;\hat{m}_{i,j}=(x_{i}^{(l)}-x_{j}^{(l)})\cdot\varphi_{\hat{m}}(m_{i,j}),$

(5)

where $d_{i,j}^{(l)}$ is the relative distance between node $i$ and node $j$, and $\varphi$ is a MLP.

The inter-message function works to extract messages $\mu_{i,j}$ and $\hat{\mu}_{i,j}$ between a node $i$ and its neighbor nodes $j$ from the other molecule graphs. Formally, $\forall i\in\mathcal{V}_{\mathcal{P}},j\in\mathcal{V}_{\mathcal{L}}\,or\ i\in\mathcal{V}_{\mathcal{L}},j\in\mathcal{V}_{\mathcal{P}}$:

$\displaystyle\mu_{i,j}=\varphi_{\mu}(h_{i}^{(l)},h_{j}^{(l)},d_{i,j}^{(l)}),\;\hat{\mu}_{i,j}=(x_{i}^{(l)}-x_{j}^{(l)})\cdot\varphi_{\hat{\mu}}(\mu_{i,j}).$

(6)

After extracting inter-message and intra-message, the aggregation function aggregates the neighbor messages of the node $i$. $\forall i\in\mathcal{V}_{\mathcal{P}}\cup\mathcal{V}_{\mathcal{L}}$:

	$\displaystyle m_{i}=\sum_{j\in\mathcal{N}(i)}m_{i,j},\;\hat{m}_{i}=\sum_{j\in\mathcal{N}(i)}\frac{1}{d_{i,j}^{(l)}+1}\cdot\hat{m}_{i,j},$		(7)
	$\displaystyle\mu_{i}=\sum_{j\in\mathcal{N}_{*}^{(l)}(i)}\varphi(\mu_{i,j})\cdot\mu_{i,j},\;\hat{\mu}_{i}=\sum_{j\in\mathcal{N}_{*}^{(l)}(i)}\frac{1}{d_{i,j}^{(l)}+1}\cdot\hat{\mu}_{i,j},$		(8)

where $\mathcal{N}(i)$ is the neighbor of node $i$ in the same graph, and $\mathcal{N}_{*}^{(l)}(i)$ is the set of nodes associated with node $i$ in the other graph.

Finally, the update function updates the position and features of each node:

	$\displaystyle x_{i}^{(l+1)}=\eta x_{i}^{(0)}+(1-\eta)x_{i}^{(l)}+\hat{m}_{i}+\hat{\mu}_{i},\;\forall i\in\mathcal{V}_{\mathcal{P}}\cup\mathcal{V}_{\mathcal{L}},$		(9)
	$\displaystyle h_{i}^{(l+1)}=(1-\beta)\cdot h_{i}^{(l)}+\beta\cdot\varphi(h_{i}^{(l)},m_{i},\mu_{i},h_{i}^{(0)}),\;\forall i\in\mathcal{V}_{\mathcal{P}}\cup\mathcal{V}_{\mathcal{L}},$		(10)

where $\beta$ and $\eta$ are feature skip connection weight and coordinates skip connection weight, respectively. Through such a message-passing paradigm, our Bi-EGMN layers make to update coordinates iteratively.

Lastly, as Bi-EGMN updates structures by modifying the coordinates rather than the torsional angles, as is done in methods like DiffDock [9] and other sampling-based methods, it is crucial to ensure the plausibility of bond lengths and bond angles of the updated structure $\hat{\mathcal{X}}^{\mathcal{L}}$. Therefore, fast structure correction steps, torsion alignment, and SMINA-based energy minimization are designed.

Torsion Alignment. We employ a rapid torsion alignment for the updated structure. The target of this alignment is to align the input structure $\mathcal{X}^{\mathcal{L}}$ with the updated structures $\hat{\mathcal{X}}^{\mathcal{L}}$ by rotating its torsional bonds. Formally, let $(b_{i},c_{i})$ denote a $i$-th rotatable bond, where $b_{i}$ and $c_{i}$ are the starting and ending atoms of the bond, respectively. We randomly select a neighboring atom $a_{i}$ of $b_{i}$ and a neighbor atom $d_{i}$ of $c_{i}$ to calculate the dihedral angle $\hat{\delta}_{i}=\angle(a_{i}b_{i}c_{i},b_{i}c_{i}d_{i})$ based on updated structure coordinates $\hat{\mathcal{X}}^{\mathcal{L}}$. Subsequently, we rotate the rotatable bond $(b_{i},c_{i})$ of input structures to match its dihedral angle $\delta_{i}$ the same as $\hat{\delta}_{i}$. This simple operation can be implemented efficiently using RDKit. After all rotatable bonds have been rotated, we align the rotated input structure to the updated structures to obtain the torsionally aligned structure $\hat{\mathcal{X}}^{\mathcal{L}\mathcal{T}}$. This process ensures the plausibility of bond lengths and bond angles in the torsionally aligned structure $\hat{\mathcal{X}}^{\mathcal{L}\mathcal{T}}$.

Energy Minimization. To further enhance the reliability of DeltaDock, we implement an energy minimization on the torsionally aligned structure $\hat{\mathcal{X}}^{\mathcal{L}\mathcal{T}}$, when an inter-molecular steric clash between the protein and ligand is detected. This energy minimization is conducted using SMINA [19], as it is a highly efficient tool for this process compared with specialized energy minimization tool OpenMM [32] (details see Appendix.A.5). The output structure of this process is $\hat{\mathcal{X}}^{\mathcal{L}{{}^{\prime}}}$.

The training object $L$ is a contrastive object defined before (Eq. 4). For a protein and its candidate pockets set $S=\{\varsigma_{1},\varsigma_{2},...\}$, the positive pair is the target pocket-ligand pair and the negative pairs are other pocket-ligand pairs. The pocket-ligand pairs across different proteins are not used. When training, we calculate the minimum center distance ($DCC_{min}$) between all candidate pockets and the ligand. If $DCC_{min}\leq 5.0$ Å, we add the ligand center into $S$ to assert the existence of positive pairs for every protein (details see Appendix.B.3.1).

We design a physics-informed loss function for the Bi-EGMN module for training. The coordinates $X^{a,\mathcal{L}}$ and $X^{r,\mathcal{L}}$ output by the last layer of atomic level and residue level are both employed in the computation of this loss. Formally, the loss function can be expressed as follows:

$$L=L_{inter}+\lambda_{1}L_{intra}+\lambda_{2}L_{vdw}+\lambda_{3}L_{bound},$$

(11)

where $\lambda$ are weight hyper-parameters. Among the four components, inter-distance map loss $L_{inter}$ is responsible for the RMSD accuracy. Other three items, namely intra-distance map loss $L_{intra}$, vdw constraint loss $L_{vdw}$, and bound matrix constraint loss $L_{bound}$ are employed for physical validity. When training and inferencing, we follow previous work [33] and employ the recycling strategy (details see Appendix.B.3.2).

As demonstrated in Table.1, DeltaDock outperforms all baseline methods. Specifically, DeltaDock achieves a remarkable success rate of 47.4% (where RMSD < 2.0 Å), surpassing the previous SOTA GDL method, DiffDock, which has a success rate of 36.0%. Recent GPU-accelerated docking methods have also made significant progress in blind docking. However, when compared to DSDP, which is the top-performing sampling-based method in the PDBbind test set, DeltaDock still exhibits superior performance across all metrics. Notably, as elucidated in Section 3.3, DeltaDock employs the same sampling algorithm as DSDP for generating the initial structure. Yet, our framework allows DeltaDock to significantly outperform DSDP.

Beyond accuracy, efficiency is a critical performance measure for molecular docking methods. As indicated in Table 1, DeltaDock maintains a competitive level of efficiency, despite the inclusion of an energy minimization operation to enhance accuracy and reliability. Molecular docking methods invariably face a trade-off between efficiency and accuracy. However, the data presented in Table 1 suggest that DeltaDock could serve as a viable tool for practical applications, balancing these two crucial aspects effectively.

Most existing GDL methods, such as DiffDock and EquiBind, are primarily designed for blind docking scenarios and are not inherently suited for site-specific docking tasks. However, DeltaDock seamlessly integrates blind docking and site-specific docking settings. In this context, the pocket is directly provided, eliminating the need for pocket selection via CPLA. The performance of DeltaDock in site-specific docking is illustrated in Fig.2. When supplied with predefined binding sites, traditional sampling methods exhibit a significant improvement in results. For instance, the docking success rate of VINA escalates from 10.3% to 45.0%. Despite this enhancement, DeltaDock consistently surpasses all baselines. Previous research suggested that while GDL docking methods excel at pocket searching, traditional methods tend to outperform GDL models in site-specific docking tasks [35]. However, as evidenced by the results presented in Table.1 and Fig.2, DeltaDock exhibits superior performance in both blind and site-specific docking scenarios, demonstrating its versatility and robustness in handling diverse docking settings.

Beyond the overall docking performance, the pocket-ligand alignment and iterative refinement results are explored (Fig. 3). As depicted in the figure, CPLA predicts significantly more accurate pockets than other methods and Bi-EGMN can diminish the discrepancy between ground-truth structures and input structures. Generally, the PDBbind test set poses a more significant challenge to Bi-EGMN than the PoseBusters dataset. And for CPLA, PoseBusters dataset is more challenging otherwise. The consistent good performance on the two datasets demonstrates the effectiveness and generalization capacity of CPLA and Bi-EGMN.

In this section, ablation studies are conducted to assess the contributions of different components. We first ablate the whole CPLA or Bi-EGMN, and then the residue-level or the atom-level in Bi-EGMN (see Appendix. B.4 for implement details). As illustrated in Table 2, it becomes clear that each component, encompassing CPLA and the bi-level strategy in Bi-EGMN, plays a significant role in enhancing the overall performance of DeltaDock. Due to the space limitation, a full ablation study can be found in Appendix. C.3.