Reprogramming Pretrained Target-Specific Diffusion Models for Dual-Target Drug Design

To collect drug combination pairs, we start from DrugCombDB¹¹1http://drugcombdb.denglab.org/main [35]. DrugCombDB is a comprehensive database devoted to the curation of drug combinations from various data sources including high-throughput screening (HTS) assays, manual curations from the literature, FDA Orange Book and external databases. DrugCombDB comprises a total of 448,555 combinations of drugs, encompassing 2,887 unique drugs and 124 human cancer cell lines. Particularly, DrugCombDB has more than 6,000,000 quantitative dose responses, from which we determine whether a drug combination is synergistic or not. Specifically, a drug combination with positive zero interaction potency (ZIP), Bliss, Loewe and the highest single agent (HSA) scores simultaneously in at least one cell line is supposed to be a synergistic one. Please refer to Appendix C for a comprehensive understanding of these scores.

After collecting synergistic drug combinations, we need to collect other necessary information (e.g., SMILES and targets) according to their drug names provided by DrugCombDB. Before this procedure, we collect synonyms and cross-matching ID (e.g., CAS Number and ChEBI ID) mainly from DrugBank²²2https://go.drugbank.com/ [30] and Therapeutic Target Database (TTD)³³3https://idrblab.net/ttd/ [67]. This step facilitates comprehensive literature reviews, ensuring that all relevant data sources that may use alternate names for a drug is considered. We then collect SMILES or structures (if possible) also mainly from DrugBank and TTD. To identify drug targets, we also utilize DrugBank and TTD as the primary data sources, and supplement these with manual curation from the literature (e.g., [28]). For drugs for which we cannot find either SMILES or targets, we exclude them from our previously collected dataset of positive drug combinations.

For certain drug-target pairs, we incorporate their complex structures directly into our dataset if they are available in PDBBind [38], a repository of protein-ligand binding structures sourced from the Protein Data Bank (PDB) [2]. For drug-target pairs not present in PDBBind, we initially attempt to retrieve the target structures from PDB; if the structures are unavailable, we then source them from the AlphaFold Protein Structure Database (AlphaFold DB) [59] and exclude those whose confidence scores, referred to as pLDDT which provide an assessment of the structures predicted by AlphaFold 2 [29], are less than 70. For these protein targets with structures from PDB or AlphaFold DB, we first utilize P2Rank [31], a program that precisely predicts ligand-binding pockets from a protein structure, to find the most possible pocket given the target structure, and use AutoDock Vina [12] to obtain the protein-ligand complex structures. For each drug, there may exist more than one targets, in which case we use AutoDock Vina to measure the binding affinity and selected the target with the best binding affinity. Finally, we obtain 12,917 postive drug combinations with protein-ligand complex structures, among which there are 438 unique drugs. The 12,917 pairs of targets can be used to evaluate the ability of methods for dual-target drug design. And the binding ligands can be used for reference molecules.

In this following, we denote the type of an atom as ${\bm{v}}\in\mathbb{R}^{K}$ and the coordinate of an atom as ${\bm{x}}\in\mathbb{R}^{3}$, where $K$ is the number of atom types of our interest. For single-target drug design, given a protein binding site denoted as a set of atoms ${\mathcal{P}}=\{({\bm{x}}_{P}^{(i)},{\bm{v}}_{P}^{(i)})\}_{i=1}^{N_{P}}$, where $N_{P}$ is the number of protein atoms, our goal is to generate binding molecules ${\mathcal{M}}=\{({\bm{x}}^{(i)}_{L},{\bm{v}}_{L}^{(i)})\}_{i=1}^{N_{M}}$. For brevity, we denote the molecule as $M=[{\mathbf{x}},{\mathbf{v}}]$, where $[\cdot,\cdot]$ denotes the concatenation operator and ${\mathbf{x}}\in\mathbb{R}^{N_{M}\times 3},{\mathbf{x}}\in\mathbb{R}^{N_{M}\times K}$ denote the coordinates in 3D space and one-hot atom types, respectively. So we can use generative models to model the conditional distribution $p(M|{\mathcal{P}})$.

In the forward diffusion process of the diffusion model, noises are gradually injected into the data sample (i.e., small molecule $M_{0}\sim p(M|{\mathcal{P}})$) and lead to a sequence of latent variable $M_{1},M_{2},\ldots,M_{T}$. The final distribution $p(M_{T}|{\mathcal{P}})$, also known as prior distribution, is approximately standard normal distribution for atom positions and uniform distribution for atom types. The reverse generative process learns to recover data distribution from the noise distribution with a neural network parameterized by ${\bm{\theta}}$. The forward and reverse processes are both Markov chains defined as follows:

$\displaystyle q(M_{1:T}|M_{0},{\mathcal{P}})=\prod_{t=1}^{T}q(M_{t}|M_{t-1},{\mathcal{P}})\quad\text{and}\quad p_{\bm{\theta}}(M_{0:T-1}|M_{T},{\mathcal{P}})=\prod_{t=1}^{T}p_{\bm{\theta}}(M_{t-1}|M_{t},{\mathcal{P}}).$

(1)

More specifically, the forward transition kernel in Guan et al. [19] are defined as follows:

$\displaystyle q(M_{t}|M_{t-1},{\mathcal{P}})={\mathcal{N}}({\mathbf{x}}_{t};\sqrt{1-\beta_{t}}{\mathbf{x}}_{t-1},\beta_{t}{\bm{I}})\cdot{\mathcal{C}}({\mathbf{v}}_{t}|(1-\beta_{t}){\mathbf{v}}_{t-1}+\beta_{t}/K),$

(2)

where $\{\beta_{t}\}_{t=1}^{T}$ are fixed noise schedule. The above diffusion process can be efficiently sampled directly from time step $0$ to $t$ as follows:

$\displaystyle q({\mathbf{x}}_{t}|{\mathbf{x}}_{0})={\mathcal{N}}({\mathbf{x}}_{t};\sqrt{\bar{\alpha}_{t}}{\mathbf{x}}_{0},(1-\bar{\alpha}_{t}){\bm{I}})\quad\text{and}\quad q({\mathbf{v}}_{t}|{\mathbf{v}}_{0})={\mathcal{C}}({\mathbf{v}}_{t}|\bar{\alpha}_{t}{\mathbf{v}}_{0}+(1-\bar{\alpha}_{t}/K),$

(3)

where $\alpha_{t}:=1-\beta_{t}$ and $\bar{\alpha}_{t}:=\prod_{s=1}^{t}\alpha_{s}$. The posterior can be easily computed via Bayes theorem as:

$\displaystyle q({\mathbf{x}}_{t-1}|{\mathbf{x}}_{t},{\mathbf{x}}_{0})={\mathcal{N}}({\mathbf{x}}_{t-1};\tilde{{\bm{\mu}}}_{t}({\mathbf{x}}_{t},{\mathbf{x}}_{0}),\tilde{\beta}_{t}{\bm{I}})\quad\text{and}\quad q({\mathbf{v}}_{t}|{\mathbf{v}}_{0})={\mathcal{C}}({\mathbf{v}}_{t}|\bar{\alpha}_{t}{\mathbf{v}}_{0}+(1-\bar{\alpha}_{t})/K),$

(4)

where $\tilde{\beta}_{t}=\frac{1-\bar{\alpha}_{t-1}}{1-\bar{\alpha}_{t}}\beta_{t}$, $\tilde{\mu}_{t}({\mathbf{x}}_{t},{\mathbf{x}}_{0})=\frac{\sqrt{\bar{\alpha}_{t-1}}\beta_{t}}{1-\bar{\alpha}_{t}}$, $\tilde{{\bm{c}}}_{t}({\mathbf{v}}_{t},{\mathbf{v}}_{0})={\bm{c}}^{*}/\sum_{k=1}^{K}c^{*}_{k}$ and ${\bm{c}}^{*}({\mathbf{v}}_{t},{\mathbf{v}}_{0})=[\alpha_{t}{\mathbf{v}}_{t}+(1-\alpha_{t})/K]\odot[\bar{\alpha}_{t-1}{\mathbf{v}}_{0}+(1-\bar{\alpha}_{t-1})/K]$.

Accordingly, the reverse transition kernel are defined as follows:

$\displaystyle p_{\bm{\theta}}(M_{t-1}|M_{t},{\mathcal{P}})={\mathcal{N}}({\mathbf{x}}_{t-1};{\bm{\mu}}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}}),\sigma_{t}^{2}{\bm{I}})\cdot{\mathcal{C}}({\mathbf{v}}_{t-1}|{\mathbf{c}}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})).$

(5)

Guan et al. [19] use SE(3)-equivariant neural networks [52, 18] to parameterize ${\bm{\mu}}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})$ and ${\mathbf{c}}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})$. More specifically, the $[{\mathbf{x}}_{0},{\mathbf{v}}_{0}]$ are first predicted using neural network $f_{\bm{\theta}}$, i.e., $[\hat{{\mathbf{x}}}_{0},\hat{{\mathbf{v}}}_{0}]=f_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})$ and then substitute in the posterior as in Equation 4. At the $l$-th layer of $f_{\bm{\theta}}$, the hidden embedding ${\mathbf{h}}$ and coordinates ${\mathbf{x}}$ of each atom are updated alternately as follows:

	$\displaystyle{\mathbf{h}}^{l+1}_{i}={\mathbf{h}}^{l}_{i}+\sum_{j\in{\mathcal{V}},i\neq j}f_{{\bm{\theta}}_{h}}({\mathbf{h}}^{l}_{i},{\mathbf{h}}^{l}_{j},d^{l}_{ij},{\mathbf{e}}_{ij}),$		(6)
	$\displaystyle{\mathbf{x}}^{l+1}_{i}={\mathbf{x}}^{l}_{i}+\sum_{j\in{\mathcal{V}},i\neq j}({\mathbf{x}}^{l}_{i}-{\mathbf{x}}^{l}_{j})f_{{\bm{\theta}}_{x}}({\mathbf{h}}^{l+1}_{i},{\mathbf{h}}^{l+1}_{j},d^{l}_{ij},{\mathbf{e}}_{ij})\cdot\mathbf{1}_{\text{ligand}},$		(7)

where ${\mathcal{V}}$ is a k-nearest neighbors (knn) graph, $d_{ij}=\|{\mathbf{x}}_{i}-{\mathbf{x}}_{j}\|$ is the Euclidean distance between two atoms $i$ and $j$, ${\mathbf{e}}_{ij}$ is an additional feature that indicates the connection is between protein atoms, ligand atoms or protein atom and ligand atom, and $\mathbf{1}_{\text{ligand}}$ is a mask for ligand nodes since only coordinates of ligand atoms are supposed to be updated.

The diffusion model is trained to minimize the KL-divergence between the ground-truth posterior $q(M_{t-1}|M_{0},M_{t},{\mathcal{P}})$ and the estimated posterior $p_{\bm{\theta}}(M_{t-1}|M_{t},{\mathcal{P}})$. After being trained, given a specific pocket, the ligand molecule can be generated by first sampling from prior distribution and sequentially applying the reverse generative process defined above.

The goal of dual-target drug design is to design a ligand molecule $M$ that can bind to both given pocket ${\mathcal{P}}_{1}$ and pocket ${\mathcal{P}}_{2}$. The problem can be also formulated as a generative task which models the conditional distribution $p(M|{\mathcal{P}}_{1},{\mathcal{T}}{\mathcal{P}}_{2})$. Notebly, we introduce a transformation operator ${\mathcal{T}}$ here. This is because protein pockets exhibit a wide variety of shapes and chemical characteristics and it is necessary to achieve spatial alignment of the dual pockets when modeling the conditional distribution with both of them as conditions. To maintain a neat but siginificant setting, we restrict the transformation ${\mathcal{T}}$ to encompass solely translations ${\mathcal{T}}_{T}$ and rotations ${\mathcal{T}}_{R}$, i.e., ${\mathcal{T}}={\mathcal{T}}_{T}\circ{\mathcal{T}}_{R}$. To be more precise, $({\mathcal{T}}_{T}\circ{\mathcal{T}}_{R}){\mathcal{P}}_{2}=\{({\bm{R}}{\bm{x}}_{P_{2}}^{(i)}+{\bm{t}}),{\bm{v}}^{(i)}_{P_{2}}\}_{i=1}^{N_{P_{2}}}$, where ${\bm{R}}\in\text{SO(3)}$ represents the rotation and ${\bm{t}}\in\mathbb{R}^{3}$ represents the translation.

Protein pockets can have intricate 3D structures with varying depths, widths, and surface contours, and distinct chemical properties, including differences in surface electron potentials, hydrophobicity, and the distribution of functional groups. The complex nature of pocket characteristics hinder us from directly aligning two pockets. Nevertheless, the binding mode is determined by the complex nature of pocket information, thus the protein-binding priors can effectively summarize the essential information needed for aligning the two pockets. This motivate us to propose to use a ligand molecule as a prober to implicitly reflect the spatial arrangement of the two pockets and then align them with the protein-ligand binding priors. More specifically, we first dock a ligand molecule to pocket ${\mathcal{P}}1$ and pocket ${\mathcal{P}}2$ separately. We can then compute ${\bm{R}}$ and ${\bm{t}}$ by aligning the two docked poses of the ligand molecule. Experiments have demonstrated that even ligand molecules capable of approximate binding to the two pockets can effectively indicate a specific spatial alignment between them. Further details are provided in Section 4.

Inspired by compositional visual generation [9, 37], we can model the dual-target drug design with the following composed distribution:

$\displaystyle p(M|{\mathcal{P}}_{1},{\mathcal{T}}{\mathcal{P}}_{2})\propto p_{\bm{\theta}}(M|{\mathcal{P}}_{1})p_{\bm{\theta}}(M|{\mathcal{T}}{\mathcal{P}}_{2}).$

(8)

Following Liu et al. [37], for the atom position prediction, we can reparameterize ${\bm{\mu}}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})$ with ${\mathbf{x}}_{t}-\bm{\epsilon}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})$, so that we can rewrite the transition kernel of the reverse generative process as follows:

$\displaystyle p_{\bm{\theta}}({\bm{x}}_{t-1}|{\bm{x}}_{t},{\mathcal{P}})={\mathcal{N}}({\mathbf{x}}_{t-1};{\bm{\mu}}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}}),\sigma_{t}^{2}{\bm{I}})={\mathcal{N}}({\mathbf{x}}_{t-1};{\mathbf{x}}_{t}-\bm{\epsilon}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}}),\sigma_{t}^{2}{\bm{I}}).$

(9)

The reversed transition kernel corresponds to a step as follow:

$\displaystyle{\mathbf{x}}_{t-1}={\mathbf{x}}_{t}-\bm{\epsilon}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})+{\mathcal{N}}({\bm{0}},\sigma_{t}^{2}{\bm{I}}).$

(10)

where $\bm{\epsilon}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})$ can be viewed as a drift term and ${\mathcal{N}}({\bm{0}},\sigma_{t}^{2}{\bm{I}})$ can be viewed as a diffusion term. As Liu et al. [37] points out, this is analogous to the Langevin dynamics [11] that used to sample from Energy-Based Models (EBMs) [10, 55], which can be formulated as follows:

$\displaystyle{\mathbf{x}}_{t-1}={\mathbf{x}}_{t}-\frac{\lambda}{2}\nabla_{{\mathbf{x}}}E_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}})+{\mathcal{N}}({\bm{0}},\sigma_{t}^{2}{\bm{I}}).$

(11)

The sampling procedure produces samples from the probability density $p_{\bm{\theta}}({\mathbf{x}}|{\mathcal{P}})\propto\exp{(-E_{\bm{\theta}}({\mathbf{x}}|{\mathcal{P}}))}$ where $E_{\bm{\theta}}({\mathbf{x}}|{\mathcal{P}})$ is a energy function parameterized by model ${\bm{\theta}}$. Thus, accordingly, the composed distribution of atom positions can be written as:

$\displaystyle p({\mathbf{x}}|{\mathcal{P}}_{1},{\mathcal{T}}{\mathcal{P}}_{2})\propto p_{\bm{\theta}}({\mathbf{x}}|{\mathcal{P}}_{1})p_{\bm{\theta}}({\mathbf{x}}|{\mathcal{T}}{\mathcal{P}}_{2})\propto\exp\Big{(}-\big{(}E_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}}_{1})+E_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{T}}{\mathcal{P}}_{2})\big{)}\Big{)}.$

This composed distribution corresponds to Langevin dynanmics as follows:

$\displaystyle{\mathbf{x}}_{t-1}={\mathbf{x}}_{t}-\frac{\lambda}{2}\nabla_{\mathbf{x}}\Big{(}-\big{(}E_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}}_{1})+E_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{T}}{\mathcal{P}}_{2})\big{)}\Big{)}.$

(12)

Accordingly, each step in the compositional reverse generative sampling process can be defined as:

$\displaystyle{\mathbf{x}}_{t-1}={\mathbf{x}}_{t}-\eta\Big{(}\big{(}\bm{\epsilon}_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}}_{1})+\epsilon_{\bm{\theta}}([{\mathbf{x}}_{t},{\mathbf{v}}_{t}],t,{\mathcal{P}}_{2})\Big{)},$

(13)

where we additionally introduce a hyperparameter $\eta$ here to control the strength of the drift. In theory, this is equivalent to making a more flexible assumption that $p({\mathbf{x}}|{\mathcal{P}}_{1},{\mathcal{T}}{\mathcal{P}}_{2})\propto[p_{\bm{\theta}}({\mathbf{x}}|{\mathcal{P}}_{1})p_{\bm{\theta}}({\mathbf{x}}|{\mathcal{T}}{\mathcal{P}}_{2})]^{\eta}$. In this case the transition kernel is defined as $p_{\bm{\theta}}({\mathbf{x}}_{t-1}|{\mathbf{x}}_{t},{\mathcal{P}}_{1},{\mathcal{T}}{\mathcal{P}}_{2})=[p_{\bm{\theta}}({\mathbf{x}}_{t-1}|{\mathbf{x}}_{t},{\mathcal{P}}_{1})p_{\bm{\theta}}({\mathbf{x}}_{t-1}|{\mathbf{x}}_{t},{\mathcal{T}}{\mathcal{P}}_{2})]^{\eta}$. In practice, we set $\eta=1/2$ by default. This composition operation is equivalent to averaging two $\hat{{\mathbf{x}}}_{0}$ predicted on two complex graphs, i.e. ${\mathcal{V}}_{1}$ and ${\mathcal{V}}_{2}$. Similarly, for the atom types which are discrete variables, we can also compose the transition kernel as follows: $p_{\bm{\theta}}({\mathbf{v}}_{t-1}|{\mathbf{v}}_{t},{\mathcal{P}}_{1},{\mathcal{T}}{\mathcal{P}}_{2})\propto p_{\bm{\theta}}({\mathbf{v}}_{t-1}|{\mathbf{v}}_{t},{\mathcal{P}}_{1})p_{\bm{\theta}}({\mathbf{v}}_{t-1}|{\mathbf{v}}_{t},{\mathcal{T}}{\mathcal{P}}_{2})$. Note that $p_{\bm{\theta}}({\mathbf{v}}_{t-1}|{\mathbf{v}}_{t},{\mathcal{P}}_{1})$ is categorical distribution and its dimension (i.e., the number of atom types of our interest) is $K$, which is small in practice. So $p_{\bm{\theta}}({\mathbf{v}}_{t-1}|{\mathbf{v}}_{t},{\mathcal{P}}_{1},{\mathcal{T}}{\mathcal{P}}_{2})$ can be computed analytically. We name the compositional reverse sampling with composed transition kernel (i.e., composed drift) as CompDiff.

We further improve the compositional reverse sampling by introducing the composition into each layer of the equivariant neural network in the pretrained diffusion model. For brevity, we denote the SE(3)-equivariant message at the $l$-th layer for the $i$-th atom of the complex graph ${\mathcal{V}}_{v}$ ($v=1,2$) as introduced in Equation 7 as follows:

	$\displaystyle\Delta{\mathbf{h}}^{l}_{i}({\mathcal{V}}_{v}):=\sum_{j\in{\mathcal{V}}_{v},i\neq j}f_{{\bm{\theta}}_{h}}({\mathbf{h}}^{l}_{i},{\mathbf{h}}^{l}_{j},d^{l}_{ij},{\mathbf{e}}_{ij}),$		(14)
	$\displaystyle\Delta{\mathbf{x}}^{l}_{i}({\mathcal{V}}_{v}):=\sum_{j\in{\mathcal{V}}_{v},i\neq j}({\mathbf{x}}^{l}_{i}-{\mathbf{x}}^{l}_{j})f_{{\bm{\theta}}_{x}}({\mathbf{h}}^{l+1}_{i},{\mathbf{h}}^{l+1}_{j},d^{l}_{ij},{\mathbf{e}}_{ij})\cdot\mathbf{1}_{\text{ligand}}.$		(15)

The above SE(3)-equivariant message can also be interpreted as drift in 3D and latent space. Thus we can also compose them as follows:

$\displaystyle{\mathbf{h}}^{l+1}_{i}={\mathbf{h}}^{l}_{i}+\big{(}\Delta{\mathbf{h}}^{l}_{i}({\mathcal{V}}_{1})+\Delta{\mathbf{h}}^{l}_{i}({\mathcal{V}}_{2})\big{)}/2\quad\text{and}\quad{\mathbf{x}}^{l+1}_{i}={\mathbf{x}}^{l}_{i}+(\Delta{\mathbf{x}}^{l}_{i}({\mathcal{V}}_{1})+\Delta{\mathbf{x}}^{l}_{i}({\mathcal{V}}_{2}))/2.$

(16)

We name the compositional reverse sampling with the above SE(3)-equivariant message at each layer as DualDiff. The proof of SE(3)-equivariance can be found in Appendix G. This more meticulous composition is supposed to lead to higher-quality samples.

We use our dataset introduced in Section 3.1. All 12,917 pairs of targets (including 438 unique targets) are used for evaluation. For each target, there is an associated reference molecule that can be considered a benchmark for high-quality ligand molecules and utilized in linker design methods.

We compare our method with various baselines: Pocket2Mol [45] generates 3D molecules atom by atom in an autoregressive manner given a specific protein binding site. TargetDiff [19] is a diffusion-based method which generates atom coordinates and atom types in a non-autoregressive way. Note that Pocket2Mol and TargetDiff are both proposed for structure-based single-target drug design. DiffLinker [25] is a diffusion-based model for linker design with given fragment poses. LinkerNet [17] is a diffusion-based model for co-designing molecular fragment poses and the linker. DiffLinker and LinkerNet are repurposed for dual-target drug design as we introduced in Section 3.3. The code is available at https://github.com/zhouxiangxin1998/DualDiff.

We evaluate generated ligand molecules from the perspectives of target binding affinity and molecular properties. We employ AutoDock Vina [12] to estimate the target binding affinity, following Peng et al. [45], Guan et al. [19]. We first evaluate Pocket2Mol and TargetDiff under the single-target setting as a preliminary verification (see Appendix E). For dual-target drug design, we use each method to design 10 molecules for each pair of targets, denoted as ${\mathcal{P}}_{1}$ and ${\mathcal{P}}_{2}$. (For reference molecule, Pocet2Mol and TargetDiff, the ligand molecule is generated for ${\mathcal{P}}_{1}$ but the target binding affinity is evaluated on both ${\mathcal{P}}_{1}$ and ${\mathcal{P}}_{2}$.) We then collect all generated molecules across 12,917 pairs of targets and report the mean and median (denoted as “Avg.” and “Med.” respectively) of affinity-related metrics (P-1 Vina Dock, P-2 Vina Dock, Max Vina Dock, and Dual High Affinity) and property-related metrics (drug-likeness QED [3], synthesizability SA [14], and diversity). Vina Dock incorporates a re-docking step to assess the highest binding affinity achievable. Here we introduce P-1 Vina Dock and P-2 Vina Dock to represent the Vina Dock score evaluated on ${\mathcal{P}}_{1}$ and ${\mathcal{P}}_{2}$, respectively. Besides, we introduce Max Vina Dock, which represents the maximum Vina Dock of a given molecule towards ${\mathcal{P}}_{1}$ and ${\mathcal{P}}_{2}$. The Vina Dock will be low if and only if the molecule can bind to both targets simultaneously, which is the goal of dual-target drug design. Additionally, we report Dual High Affinity (abbreviated as Dual High Aff.) which represents the proportion of generated molecules that exhibit binding affinity that exceeds that of the reference molecules on both the respective targets. This reflects the success rate in achieving higher binding affinities simultaneously on both targets in dual-target drug design. We also evaluate the RMSD between docked poses towards dual targets.

We perform different methods to align the pockets of dual targets for CompDiff and DualDiff. See the results in Table 2. Naively, we can align two pockets by their geometric centers (denoted as “-Center”). In our method, we propose to align pockets with protein-ligand binding priors. We select the ligand with minimum RMSD between docked poses (resp. minimum sum of Vina Dock scores) towards dual targets as the anchor to align the dual targets, which is denoted as “-RMSD” (resp. “-Score”). DualDiff-Score achieved the best performance among all variants, demonstrating the effectiveness of the alignment method.

We conduct ablation on different strategies of identifying fragments for DiffLinker and LinkerNet for dual-target drug design. Since we use Vina Dock to select fragments, we have tried different box sizes for docking, i.e., 5Å and 8Å. For DiffLinker, since the relative poses of fragments are required as input, we dock fragments derived from ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ to target ${\mathcal{P}}_{1}$ (or ${\mathcal{P}}_{2}$). We then apply DiffLinker both with and without considering the pocket. The corresponding methods are denoted as “DiffLinker-5/-8/-pocket-5/-pocket-8”. For LinkerNet, the relative poses of fragments are not required. So we try two setting: dock fragments from ${\mathcal{M}}_{1}$ and ${\mathcal{M}}_{2}$ to target ${\mathcal{P}}_{1}$ (or ${\mathcal{P}}_{2}$); dock fragments from ${\mathcal{M}}_{1}$ (resp. ${\mathcal{M}}_{2}$) to target ${\mathcal{P}}_{1}$ (resp. ${\mathcal{P}}_{2}$). The difference is that all fragments are docked to the same pocket in the former setting while the fragments from two ligand molecules are docked to their respective pockets. The corresponding methods are denoted as “LinkerNet-5/-8/-self-5/-self-8”. The results are reported in Table 3. The results show that a larger docking box size allows for better selection of fragments. As we expected, DiffLinker with consideration of pockets achieves better performance. Among all variants, “LinkerNet-self-8” achieves best performance, which implies that adjusting relative poses of fragments may play a crutial role in linker deisgn. This feature of LinkerNet allows for selecting the most important fragments for each pockets of the dual targets.