ShEPhERD: Diffusing shape, electrostatics, and pharmacophores for bioisosteric drug design

## Abstract

Abstract Engineering molecules to exhibit precise 3D intermolecular interactions with their environment forms the basis of chemical design.
In ligand-based drug design, bioisosteric analogues of known bioactive hits are often identified by virtually screening chemical libraries with shape, electrostatic, and pharmacophore similarity scoring functions.
We instead hypothesize that a generative model which learns the joint distribution over 3D molecular structures and their interaction profiles may facilitate 3D interaction-aware chemical design.
We specifically design ShEPhERD 1 1 1 ShEPhERD : Sh ape, E lectrostatics, and Ph armacophore E xplicit R epresentation D iffusion , an SE(3)-equivariant diffusion model which jointly diffuses/denoises 3D molecular graphs and representations of their shapes, electrostatic potential surfaces, and (directional) pharmacophores to/from Gaussian noise. Inspired by traditional ligand discovery, we compose 3D similarity scoring functions to assess ShEPhERD ’s ability to conditionally generate novel molecules with desired interaction profiles. We demonstrate ShEPhERD ’s potential for impact via exemplary drug design tasks including natural product ligand hopping, protein-blind bioactive hit diversification, and bioisosteric fragment merging.

## 1 Introduction

Designing new molecules to attain specific functions via physicochemical interactions with their environment is a foundational task across the chemical sciences.
For instance, early-stage drug discovery often involves tuning the 3D shape, electrostatic potential (ESP) surface, and non-covalent interactions of small-molecule ligands to promote selective binding to a protein target (Bissantz et al., ; Huggins et al., ).
Homogeneous catalyst design requires developing organometallic, organic, or even peptidic catalysts that stabilize reactive transition states via specific noncovalent interactions (Raynal et al., ; Fanourakis et al., ; Knowles & Jacobsen, ; Wheeler et al., ; Toste et al., ; Metrano & Miller, ).
Supramolecular chemistry similarly optimizes host-guest interactions for applications across photoresponsive materials design, biomedicine, and structure-directed zeolite synthesis (Qu et al., ; Stoffelen & Huskens, ; Corma et al., ).

This essential challenge of designing molecules with targeted 3D interactions manifests across myriad tasks in ligand-based drug design. In medicinal chemistry, bioisosteric scaffold hopping aims to swap-out substructures within a larger molecule while preserving bioactivity (Langdon et al., ). Often, the swapped scaffolds share biochemically-relevant features such as electrostatics or pharmacophores, which describe both non-directional (hydrophobic, ionic) and directional (hydrogen bond acceptor/donor, aromatic $\pi$-$\pi$, halogen bonding) non-covalent interactions. When scaffold hopping extends to entire ligands, “ligand hopping” can help identify synthesizable analogues of complex natural products that mimic their 3D biochemical interactions (Grisoni et al., ). In hit expansion, ligand hopping is also used to diversify known bioactive hits by proposing alternative actives, ranging from topologically-similar “me-too” compounds to distinctly new chemotypes (Schneider et al., ; Wermuth, ). Notably, ligand hopping does not require knowledge of the protein target. Lastly, bioisosteric scaffold hopping extends beyond altering individual molecules;
Wills et al. () used “bioisosteric fragment merging” to replace a set of fragments that independently bind a protein with one ligand that captures the fragments’ aggregate 3D binding interactions.

Here, we consider three archetypal tasks in bioisosteric drug design: 3D similarity-constrained ligand hopping, protein-blind bioactive hit diversification, and bioisosteric fragment merging (Fig. 1). The unifying theme across these tasks is identifying new molecular structures which are chemically dissimilar from known matter in terms of their molecular graphs, but which are highly similar with respect to their 3D intermolecular interaction profiles.
To find such bioisosteric analogues, traditional design campaigns will virtually screen chemical libraries with 3D similarity scoring functions to query molecules’ 3D shape, electrostatic, and/or pharmacophore similarity with respect to a reference molecule (Oprea & Matter, ; Rush et al., ; Zavodszky et al., ; Sanders et al., ). However, similarity-based virtual screening has acute drawbacks: it cannot explore beyond known chemical libraries, it is inefficient by virtue of being undirected, and it can be quite slow when multiple geometries of conformationally-flexible molecules must be scored.

Figure: Figure 1: We introduce ShEPhERD, a diffusion model that jointly generates 3D molecules and their shapes, electrostatics, and pharmacophores. By explicitly modeling 3D molecular interactions, ShEPhERD can be applied across myriad challenging ligand-based drug design tasks including natural product ligand hopping, bioactive hit diversification, and bioisosteric fragment merging.
Refer to caption: extracted/5947291/figures/shepherd_tasks_v11.jpg

We instead develop a broadly applicable generative modeling framework that enables efficient 3D bioisosteric molecular design. We are motivated by multiple observations:
(1) As elaborated above, ligand-based drug discovery requires designing novel molecular structures that form specific 3D interactions.
(2) Harris et al. () have found that 3D generative models for structure-based drug design struggle to generate ligands that form biochemical interactions (e.g., hydrogen bonds) with the protein, despite training on protein-ligand complexes. This suggests the need for new strategies which explicitly model molecular interactions.
(3) Numerous chemical design scenarios beyond drug design require engineering the physicochemical interactions of small molecules. But, such settings are often data-restricted, necessitating zero/few-shot generative approaches. To address these challenges, we introduce ShEPhERD, a 3D generative model which learns the relationship between the 3D chemical structures of small molecules and their shapes, electrostatics, and pharmacophores (henceforth collectively called “interaction profiles”) in context-free environments. Specifically:

- •
We define explicit point cloud-based representations of molecular shapes, ESP surfaces, and pharmacophores that are amenable to symmetry-preserving SE(3)-equivariant diffusion modeling.
- •
We formulate a joint denoising diffusion probabilistic model (DDPM) that learns the joint distribution over 3D molecules (atom types, bond types, coordinates) and their 3D shapes, ESP surfaces, and pharmacophores. In addition to diffusing/denoising attributed point clouds for shape and electrostatics, we natively model the directionality of pharmacophores by diffusing/denoising vectors on the unit sphere. We sample from specific interaction-conditioned distributions via inpainting.
- •
Inspired by virtual screening, we craft shape, electrostatic, and pharmacophore similarity functions to score (1) the self-consistency of jointly generated molecules and interaction profiles, and (2) the similarity between conditionally generated molecules and target interaction profiles. We show that ShEPhERD can generate diverse 3D molecular structures with substantially enriched interaction-similarity to target profiles, even upon geometric relaxation with semi-empirical DFT.
- •
After training on drug-like datasets, we demonstrate in silico that ShEPhERD can facilely design small-molecule mimics of natural products via ligand hopping, diversify bioactive hits while preserving protein-binding modes, and merge fragments into bioisosteric ligands, all out-of-the-box.

We anticipate that ShEPhERD will prove immediately useful for ligand-based drug design campaigns that require the de novo design of new molecules with precise 3D interactions. However, we stress ShEPhERD’s general applicability to other areas of interaction-aware molecular design, such as organocatalyst design. We especially envision that ShEPhERD will be extended to model other structural characteristics beyond the shape, electrostatic, and pharmacophore profiles treated here.

## 2 Related Work

3D similarity scoring for ligand-based drug design. Ligand-based drug design commonly applies shape, electrostatic, and/or pharmacophore similarity scoring functions to screen for molecules with similar 3D interactions as a known bioactive molecule (Fiedler et al., ; Rush et al., ; Jackson et al., ).
Shape similarity functions typically use atom-centered Gaussians to compute the volumetric overlap between active and query molecules (Grant & Pickup, ; Grant et al., ). Many methods additionally attribute scalar (e.g., charge) or categorical (e.g., pharmacophore type) features to these Gaussians in order to score electrostatic (Good et al., ; Bolcato et al., ; OpenEye, Cadence Molecular Sciences, ) or pharmacophore similarity (Sprague, ; Dixon et al., ; Taminau et al., ; Wahl, ). To better capture intermolecular interactions,
other methods quantify the Coulombic or pharmacophoric “potential” felt by a chemical probe at points near the molecule’s surface (Cheeseright et al., ; Vainio et al., ; Cleves et al., ). We similarly use surface representations of molecular shape and electrostatic interactions, but combine this with volumetric point-cloud and vector representations of pharmacophores to model directional interactions
such as hydrogen bonding and non-covalent $\pi$-effects (Wahl, ). Notably, we develop our own 3D similarity scoring functions which natively operate on our chosen representations.

Symmetry-preserving generation of molecules in 3D.
Many generative models for small molecules have been developed to sample chemical space (Gómez-Bombarelli et al., ; Segler et al., ; Jin et al., ; Vignac et al., ; Bilodeau et al., ; Anstine & Isayev, ). Our work is most related to 3D approaches which directly generate molecules in specific conformations. These works usually preserve translational and rotational symmetries of atomistic systems by generating structures with E(3)- or SE(3)-invariant internal coordinates (Gebauer et al., ; Luo & Ji, ; Roney et al., ), or more recently, by generating Euclidean coordinates with equivariant networks (Vignac et al., ; Peng et al., ; Irwin et al., ). We generate molecular structures with an equivariant DDPM (Ho et al., ; Hoogeboom et al., ), a strategy that has seen use in protein-conditioned ligand and linker design (Igashov et al., ; Schneuing et al., ). Besides molecular structures, we jointly diffuse/denoise explicit representations of the molecule’s shape, ESP surface, and pharmacophores with their vector directions. To do so, we employ spherical harmonics-based SE(3)-equivariant Euclidean Neural Networks (E3NNs) (Geiger & Smidt, ) to encode/decode coordinates, vectors, and scalar features across our heterogeneous representations.

3D interaction-aware molecular generation.
Prior works have partially explored shape- and pharmacophore-aware molecular generative design. Multiple methods generate either (1D) SMILES or (2D) molecular graphs conditioned on 3D representations of target shapes and/or pharmacophores (Skalic et al., ; Imrie et al., ; Zhu et al., ; Xie et al., ), or use 3D similarity scores to fine-tune unconditional generative models (Papadopoulos et al., ; Neeser et al., ), which we include as a baseline (App. A.1). Since these generative models do not directly predict 3D structures, they require conformer generation as a post-processing step. On the other hand, Adams & Coley () found that generating structures natively in 3D yields more chemically diverse molecules with higher 3D shape-similarity compared to a competing shape-conditioned 1D approach. Multiple other methods have since been developed for shape-conditioned 3D generation (Chen et al., ; Lin et al., ; Le et al., ).
Regarding pharmacophores, Ziv et al. () applied a pretrained DDPM to inpaint 3D molecules given fixed N and O atoms, constrained to be hydrogen bond donors and acceptors, respectively. Yet, they neglect other HBD/HBA definitions and ignore important non-covalent interactions like aromatic $\pi$-$\pi$ interactions, hydrophobic effects, and halogen bonding.
Electrostatics-aware molecular generation has been considered only by Bolcato et al. (),
via exchanging pre-enumerated chemical fragments with similar electrostatics. But, they do not conditionally generate 3D molecules given a global ESP surface.
Our work unifies and extends prior work on interaction-aware 3D generative design by comprehensively modeling shape, electrostatics, and arbitrary pharmacophores (including their directionality) in a single general framework.

3D generative structure-based drug design (SBDD).
Tangential to our work are models that generate ligands inside a protein pocket
by training on protein-ligand complexes with (Lee et al., ; Sako et al., ; Zhung et al., ; Wang et al., ; Huang et al., ) or without (Peng et al., ; Schneuing et al., ; Guan et al., )
explicit encodings of their interactions. In contrast, we consider context-free molecular design. By modeling the interaction profiles of ligands only, ShEPhERD is more general than (yet still applicable to) SBDD.
We also retain the freedom to train on arbitrary chemical spaces that may be larger, denser, or more diverse than the ligands in the PDB.

## 3 Methodology

Figure: Figure 2: (Top) Visualization of our 3D point-cloud similarity scoring functions, which evaluate surface/shape, electrostatic, and pharmacophore similarity via weighted Gaussian overlaps. (Middle) ShEPhERD’s denoising network architecture, which uses SE(3)-equivariant neural networks to (1) embed noisy input 3D molecules and their interaction profiles into scalar and vector node features, (2) jointly interact the node features, and (3) denoise the input states. (Bottom) ShEPhERD uses inpainting to sample chemically diverse 3D molecules that exhibit desired interaction profiles.
Refer to caption: extracted/5947291/figures/shepherd_methods_v3.jpg

We seek to sample chemically diverse molecular structures from 3D similarity-constrained chemical space, where the constraints are defined by the 3D shape, electrostatics, and/or pharmacophores (“interaction profiles”) of a reference 3D molecular system. To do so, we develop ShEPhERD, a DDPM that learns the joint distribution over 3D molecular graphs and their interaction profiles. At inference, we sample from specific interaction-conditioned distributions over 3D chemical space via inpainting.
Below, we formally define our chosen representations for 3D molecules and their interaction profiles, and also define the 3D similarity scoring functions used in our evaluations. We then detail our joint DDPM and sampling protocols. Fig. 2 visualizes our overall methodology.

### 3.1 Defining representations of molecules and their interaction profiles

3D molecular graph. We define an organic molecule with $n_{1}$ atoms in a specific conformation as a 3D molecular graph $\bm{x}_{1}=(\bm{a},\bm{f},\bm{C},\bm{B})$ where $\bm{a}\in\{\text{H},\text{C},\text{O},...\}^{n_{1}}$ lists the atom types ($N_{a}$ options); $\bm{f}\in\{-2,-1,0,1,2\}^{n_{1}}$ specifies the formal charge of each atom; $\bm{C}\in\mathbb{R}^{n_{1}\times 3}$ specifies the atomic coordinates; and $\bm{B}\in\{0,1,2,3,1.5\}^{{n_{1}}\times{n_{1}}}$ is the adjacency matrix specifying the covalent bond order between pairs of atoms. For modeling purposes, the categorical variables $\bm{a},\bm{f},\bm{B}$ are represented as one-hot continuous features so that $\bm{a}\in\mathbb{R}^{{n_{1}}\times N_{a}}$, $\bm{f}\in\mathbb{R}^{{n_{1}}\times 5}$, and $\bm{B}\in\mathbb{R}^{{n_{1}}\times{n_{1}}\times 5}$.

Shape. We define the shape $\bm{x}_{2}=\bm{S}_{2}\in\mathbb{R}^{n_{2}\times 3}$ as a point cloud with $n_{2}$ points sampled from the solvent-accessible surface of $\bm{x}_{1}$. $n_{2}$ is fixed regardless of $n_{1}$. We use surface points to better decouple the shape representation from the molecule’s exact coordinates or 2D graph (App. A.5).

Electrostatic potential (ESP) surface. We represent the ESP surface $\bm{x}_{3}=(\bm{S}_{3},\bm{v})$ by the Coulombic potential at each point in a surface point cloud $\bm{S}_{3}\in\mathbb{R}^{n_{3}\times 3}$, with $\bm{v}\in\mathbb{R}^{n_{3}}$.
We compute $\bm{v}$ from the per-atom partial charges of $\bm{x}_{1}$, which are computed via semi-empirical DFT (App. A.5).

Pharmacophores. The set of $n_{4}$ pharmacophores are represented as $\bm{x}_{4}=(\bm{p},\bm{P},\bm{V})$ where $\bm{p}\in\mathbb{R}^{n_{4}\times N_{p}}$ is a matrix of one-hot encodings of the $N_{p}$ pharmacophore types; $\bm{P}\in\mathbb{R}^{n_{4}\times 3}$ is the pharmacophores’ coordinates; and $\bm{V}\in\{\mathbb{S}^{2},\bm{0}\}^{n_{4}}$ contains the relative unit or zero vectors specifying the directionality of each pharmacophore.
The directional pharmacophores include hydrogen bond acceptors (HBA) and donors (HBD), aromatic rings, and halogen-carbon bonds. Directionless pharmacophores ($\bm{V}[k]=\bm{0}$) include hydrophobes, anions, cations, and zinc binders.
App. A.5 details how we extract pharmacophores from a molecule by adapting known SMARTS patterns.

### 3.2 Shape, electrostatic, and pharmacophore similarity scoring functions

To formulate our scoring functions, we first define a point cloud $\bm{Q}\in\mathbb{R}^{n_{Q}\times 3}$ where each point $\bm{r}_{k}$ is assigned an isotropic Gaussian in $\mathbb{R}^{3}$ (Grant & Pickup, ; Grant et al., ). We measure the Tanimoto similarity $\text{sim}\in[0,1]$ between two point clouds $\bm{Q}_{A}$ and $\bm{Q}_{B}$ using first order Gaussian overlap $O_{A,B}=\sum_{a\in\bm{Q}_{A}}\sum_{b\in\bm{Q}_{B}}w_{a,b}(\frac{\pi}{2\alpha})
^{\frac{3}{2}}\exp{(-\frac{\alpha}{2}\|\bm{r}_{a}-\bm{r}_{b}\|^{2})}$ and $\text{sim}^{*}(\bm{Q}_{A},\bm{Q}_{B})=\frac{O_{A,B}}{O_{A,A}+O_{B,B}-O_{A,B}}$ where $\alpha$ is a Gaussian width and $w_{a,b}$ is a weighting factor.
Note that 3D similarities are sensitive to SE(3) transformations of $\bm{Q}_{A}$ with respect to $\bm{Q}_{B}$. We characterize the similarity at their optimal alignment by $\text{sim}(\bm{Q}_{A},\bm{Q}_{B})=\max_{\bm{R,t}}{\text{sim}^{*}(\bm{RQ}_{A}^{
T}+\bm{t},\bm{Q}_{B}})$ where $\bm{R}\in SO(3)$ and $\bm{t}\in T(3)$. We align using automatic differentiation. Note that $n_{Q_{A}}$ need not equal $n_{Q_{B}}$.

Shape scoring. The volumetric similarity between two atomic point clouds $\bm{C}_{A}$ and $\bm{C}_{B}$ is defined as $\text{sim}^{*}_{\text{vol}}(\bm{C}_{A},\bm{C}_{B})$ with $w_{a,b}=2.7$ and $\alpha=0.81$ (Adams & Coley, ).
We newly define the surface similarity between two surfaces $\bm{S}_{A}$ and $\bm{S}_{B}$ as $\text{sim}^{*}_{\text{surf}}(\bm{S}_{A},\bm{S}_{B})$ with $w_{a,b}=1$, and $\alpha=\Psi(n_{2})$. Here, $\Psi$ is a function fitted to $\text{sim}^{*}_{\text{vol}}$ depending on the choice of $n_{2}$ (App. A.6).

Electrostatic scoring. We define the similarity between two electrostatic potential surfaces $\bm{x}_{3,A}$ and $\bm{x}_{3,B}$ as $\text{sim}_{\text{ESP}}^{*}(\bm{x}_{3,A},\bm{x}_{3,B})$ with $w_{a,b}=\exp{(-\frac{\|\bm{v}_{A}[a]-\bm{v}_{B}[b]\|^{2}}{\lambda})}$, $\alpha=\Psi(n_{3})$, and $\lambda=\frac{0.3}{(4\pi\epsilon_{0})^{2}}$. $\epsilon_{0}$ is the permittivity of vacuum with units $e^{2}(\text{eV}\cdot\text{\AA})^{-1}$.
Inspired by Hodgkin & Richards () and Good (), we use the difference between electrostatic potentials to increase sensitivity to their respective magnitudes (i.e., rather than simply comparing signs).

Pharmacophore scoring. We define the similarity between two sets of pharmacophores $\bm{x}_{4,A}$ and $\bm{x}_{4,B}$ with
$\text{sim}_{\text{pharm}}^{*}(\bm{x}_{4,A},\bm{x}_{4,B})=\frac{\sum_{m\in
\mathcal{M}}O_{A,B;m}}{\sum_{m\in M}O_{A,A;m}+O_{B,B;m}-O_{A,B;m}}$. Here, $\mathcal{M}$ is the set of all pharmacophore types ($|\mathcal{M}|=N_{p}$). $w_{a,b}=1$ if $m$ is non-directional, or a scaling of vector cosine similarity $w_{a,b;m}=\frac{\bm{V}[a]_{m}^{\top}\bm{V}[b]_{m}+2}{3}$ if $m$ is directional (Wahl, ). $\alpha_{m}=\Omega(m)$ where $\Omega$ maps each pharmacophore type to a Gaussian width (App. A.6). We take the absolute value of $\bm{V}[a]_{m}^{\top}\bm{V}[b]_{m}$ for aromatic groups as we assume their $\pi$ interaction effects are symmetric across their plane.

### 3.3 Joint diffusion of molecules and their interaction profiles with ShEPhERD

ShEPhERD follows the DDPM paradigm (Ho et al., ; Sohl-Dickstein et al., ) to decompose the joint distribution over the tuple $\bm{X}=(\bm{x}_{1},\bm{x}_{2},\bm{x}_{3},\bm{x}_{4})$ as $P_{\text{data}}(\bm{X}):=P(\bm{X}^{(0)})=P(\bm{X}^{(T)})\prod_{t=1}^{T}P(\bm{X
}^{(t-1)}|\bm{X}^{(t)})$, where $P_{\text{data}}$ is the data distribution, $P(\bm{X}^{(T)})$ is (roughly) a Gaussian prior, and $P(\bm{X}^{(t-1)}|\bm{X}^{(t)})$ are Markov transition distributions learnt by a neural network.
This network is trained to reverse a forward-noising process $P(\bm{X}^{(t)}|\bm{X}^{(t-1)})=N(\alpha_{t}\bm{X}^{(t-1)},\sigma_{t}^{2}\bm{I})$ which gradually corrupts data $\bm{X}$ into Gaussian noise $\bm{X}^{(T)}$ according to a variance preserving noise schedule given by $\sigma_{t}^{2}$ and $\alpha_{t}=\sqrt{1-\sigma_{t}^{2}}$ for $t=1,...,T$. See Ho et al. () for preliminaries on DDPMs. Here, we describe the forward and reverse processes of ShEPhERD’s joint DDPM.

Forward noising processes. We follow Hoogeboom et al. () to forward-noise the 3D molecule $\bm{x}_{1}=(\bm{a},\bm{f},\bm{C},\bm{B})$. For $\bm{a}\in\mathbb{R}^{n_{1}\times N_{a}}$, we use Gaussian noising where $\bm{a}^{(t)}=\alpha_{t}\bm{a}^{(t-1)}+\sigma_{t}\bm{\epsilon}$ for $\bm{\epsilon}\sim N(\bm{0},\bm{I})$. The processes for $\bm{f}$ and $\bm{B}$ are similar, but we symmetrize the upper/lower triangles of $\bm{B}^{(t)}$. We apply isotropic noise to $\bm{C}\in\mathbb{R}^{{n_{1}}\times 3}$, but center the noise at $\bm{0}$ to ensure translational invariance: $\bm{C}^{(t)}=\alpha_{t}\bm{C}^{(t-1)}+\sigma_{t}(\bm{\epsilon}-\frac{1}{n_{1}}
\sum_{k=1}^{n_{1}}\bm{\epsilon}[k])$ for $\bm{\epsilon}[k]\sim N(\bm{0},\bm{I}_{3})$ and $\bm{\epsilon}\in\mathbb{R}^{n_{1}\times 3}$.

For the molecular shape $\bm{x}_{2}=\bm{S}_{2}\in\mathbb{R}^{{n_{2}}\times 3}$, we also forward-noise with isotropic Gaussian noise: $\bm{S}_{2}^{(t)}=\alpha_{t}\bm{S}_{2}^{(t-1)}+\sigma_{t}\bm{\epsilon}$ for $\bm{\epsilon}\sim N(\bm{0},\bm{I})$. We do not subtract the noise’s center of mass (COM), though; this ensures that the model can learn to denoise $\bm{x}_{2}$ such that it is centered with respect to $\bm{x}_{1}$. For the ESP surface $\bm{x}_{3}=(\bm{S}_{3},\bm{v})$, we forward-noise the surface $\bm{S}_{3}\in\mathbb{R}^{n_{3}\times 3}$ in the same manner as $\bm{S}_{2}$. We forward-noise $\bm{v}_{3}\in\mathbb{R}^{n_{3}}$ in the typical way: $\bm{v}^{(t)}=\alpha_{t}\bm{v}^{(t-1)}+\sigma_{t}\bm{\epsilon}$ for $\bm{\epsilon}\sim N(\bm{0},\bm{I})$.

For the pharmacophores $\bm{x}_{4}=(\bm{p},\bm{P},\bm{V})$, we forward-noise their types $\bm{p}$ just like the atom types $\bm{a}$, and their positions $\bm{P}$ just like the shape $\bm{S}_{2}$. Diffusing the vector directions $\bm{V}\in\{\mathbb{S}^{2},\bm{0}\}^{n_{4}}$ is complicated since some pharmacophores are directionless. To unify their treatment, we interpret the pharmacophore vectors as Euclidean points in $\mathbb{R}^{3}$ (e.g., only noiseless vectors have norm 1.0 or 0.0). We then forward-noise the vectors like any point cloud: $\bm{V}^{(t)}=\alpha_{t}\bm{V}^{(t-1)}+\sigma_{t}\bm{\epsilon}$ for $\bm{\epsilon}\sim N(\bm{0},\bm{I})$.

Whereas the above processes describe the single-step forward transition distributions, note that we can efficiently sample noised structures given any time horizon. For instance, we may directly sample
$\bm{a}^{(t)}=\overline{\alpha}_{t}\bm{a}^{(0)}+\overline{\sigma}_{t}\bm{\epsilon}$ for $\bm{\epsilon}\sim N(\bm{0},\bm{I})$, where $\overline{\alpha}_{t}=\prod_{s=1}^{t}\alpha_{s}$ and $\overline{\sigma}_{t}=\sqrt{1-\overline{\alpha}^{2}_{t}}$.

Reverse denoising process. Starting from any $\bm{X}^{(t)}$ (but typically pure noise $\bm{X}^{(T)}\sim N(\bm{0},\bm{I})$), the DDPM iteratively denoises $\bm{X}^{(t)}$ by stochastically interpolating towards a predicted clean structure $\hat{\bm{X}}^{(t)}\approx\bm{X}^{(0)}\sim P_{\text{data}}(\bm{X}^{(0)}|\bm{X}^
{(t)})$
resembling true samples from the data distribution. We follow the DDPM formulation where rather than predicting $\bm{X}^{(0)}$ directly, the network predicts the true noise $\bm{\hat{\epsilon}}^{(t)}\approx\bm{{\epsilon}}$ that, when applied to data $\bm{X}^{(0)}$, yields $\bm{X}^{(t)}=\overline{\alpha}_{t}\bm{X}^{(0)}+\overline{\sigma}_{t}\bm{\epsilon}$. In this case, the single-step denoising update can be derived as:
$\bm{X}^{(t-1)}=\frac{1}{\alpha_{t}}\bm{X}^{(t)}-\frac{\sigma_{t}^{2}}{\alpha_{
t}\overline{\sigma}_{t}}\bm{\hat{\epsilon}}^{(t)}+\frac{\sigma_{t}\overline{
\sigma}_{t-1}}{\overline{\sigma}_{t}}\bm{\epsilon}^{\prime}$,
where the additional noise $\bm{\epsilon}^{\prime}\sim N(\bm{0},\bm{I})$ (set to $\bm{\epsilon}^{\prime}=\bm{0}$ for $t=1$) makes each denoising step stochastic.

ShEPhERD employs a single denoising network $\bm{\eta}$ that is trained to jointly predict the noises
$\bm{\hat{\epsilon}}_{1}^{(t)},\bm{\hat{\epsilon}}_{2}^{(t)},\bm{\hat{\epsilon}
}_{3}^{(t)},\bm{\hat{\epsilon}}_{4}^{(t)}=\bm{\eta}(\bm{x}_{1}^{(t)},\bm{x}_{2
}^{(t)},\bm{x}_{3}^{(t)},\bm{x}_{4}^{(t)},t)$
,
where $\bm{\hat{\epsilon}}_{1}^{(t)}=(\bm{\hat{\epsilon}_{a}}^{(t)},\bm{\hat{\epsilon
}_{f}}^{(t)},\bm{\hat{\epsilon}_{C}}^{(t)},\bm{\hat{\epsilon}_{B}}^{(t)})$,
$\bm{\hat{\epsilon}}_{2}^{(t)}=(\bm{\hat{\epsilon}}_{\bm{S}_{2}}^{(t)})$,
$\bm{\hat{\epsilon}}_{3}^{(t)}=(\bm{\hat{\epsilon}}_{\bm{S}_{3}}^{(t)},\bm{\hat
{\epsilon}}_{\bm{v}}^{(t)})$, and
$\bm{\hat{\epsilon}}_{4}^{(t)}=(\bm{\hat{\epsilon}_{p}}^{(t)},\bm{\hat{\epsilon
}_{P}}^{(t)},\bm{\hat{\epsilon}_{V}}^{(t)})$. At inference, we jointly apply the denoising updates to obtain $\bm{X}^{(t-1)}$. When computing $\bm{x}_{1}^{(t-1)}$, we remove the COM from $\bm{\hat{\epsilon}_{C}}^{(t)}$ and the extra noise $\bm{{\epsilon}_{C}}^{\prime}$.

The forward and reverse processes of our joint DDPM are designed to be straightforward to make ShEPhERD flexible:
One may freely adjust the exact representations of the shape, electrostatics, or pharmacophores as long as they can be represented as a point-cloud with one-hot, scalar, and/or vector attributes.
One can also directly model specific marginal distributions (e.g., $P(\bm{x}_{1},\bm{x}_{2})$, $P(\bm{x}_{1},\bm{x}_{3})$, or $P(\bm{x}_{1},\bm{x}_{4})$) by simply modeling a subset of the variables $\{\bm{x}_{2},\bm{x}_{3},\bm{x}_{4}\}$.
Finally, ShEPhERD can be easily extended to model other structural features or interaction profiles beyond those considered here by defining their explicit structural representations and forward/reverse processes.

Denoising network design. We design ShEPhERD’s denoising network $\bm{\eta}$ to satisfy three criteria:

- •
Symmetry-preserving: The noise predictions $\bm{\hat{\epsilon}}_{1}^{(t)},\bm{\hat{\epsilon}}_{2}^{(t)},\bm{\hat{\epsilon}
}_{3}^{(t)},\bm{\hat{\epsilon}}_{4}^{(t)}$ are T(3)-invariant and SO(3)-equivariant with respect to global SE(3)-transformations of $\bm{X}^{(t)}$ in order to (1) efficiently preserve molecular symmetries, and (2) ensure $\bm{x}_{1}^{(t)},\bm{x}_{2}^{(t)},\bm{x}_{3}^{(t)},$ and $\bm{x}_{4}^{(t)}$ remain aligned during denoising.
- •
Expressive: $\bm{\eta}$ captures both local and global relationships between $\bm{x}_{1}^{(t)},\bm{x}_{2}^{(t)},\bm{x}_{3}^{(t)},$ and $\bm{x}_{4}^{(t)}$.
- •
General: To promote applications across chemical design, $\bm{\eta}$ accommodates other definitions of shape/electrostatics/pharmacophores and permits incorporating other structural interactions, too.

To achieve these criteria, we design $\bm{\eta}$ to have three components: (1) a set of embedding modules which equivariantly encode the heterogeneous $\bm{x}_{i}$ into a uniform sets of latent $l{=}0$ (scalar) and $l{=}1$ (vector) node representations; (2) a joint module which locally and globally interacts these latent node representations; and (3) a set of denoising modules which predict $\bm{\hat{\epsilon}}_{i}$ for each $\bm{x}_{i}$.

The embedding modules use SE(3)-equivariant E3NNs (we choose to use expressive EquiformerV2 modules (Liao et al., )) to individually encode the 3D structures of $\bm{x}_{1}^{(t)},\bm{x}_{2}^{(t)},\bm{x}_{3}^{(t)},$ and $\bm{x}_{4}^{(t)}$ into latent codes $(\bm{z}^{(t)}_{i},\bm{\tilde{z}}^{(t)}_{i})=\phi_{i}(\bm{x}_{i}^{(t)},t)\
\forall\ i\in[1,2,3,4]$. Here, $\bm{z}_{i}\in\mathbb{R}^{n_{i}\times d}$ are invariant scalar representations of each node (e.g., atom, point, or pharmacophore), and $\bm{\tilde{z}}_{i}\in\mathbb{R}^{n_{i}\times 3\times d}$ are equivariant vector representations of each node. To make each system $\bm{x}_{i}$ sensitive to relative translations with respect to $\bm{x}_{1}^{(t)}$, we also include a virtual node that is positioned at the center of mass of $\bm{x}_{1}^{(t)}$, and which remains unnoised. Prior to 3D message passing with the E3NNs, scalar atom/point/pharmacophore features (e.g., $\bm{a}^{(t)}$, $\bm{f}^{(t)}$, $\bm{v}^{(t)}$, $\bm{p}^{(t)}$) are embedded into $l{=}0$ node features, and vector features (e.g., $\bm{V}^{(t)}$) are directly assigned as $l{=}1$ node features. The pairwise bond representations $\bm{B}^{(t)}$ for $\bm{x}_{1}^{(t)}$ are also embedded as $l{=}0$ edge attributes. Finally, for each $\phi_{i}$, we embed sinusoidal positional encodings of the time step $t$ and add these to all the $l{=}0$ node embeddings.

The joint module consists of two steps to locally and globally interact the joint latent variables:

(1) We collate the coordinates of $\bm{x}_{1}^{(t)},\bm{x}_{2}^{(t)},\bm{x}_{3}^{(t)},$ and $\bm{x}_{4}^{(t)}$ to form a heterogeneous 3D graph where the nodes are attributed with their corresponding latent features $(\bm{z}^{(t)}_{i},\bm{\tilde{z}}^{(t)}_{i})\ \forall\ i\in[1,2,3,4]$. We then encode this heterogeneous 3D graph with another E3NN (EquiformerV2) module $\phi_{joint}^{local}$ and residually update the nodes’ latent features: $(\bm{z}^{(t)}_{i},\bm{\tilde{z}}^{(t)}_{i})\mathrel{{+}{=}}\phi_{joint}^{local
}\bigg{(}(\bm{x}_{i}^{(t)},\bm{z}^{(t)}_{i},\bm{\tilde{z}}^{(t)}_{i})\ \forall
\ i\in[1,2,3,4]\bigg{)}$.

(2) We sum-pool the updated $l{=}1$ node features across each sub-graph, concatenate, and then embed with an equivariant feed-forward network $\phi_{joint}^{global}$ to obtain a global $l{=}1$ code describing the overall system: $\bm{\tilde{z}}_{joint}^{(t)}=\phi_{joint}^{global}\bigg{(}\text{Cat}\bigg{[}
\bigg{(}\sum_{k=1}^{n_{i}}\bm{\tilde{z}}^{(t)}_{i}[k]\bigg{)}\ \forall\ i\in[1
,2,3,4]\bigg{]}\bigg{)}$. We then apply equivariant tensor products between $\bm{\tilde{z}}_{joint}^{(t)}\in\mathbb{R}^{d\times 3}$ and an $l{=}0$ embedding of $t$. This yields $l{=}0$ and $l{=}1$ global latent features $(\bm{z}_{joint}^{(t)},\bm{\tilde{z}}_{joint}^{(t)})$, which are residually added to the node representations $(\bm{z}^{(t)}_{i},\bm{\tilde{z}}^{(t)}_{i})$.

The denoising modules predict the noises $\bm{\hat{\epsilon}}_{i}$ from the node-level features $(\bm{z}^{(t)}_{i},\bm{\tilde{z}}^{(t)}_{i})$.
For $\bm{x}_{2}$, $\bm{x}_{3}$, and $\bm{x}_{4}$, the scalar noises ($\bm{\hat{\epsilon}}_{\bm{v}}^{(t)}$, $\bm{\hat{\epsilon}}_{\bm{p}}^{(t)}$) are predicted from the corresponding ($l{=}0$) $\bm{z}^{(t)}_{i}$ codes using multi-layer perceptrons (MLPs),
whereas the vector noises ($\bm{\hat{\epsilon}}_{\bm{S}_{2}}^{(t)}$, $\bm{\hat{\epsilon}}_{\bm{S}_{3}}^{(t)}$, $\bm{\hat{\epsilon}_{P}}^{(t)}$, $\bm{\hat{\epsilon}_{V}}^{(t)}$) are predicted from the $(l{=}1)$ $\bm{\tilde{z}}^{(t)}_{i}$ codes using equivariant feed-forward networks (“E3NN-style” coordinate predictions).
For $\bm{x}_{1}$, $\bm{\hat{\epsilon}_{a}}^{(t)}$ and $\bm{\hat{\epsilon}_{f}}^{(t)}$ are predicted from $\bm{z}^{(t)}_{1}$ with simple MLPs. $\bm{\hat{\epsilon}_{B}}^{(t)}$ is predicted from pairs $(\bm{{z}}^{(t)}_{1}[k],\bm{{z}}^{(t)}_{1}[j])$ using a permutation-invariant MLP (App. A.7.1). $\bm{\hat{\epsilon}_{C}}^{(t)}$ is predicted from $(\bm{z}^{(t)}_{1},\bm{\tilde{z}}^{(t)}_{1})$ using E3NN-style and EGNN-style (Satorras et al., ) coordinate predictions (App. A.7.1).

Training. We train the denoising network $\bm{\eta}$ with unweighted L2 regression losses $l_{i}=||\bm{\hat{\epsilon}}_{i}-\bm{{\epsilon}}_{i}||^{2}$ between the predicted and true noises.
Importantly, our framework permits us to train models which directly learn certain marginal distributions. In our experiments, we train models to learn $P(\bm{x}_{1})$, $P(\bm{x}_{1},\bm{x}_{2})$, $P(\bm{x}_{1},\bm{x}_{3})$, $P(\bm{x}_{1},\bm{x}_{4})$, and $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$. Note that since $\bm{x}_{3}$ defines an (attributed) surface $\bm{S}_{3}$, jointly modeling $(\bm{x}_{2},\bm{x}_{3})$ is redundant; $\bm{x}_{3}$ implicitly models $\bm{x}_{2}$. App. A.7 provides details on training protocols, choice of noise schedules, feature scaling, and model hyperparameters.

Sampling. For unconditional generation, we first sample $\bm{X}^{(T)}$ from isotropic Gaussian noise, and then denoise for $T$ steps to sample $(\bm{x}_{1}^{(0)},\bm{x}_{2}^{(0)},\bm{x}_{3}^{(0)},\bm{x}_{4}^{(0)})$.
We then argmax $\bm{a}^{(0)}$, $\bm{f}^{(0)}$, $\bm{B}^{(0)}$, and $\bm{p}^{(0)}$ to obtain discrete atom/bond/pharm. types, and round each $\bm{V}^{(0)}[k]$ to have norm 1.0 or 0.0.
To help break the spherical symmetry of the coordinates of $\bm{X}^{(t)}$ at early time steps, we find it useful to strategically add extra noise to the reverse process (App. A.8). For conditional generation, we use inpainting to sample $\bm{x}_{1}^{(0)}$ conditioned on target interaction profiles. Namely, we first simulate the forward-noising of the target profiles $(\bm{x}^{*}_{2},\bm{x}^{*}_{3},\bm{x}^{*}_{4})\rightarrow(\bm{x}^{*}_{2},\bm{x
}^{*}_{3},\bm{x}^{*}_{4})^{(t)}\ \forall\ t\in[1,T]$. Then, during the reverse denoising process, we replace the ShEPhERD-denoised $(\bm{x}_{2}^{(t)},\bm{x}_{3}^{(t)},\bm{x}_{4}^{(t)})$ with the noisy target $(\bm{x}^{*}_{2},\bm{x}^{*}_{3},\bm{x}^{*}_{4})^{(t)}$. Like other molecule DDPMs, we must specify $n_{1}$ and $n_{4}$ for each sample.

## 4 Experiments

We train and evaluate ShEPhERD using two new datasets. Our first dataset (ShEPhERD-GDB17) contains 2.8M molecules sampled from medicinally-relevant subsets of GDB17 (Ruddigkeit et al., ; Awale et al., ; Bühlmann & Reymond, ). Each molecule contains ${\leq}$17 non-hydrogen atoms with element types in {H, C, N, O, S, F, Cl, Br, I}, and includes one conformation optimized with GFN2-xTB in the gas phase (Bannwarth et al., ). Our second dataset (ShEPhERD-MOSES-aq) contains 1.6M drug-like molecules from MOSES (Polykovskiy et al., ) with up to 27 non-hydrogen atoms. Each molecule contains one conformation optimized with GFN2-xTB in implicit water. In all experiments, hydrogens are treated explicitly. Whereas we use ShEPhERD-GDB17 to evaluate ShEPhERD in straightforward unconditional and conditional generation settings, we use ShEPhERD-MOSES-aq to challenge ShEPhERD to design drug-like analogues of natural products, to diversify bioactive hits, and to merge fragments from a fragment screen.

Figure: Figure 3: (Left) Self-consistency of jointly generated 3D molecules and their shapes, ESP surfaces, or pharmacophores, as assessed via 3D similarity between the generated vs. true interaction profiles of the generated molecules. Shape and ESP consistency are bounded due to randomness in surface sampling. (Right) Distributions of 3D interaction similarities between ShEPhERD-generated or dataset-sampled 3D molecules (post-relaxation and realignment) and 100 target molecules from ShEPhERD-GDB17, including the top-1 scores given 20 samples per target. ShEPhERD generates 3D molecules with low graph similarity to the target molecule, and which have stable geometries as measured by heavy-atom RMSD upon xTB-relaxation. Also shown are top-scoring samples overlaid on their target profiles (surfaces are upsampled for visualization), labeled with 3D similarity scores.
Refer to caption: extracted/5947291/figures/GDB_results_v6.jpg

Unconditional joint generation. We first evaluate ShEPhERD’s ability to jointly generate 3D molecules and their interaction profiles in a self-consistent way. Namely, a well-trained model that learns the joint distribution should generate interaction profile(s) that match the true interaction profile(s) of the generated molecule. Fig. 3 reports distributions of the 3D similarity between generated and true interaction profiles across 1000 samples (with $n_{1}\in[11,60]$, $n_{4}\sim P_{\text{data}}(n_{4}|n_{1})$) from models trained on ShEPhERD-GDB17 to learn $P(\bm{x}_{1},\bm{x}_{2})$, $P(\bm{x}_{1},\bm{x}_{3})$, or $P(\bm{x}_{1},\bm{x}_{4})$.(^2^22We compute the true profiles of generated molecules after they are relaxed with xTB and realigned to the unrelaxed structures by minimizing heavy-atom RMSD; this avoids optimistically scoring generated molecules that are highly strained. Nonetheless, the average RMSD upon relaxation is ${<}$0.1 Å for all unconditional models trained on ShEPhERD-GDB17. App. A.4 reports such unconditional metrics (validity, novelty, etc.).)
When we compare against the (optimally aligned) similarities $\text{sim}_{\text{surf}}$, $\text{sim}_{\text{ESP}}$, and $\text{sim}_{\text{pharm}}$ between the true profiles of the generated molecules and those of random molecules from the dataset, we confirm that ShEPhERD’s generated profiles have substantially enriched similarities to the true profiles in all cases. Interestingly, ShEPhERD is more self-consistent when generating shapes or ESP surfaces than when generating pharmacophores. We partially attribute this performance disparity to the discrete nature of the pharmacophore representations and the requirement of specifying $n_{4}$; ShEPhERD occasionally generates molecules that have more true pharmacophores than generated pharmacophores, as demonstrated by the samples shown in Fig. 3 (i.e., some have missing HBAs).

Interaction-conditioned generation. We now evaluate ShEPhERD’s ability to sample from interaction-conditioned chemical space. To do so, we reuse the same $P(\bm{x}_{1},\bm{x}_{2})$, $P(\bm{x}_{1},\bm{x}_{3})$, and $P(\bm{x}_{1},\bm{x}_{4})$ models trained on ShEPhERD-GDB17, but use inpainting to sample from the interaction-conditioned distributions $P(\bm{x}_{1}|\bm{x}_{2})$, $P(\bm{x}_{1}|\bm{x}_{3})$, and $P(\bm{x}_{1}|\bm{x}_{4})$. Specifically, we extract the true interaction profiles from 100 random target molecules (held out from training), and use ShEPhERD to inpaint new 3D molecular structures given these target profiles. After generating 20 structures per target profile and discarding (without replacement) any invalid structures (App. A.3) or those with 2D graph similarities ${>}$0.3 to the target molecule, we relax the generated structures with xTB and compute their optimally-realigned 3D similarities to the target. Fig. 3 plots the distributions of 3D similarity scores between all valid (sample, target) pairs as well as the top-1 scores amongst the 20 samples per target. We compare against analogous similarity distributions for randomly sampled molecules from the dataset. Overall, ShEPhERD generates structures with very low graph similarity (${\geq}94\%$ of valid samples have graph similarity ${<}$0.2) but significantly enriched 3D similarities to the target, for all versions of the model. Qualitatively, ShEPhERD can generate molecular structures that satisfy very detailed target interactions, including multiple directional pharmacophores (Fig. 3).

Figure: Figure 4: (Left) Examples of ShEPhERD-generated analogues of natural product targets, labeled by SA score, ESP similarity, and pharmacophore similarity to the target. Similarities are computed after xTB-relaxation and ESP-optimal realignment. (Right) Distributions of Vina scores for ${\leq}$500 samples from ShEPhERD when conditioning on the bound or lowest-energy pose of co-crystal PDB ligands across 7 proteins. We compare against the Vina scores of 10K virtually screened molecules from ShEPhERD-MOSES-aq. For 5mo4 and 7l11, we show top-scoring ShEPhERD-generated ligands (conditioned on low-energy poses), and overlay a selection from their top-10 docked poses on the PDB ligands. (Bottom) ShEPhERD’s bioisosteric fragment merging workflow. We extract the ESP surface and pharmacophores of 13 fragments from a fragment screen, and show ShEPhERD-generated ligands with low SA score and high 3D similarities to the fragments’ interaction profiles.
Refer to caption: extracted/5947291/figures/shepherd_demonstrations_v10.jpg

Natural product ligand hopping. Numerous clinically-approved drugs have structures derived from natural products due to their rich skeletal complexity, high 3D character, and wide range of pharmacophores that impart uniquely selective biological function. But, the structural complexity of natural products limits their synthetic tractability. As such, designing synthetically-accessible small-molecule analogues of natural products that mimic their precise 3D interactions is a preeminent task in scaffold/ligand hopping. To imitate this task, we select three complex natural products from COCONUT (Sorokina et al., ), including two large macrocycles and a fused ring system with 9 stereocenters. We then apply ShEPhERD (trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq) to generate drug-like molecules conditioned jointly on the ESP surface and pharmacophores of the lowest-energy conformer of each natural product, again via inpainting. We emphasize that these natural products are out-of-distribution compared to the drug-like molecules contained in ShEPhERD-MOSES-aq. Upon generating 2500 samples from the conditional distribution $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ for each natural product, ShEPhERD identifies small-molecule mimics that attain high ESP and pharmacophore similarity to the natural products (Fig. 4), despite having much simpler chemical structures as assessed via SA score (Ertl & Schuffenhauer, ). Crucially, ShEPhERD generates molecules with higher 3D similarity compared to 2500 molecules sampled from the dataset, and compared to 10K molecules optimized by REINVENT (Blaschke et al., ) (App. A.1).

Bioactive hit diversification. Whereas stucture-based drug design aims to design high-affinity ligands given the structure of the protein target, ligand-based drug design attempts to develop and optimize bioactive hit compounds in the absence of protein information. A common task in ligand-based drug design is diversifying the chemical structures of previously identified bioactive hits (i.e., from a phenotypic experimental screen) through interaction-preserving scaffold hopping, often as a means to reduce toxicity, increase synthesizability, or evade patent restrictions.
We simulate bioactive hit diversification by using ShEPhERD to generate analogues of 7 experimental ligands from the PDB. To evaluate their likelihood of retaining bioactivity, we use Autodock Vina (Trott & Olson, ; Eberhardt et al., ) to dock the generated ligands to their respective proteins, treating the Vina docking scores as a weak surrogate for bioactivity. To best imitate ligand-based design, we condition ShEPhERD on the ESP surface and pharmacophores of the lowest-energy conformer of each PDB ligand, rather than their bound poses (we also simulate this scenario for comparison). We then use inpainting to generate 500 samples from $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$, and dock the valid samples. Fig. 4 shows the distributions of Vina scores for ShEPhERD-generated analogues, compared against a docking screen of 10K random compounds from ShEPhERD-MOSES-aq. Despite having no knowledge of the protein targets, ShEPhERD enriches average and top-10 Vina scores in 5/7 cases. For the top-scoring generated ligands for 5mo4 and 7l11, ShEPhERD generates substantial scaffold hops that lead to diverse 2D graph structures relative to the experimental ligands. Nevertheless, upon docking, the generated molecules still explore poses that closely align with the experimental crystal poses.

Bioisosteric fragment merging. Fragment screening seeks to identify protein-ligand binding modes by analyzing how small chemical fragments bind to a protein of interest. Clusters of protein-bound fragment hits can then be analyzed to design high-affinity ligands. Multiple methods have been developed to link fragment hits into a single ligand containing the original fragments and the new linker. Recently, Wills et al. () showed that merging (not linking) the fragments to form a bioisosteric ligand (which may not contain the exact fragment hits) can diversify ligand hypotheses while preserving the fragments’ important binding interactions.
While ShEPhERD could link fragments, ShEPhERD is uniquely suited to bioisosteric fragment merging as it can condition on aggregate fragment interactions. We use ShEPhERD to merge a set of 13 fragments experimentally identified to bind to the antiviral target EV-D68 3C protease (Lithgo et al., ; Wills et al., ). We extract $n_{4}=27$ pharmacophores by clustering common motifs and selecting interactions identified by Fragalysis (Diamond Light Source, ) (App. A.2).
We also compute an aggregate ESP surface by sampling points from the surface of the overlaid fragments and averaging the fragments’ ESP contributions at each point. We then condition ShEPhERD on these profiles to sample 1000 structures ($n_{1}\in[50,89]$) from $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ via inpainting. Fig. 4 shows samples with SA $\leq 4.0$ that score in the top-10 by combined ESP and pharmacophore similarity. Visually, ShEPhERD generates structures that align well to the fragments and preserve many of their binding interactions, even though $n_{1}$ and $n_{4}$ are significantly out-of-distribution from ShEPhERD-MOSES-aq (App. A.9).

## 5 Conclusion

We introduced ShEPhERD, a new 3D molecular generative model that facilitates interaction-aware chemical design by learning the joint distribution over 3D molecular structures and their shapes, electrostatics, and pharmacophores. Empirically, ShEPhERD can sample chemically diverse molecules with highly enriched interaction-similarity to target structures, as assessed via custom 3D similarity scoring functions. In bioisosteric drug design, ShEPhERD can design small-molecule mimics of complex natural products, diversify bioactive hits while enriching docking scores despite having no knowledge of the protein, and merge fragments from experimental fragment screens into bioisosteric ligands. We anticipate that future work will creatively extend ShEPhERD to other areas of interaction-aware chemical design such as structure-based drug design and organocatalyst design.

#### Reproducibility Statement

Our main text and appendices provide all critical details necessary to understand and reproduce our work, including our training and sampling protocols. To ensure reproducibility, we make our datasets and all training, inference, and evaluation code available on Github at [https://github.com/coleygroup/shepherd](https://github.com/coleygroup/shepherd) and [https://github.com/coleygroup/shepherd-score](https://github.com/coleygroup/shepherd-score).

#### Acknowledgments

This research was supported by the Office of Naval Research under grant number ONR N00014-21-1-2195. This material is based upon work supported by the National Science Foundation Graduate Research Fellowship under Grant No. 2141064. The authors acknowledge the MIT SuperCloud and Lincoln Laboratory Supercomputing Center for providing HPC resources that have contributed to the research results reported within this paper.

## Appendix A Appendix

### A.1 Comparing ShEPhERD to REINVENT and virtual screening for natural product ligand hopping

REINVENT, a state-of-the-art baseline for generative molecular design and optimization, applies a reinforcement learning policy to iteratively update a SMILES recursive neural network (RNN) with a provided reward function (Blaschke et al., ). We applied REINVENT to the three natural product ligand hopping tasks defined in section 4, using a combination of our ESP and pharmacophore 3D similarity scoring functions (section 3.2) as REINVENT’s reward function. We followed the REINVENT implementation in the Practical Molecular Optimization (PMO) benchmark (Gao et al., ). REINVENT was pretrained on the ZINC database (Sterling & Irwin, ) and was deployed with a batch size of 64, $\sigma=500$, an experience replay of 24, and an oracle budget of $10,000$. Training was performed using an Adam optimizer with a learning rate of $5\times 10^{-4}$. Any SMILES which failed during pharmacophore scoring were assigned a score of 0.

Since REINVENT generates molecules in a 1D SMILES representation, we generate up to 5 conformers for each SMILES in order to apply our 3D similarity scoring functions as the reward function. The procedure to compute the reward for a single generated SMILES is as follows: 1) embed 5 conformers with RDKit’s EmbedMultipleConfs function, which uses ETKDG; 2) optimize each conformer with MMFF94 for a maximum of 200 steps; 3) cluster the conformers with Butina clustering using an RMSD threshold of $0.1\text{ \AA}$; 4) relax each remaining conformer with xTB in implicit water; 5) extract the ESP surface $\bm{x}_{3}$ with $n_{3}=400$ and pharmacophore profile $\bm{x}_{4}$ of each relaxed conformer; 6) align each conformer to the target natural product by optimizing our 3D ESP scoring function; 7) calculate the ESP and pharmacophore similarity scores of the ESP-aligned conformers; 8) add the ESP and pharmacophore similarity scores to obtain one combined score per conformer; and 9) take the maximum score across the conformers as the reward.

Figure: Figure 5: Distributions of 3D ESP and pharmacophore similarity to each of three natural product targets for small-molecules sampled via ShEPhERD, REINVENT, or virtual screening (VS). Similarity scores are computed after optimizing the molecular geometry with xTB and aligning the structure to the natural product by maximizing ESP similarity. For ShEPhERD and REINVENT, only valid samples with SA score $<4.5$ are included. Also visualized are examples of top-scoring samples generated by ShEPhERD and REINVENT.
Refer to caption: extracted/5947291/figures/appendix_figures/reinvent_comparisons_v1.jpg

Fig. 5 compares the distributions of ESP and pharmacophore similarity for samples obtained by (1) using ShEPhERD (trained on ShEPhERD-MOSES-aq) to sample 2500 molecules from $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ via inpainting; (2) virtually screening (VS) 2500 random 3D molecules from ShEPhERD-MOSES-aq; and (3) optimizing REINVENT with an oracle/sampling budget of 10,000. Note that REINVENT was pretrained on ZINC, and ShEPhERD’s training set (MOSES-aq) is a small subset of ZINC. Each 3D molecule generated by ShEPhERD was relaxed with xTB prior to realigning the relaxed structure to the natural product (via maximizing ESP similarity) and scoring the ESP and pharamcophore similarity of the aligned pose. 3D molecules sampled from ShEPhERD-MOSES-aq (which are already xTB-relaxed) were directly aligned to the natural product in the same manner. Samples from REINVENT were scored using the procedure outlined in the preceding paragraph. For both ShEPhERD and REINVENT, we only compare valid samples that have SA scores lower than 4.5. This means that although we initially obtain 2500 and 10000 samples from ShEPhERD and REINVENT, respectively, we only compare ${\sim}500$ samples from ShEPhERD against ${\sim}9000$ samples from REINVENT, for each case study. Despite the fewer number of samples, ShEPhERD still finds molecules beneath the SA-score threshold that score higher than molecules optimized by REINVENT, for all three natural products. Both ShEPhERD and REINVENT find much better molecules than those obtained by randomly sampling from the dataset.

We also emphasize that unlike REINVENT, ShEPhERD is not trained to directly optimize 3D similarity scores to each natural product. Indeed, we use the same ShEPhERD model with the same weights for each natural product target. Nor does ShEPhERD employ any inference optimization strategies beyond conditional generation via inpainting. Hence, although ShEPhERD is already quite capable out-of-the-box, we expect ShEPhERD to be able to find even better small-molecule mimics by combining inpainting with other inference-time optimization strategies. For instance, one could use a genetic algorithm which iteratively mutates and evolves the best-scoring samples by partial forward-noising and subsequent denoising.

### A.2 Additional details on experiments

#### A.2.1 Unconditional joint generation on ShEPhERD -GDB17

We define self-consistency as the 3D similarity between the ShEPhERD’s generated interaction profile (i.e., $\bm{x}_{2}$, $\bm{x}_{3}$, or $\bm{x}_{4}$) and the true interaction extracted from the generated molecule $\bm{x}_{1}$ after xTB relaxation and alignment to the unrelaxed pose via minimizing heavy atom RMSD. We align via heavy-atom RMSD (rather than aligning with the scoring functions) for our self-consistency evaluations in order to preserve the natively generated pose of $\bm{x_{1}}$ as much as possible. Since the $P(\bm{x}_{1},\bm{x}_{2})$ and $P(\bm{x}_{1},\bm{x}_{3})$ models employ $n_{2},n_{3}=75$, all self-consistency similarity evaluations also use 75 surface points when extracting surface and ESP interaction profiles from the (xTB-relaxed and RMSD-aligned) generated structures. Note that this differs from our other experiments involving conditional generation, which employ $n_{2},n_{3}=400$ for similarity scoring.

For each valid ShEPhERD-generated molecule, the lower bound of self-consistency is calculated by (globally) aligning a sampled molecule from ShEPhERD-GDB17 to the true interaction profile.
Since surface sampling is stochastic (App. A.5), the self-similarity between two surfaces extracted from the same $\bm{x}_{1}$ may not equal 1.0. Hence, we compare the self-consistency of ShEPhERD’s jointly generated samples to upper bounds for the surface and ESP models, calculated as the expectation of the self-similarity between repeated extractions of the same profile from the same $\bm{x}_{1}$.
We estimate these upper bounds by resampling the true profiles for each of the xTB-relaxed generated molecules 5 times and computing the average similarity between the resampled profiles. The upper bound for pharmacophore self-consistency is 1.0 since pharmacophore extraction is deterministic.

#### A.2.2 Conditional generation on ShEPhERD -GDB17

We use our 3D similarity scoring functions to evaluate the conditional ShEPhERD models. For each ShEPhERD-generated sample that is valid after xTB relaxation, we score the similarity of the relaxed molecule to the target interaction profile after optimal alignment using the relevant scoring function. For surface and ESP similarity, we use $n_{2},n_{3}=400$ to lessen the effect of stochastic surface sampling on the scores – we only notice small deviations in scores (generally ${<}0.03$) due to stochasticity (Fig. 17), which should be considered when comparing any two similarity scores. As a baseline, we also sample 20 molecules from the ShEPhERD-GDB17 dataset for each target and score their similarity after optimally aligning the sampled molecules to the target interaction profile.

#### A.2.3 Natural product ligand hopping

We selected the three natural products in Fig. 4 by manually searching the COCONUT database ([https://coconut.naturalproducts.net/](https://coconut.naturalproducts.net/)) for natural product with complex structures that are meaningfully out of distribution compared to ShEPhERD-MOSES-aq in terms of their molecular size, 3D conformation, number of pharmacophores, and/or stereochemistry. Notably, all selected natural products have SA scores exceeding 6.0 and have intricate 3D shapes.

Each natural product underwent an extensive conformer search and optimization to identify the lowest-energy conformer as evaluated with xTB in implicit water. Up to 1000 conformers of each natural product were initially enumerated with RDKit before undergoing geometry relaxation with xTB. We extracted the interaction profiles of the lowest-energy conformer following App. A.5.

For each natural product, we sampled $2500$ structures from $P(\bm{x}_{1}|\bm{x}^{*}_{3},\bm{x}^{*}_{4})$ by inpainting, where $(\bm{x}^{*}_{3},\bm{x}^{*}_{4})$ are the ESP and pharmacophore interaction profiles of the natural product. We use the version of ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq. We used $n_{3}=75$, specified $n_{4}$ based on the target $\bm{x}^{*}_{4}$, and swept over 25 values of $n_{1}$ in the range $[36,80]$ to force ShEPhERD to generate both small and larger molecules, even though smaller molecules are less likely to occupy the entire volume of the target natural product (and hence less likely to score well). We sampled 100 structures for each choice of $n_{1}$. For the valid samples, we relaxed the generated 3D molecular conformation with xTB, and realigned the relaxed structure to the natural product by optimizing our ESP similarity scoring function after resampling $\bm{x}_{3}$ for both structures at $n_{3}=400$. We used ESP-based alignment instead of pharmacophore-based alignment as ESP surfaces better capture the 3D shapes of the molecules, and because ESP-based alignments are more robust. Given the ESP-aligned relaxed structure, we scored its ESP and pharmacophore similarity to the natural product. We consider both scores in our evaluations and comparisons to REINVENT (App. A.1).

#### A.2.4 Bioactive hit diversification

We used the same PDB co-crystal structures as used in (Zheng et al., ) and use the prepared receptor files .pdbqt AutoDock Vina files provided by Therapeutic Data Commons (TDC) (Huang et al., ). These targets are 1iep (Nagar et al., ; ), 3eml (Jaakola et al., ; ), 3ny8 (Wacker et al., ; ), 4rlu (Li et al., ; Dong et al., ), 4unn (Read et al., ; Naik et al., ), 5mo4 (Cowan-Jacob, ; Wylie et al., ), and 7l11 (Deshmukh et al., ; Zhang et al., ). We use the corresponding co-crystal ligands: STI, ZMA, JRZ, HCC, QZZ, Nilotinib (NIL), and XF1.

For each experimental ligand, the lowest-energy conformer was identified by initially embedding 1000 conformers with RDKit and relaxing each with xTB in implicit water. Then we extract the interaction profiles of the lowest energy conformer following App. A.5. We use the crystal structure pose from the PDB (Berman et al., ) to also extract the interaction profiles and obtain the “PDB pose” conformation. We added hydrogens to these PDB poses using RDKit’s AddHs function with addCoords=True, but did not relax the structures with xTB. The ESP surfaces were computed from the partial charges of the native conformation (with added hydrogens) using a single point xTB calculation in implicit water.

For the two poses of each experimental ligand, we generate 500 samples from $P(\bm{x}_{1}|\bm{x}^{*}_{3},\bm{x}^{*}_{4})$ by inpainting, where $(\bm{x}^{*}_{3},\bm{x}^{*}_{4})$ are the ESP and pharmacophore interaction profiles of the ligand. We use the ShEPhERD model that was trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq. In a similar manner to the natural products, we used $n_{3}=75$ and specified $n_{4}$ based on the target $\bm{x}_{4}^{*}$. We swept over 25 values of $n_{1}$ in the range $[n_{1}^{*}-12,n_{1}^{*}+12]$ where $n_{1}^{*}$ is the total number of atoms in the PDB ligand. We generate 20 structures for each choice of $n_{1}$. For valid samples, we extract SMILES strings. Different from the other experiments, we do not use xTB for relaxation as we simply dock the valid molecules from a SMILES representation.

To evaluate valid samples, we implement an adapted version of the Vina_smiles class from the TDC oracles. For each SMILES, we 1) embed a conformer with RDKit ETKDG with a random seed of 123456789; 2) optimize with MMFF94 for 200 steps; 3) prepare the conformer for docking with Meeko (forlilab, ); 4) dock the molecule with Autodock Vina v1.2.5 (Trott & Olson, ; Eberhardt et al., ) using a random seed of 987654321 and exhaustiveness of 32; and 5) obtain a Vina score in kcal/mol, which we use as a (poor) proxy for binding affinity.

Figure: Figure 6: (Left) The true binding pose of the experimental ligand extracted from the PDB (pink) overlaid on the top-scoring re-docked pose (light blue). (Middle) Ensemble of the top-10 docked poses for the experimental ligand, ranked by Vina docking score. (Right) Top-1 Vina-scoring pose for each of the three select ShEPhERD samples highlighted in Fig. 4 (green-grey), overlaid on the top-scoring re-docked pose for the experimental ligand (light blue). The top and bottom rows show results for 5mo4’s and 7l11’s experimental co-crystal ligands, respectively. Although the top-1 docked pose for the ShEPhERD-generated ligands may not perfectly align with the top scoring pose of the experimental ligand, the ensembles of their docked poses do closely match, indicating that the experimental ligands and the ShEPhERD-diversified analogues explore common binding modes. We note that true experimental binding poses often do not match the top-1 docked poses simulated from docking programs like AutoDock Vina, and hence we show the ensembles of docked poses to account for this uncertainty.
Refer to caption: extracted/5947291/figures/appendix_figures/docking_si.jpg

Fig. 6 showcases examples for the 5mo4 and 7l11 docking targets. For each target, we show the PDB ligand in its crystal structure pose, the top-scoring redocked pose, and its ensemble of top-10 scoring poses. We then visualize the top scoring pose of three top-scoring ShEPhERD-generated ligands, overlaid on the best docking pose of the PDB ligand. Although the top-scoring docked pose of the ShEPhERD-generated ligands may not align perfectly with the top-scoring pose of the PDB ligand, it is evident that the generated ligands explore poses that engage in the same binding interactions as PDB ligands when we compare their ensembles of docked poses. Moreover, as Fig. 4 illustrates, in each case we can identify a docked pose for each of the ShEPhERD-generated ligands that closely aligns with the top pose of the PDB ligand.

Figure: Figure 7: The distributions of the top-10 scoring ShEPhERD-generated samples when conditioning on the PDB ligand’s lowest energy conformer (teal) or bound pose (light blue), for all protein targets. We compare to the distributions for the top-10 scoring compounds from the 10K docking screen (red). The Vina score of the redocked PDB ligand is denoted by the dotted line. Note that in each case, only 500 molecules were generated by ShEPhERD – in contrast to the 10K sampled in the virtual screen. Furthermore, we only docked the ${<}500$ samples for which we could extract a valid SMILES string from the generated 3D conformation.
Refer to caption: x1.png

Fig. 7 shows the distributions of the top-10 Vina docking scores for ShEPhERD-generated ligands compared to the top-10 docking scores obtained by virtually screening 10K molecules from ShEPhERD-MOSES-aq. Despite having no knowledge of the protein pocket and only generating ${<}500$ (valid) ligands per target, ShEPhERD enriches the top-1 Vina score relative to the redocked PDB ligand for all 7 targets, 3/7 targets relative to the docking screen when conditioning on the lowest energy conformer, and 5/7 targets relative to the docking screen when conditioning on the PDB ligand’s bound pose. Note that conditioning on the PDB ligand’s bound pose indirectly leaks information about the protein pocket. We emphasize that for the 2/7 cases (3ny8 and 4rlu) where the docking screen outperforms (in terms of top-1 Vina score) the ShEPhERD-generated ligands when conditioning on the PDB ligand’s bound pose, the Vina docking score of the PDB ligand is relatively low (${<}9$ kcal/mol). Since we use ShEPhERD to diversify the PDB ligand while preserving its 3D interactions – not optimize docking score – we do not necessarily expect ShEPhERD to find higher-scoring ligands compared to an unbiased docking screen which can explore larger regions of chemical space.

Figure: Figure 8: Docking results for ShEPhERD-generated molecules when conditioning on the best scoring ligand from the 10K docking screen. We show results when conditioning on the lowest energy conformer (teal) or top docking pose (light blue) of the top-scoring screened ligand. We also show the original 10K docking screen (red), the scores of the redocked PDB ligands (blue dotted line), and the scores of the conditioning molecule (orange dotted line). The top plot shows the full Vina score distributions, while the bottom plot only shows the distributions for the top-10 scoring compounds.
Refer to caption: extracted/5947291/figures/appendix_figures/docking_cond_on_screen.jpg

To investigate the potential of using ShEPhERD to optimize docking score, we also condition ShEPhERD on the lowest energy conformer and docked pose of the best scoring compound from the 10K docking screen. Results are shown in Fig. 8. With this strategy, ShEPhERD generates new ligands with higher top-1 Vina scores for 5/7 targets relative to the best scoring compounds from the 10K docking screen (which served as the conditioning information to ShEPhERD). Notably, ShEPhERD drastically improves Vina scores on the 3ny8 target compared to when conditioning on the experimental PDB ligand, shifting the mean of the top-10 scores from $-9.96$ to $-11.70$ kcal/mol (Fig. 7). This PDB ligand scores poorly (the worst among all PDB ligands). Another interesting phenomenon is that neither the 10K screened molecules nor the ShEPhERD-generated molecules (conditioned on the best from the screen) manage to surpass the docking score of 7l11’s experimental PDB ligand. However, ShEPhERD does improve this docking score when directly conditioning on the 7l11 PDB ligand (Fig. 7), highlighting ShEPhERD’s flexibility in structure-based drug design.

#### A.2.5 Bioisosteric fragment merging

We demonstrate ShEPhERD’s native capacity for bioisosteric fragment merging on the antiviral target EV-D68 3C protease. Starting from the 25 fragments hits in the catalytic pocket identified by Wills et al. (), we downselected 13 fragments that still span the P1, P2, and P3 regions of the active site through visual inspection using Fragalysis (Diamond Light Source, ) and ChimeraX (Meng et al., ). Specifically, we removed fragments that did not have common interactions, that significantly deviated from the fragment ensemble’s aggregate shape, or added no extra information. The final list of fragments were: x0147_0A, x1071_0A, x1083_0A, x1084_0A, x1140_0A, x1163_0A, x1498_0A, x1537_0A, x1594_0A, x1919_0A, x2021_0A, x2099_0A, and x2135_0A.

For each fragment, we added hydrogens using RDKit’s AddHs function with addCoords=True. We then generated an aggregate shape representation by sampling from a surface with a probe radius of 0.6. We compute partial charges with an xTB single point calculation in implicit water for each fragment separately. Then, we compute the ESP at each of the surface points for each fragment, and average the ESP contributions at each point. See App. A.5 for more details.

Extraction of the pharmacophores was a manual task based on selecting those that interacted with the pocket (as detected by Fragalysis and ChimeraX). After manual downselection, pharmacophores of the same type underwent clustering if there were at least two within a threshold distance. These threshold distances were $1.0$ Å for HBAs, $2.0$ Å for HBDs, $1.5$ Å for aromatic groups, and $2.0$ Å for hydrophobes. These were somewhat arbitrary cutoffs that yielded reasonable pharmacophore hypotheses upon visual inspection. Table 1 lists the final 27 selected pharmacophores.

**Table 1: The 27 pharmacophores used to condition ShEPhERD for the bioisosteric fragment merging demonstration. Only hydrogen bond acceptors (HBA), hydrogen bond donors (HBD), aromatic groups, and hydrophobes are contained in these 27 pharmacophores.**
| Type | Coordinates ($\mathbf{P}$) | Vector ($\mathbf{V}$) |
| --- | --- | --- |
| HBA | (-7.5155, -3.7035, -9.0635) | (0.1900, 0.7041, -0.6842) |
|  | (-6.7175, -4.8185, -8.8258) | (0.8205, -0.2287, -0.5240) |
|  | (-2.3550, -7.3310, -0.0453) | (0.7719, 0.4600, -0.4388) |
|  | (-6.5197, -4.3363, -4.1827) | (0.5106, 0.7063, -0.4905) |
|  | (-6.9450, -6.9160, -5.0510) | (-0.1783, -0.4378, -0.8812) |
|  | (-3.2780, -7.2170, 2.2860) | (-0.7287, 0.6270, -0.2755) |
|  | (-5.7270, -2.3950, -0.6220) | (0.4491, -0.6321, 0.6316) |
|  | (-1.9580, -8.8820, 3.3910) | (0.7956, -0.5564, -0.2396) |
|  | (-5.0980, -6.4300, -1.5740) | (0.4851, 0.6853, -0.5432) |
|  | (-5.8830, -1.1980, 0.0490) | (0.4090, -0.3130, 0.8572) |
|  | (-5.6690, 0.0200, -2.0780) | (0.5578, 0.5389, -0.6312) |
| HBD | (-9.4350, -2.6960, -8.2680) | (0.1103, 0.6563, -0.7464) |
|  | (-3.2780, -7.2170, 2.2860) | (0.7506, 0.2245, 0.6214) |
|  | (-1.9580, -8.8820, 3.3910) | (-0.2298, 0.1607, 0.9599) |
|  | (-6.9480, -2.0740, -1.9750) | (-0.3293, 0.3106, -0.8917) |
| Aromatic | (-7.6710, -4.5475, -8.2231) | (-0.5639, -0.4342, -0.7025) |
|  | (-11.6032, -2.0268, -6.4003) | (0.3775, 0.9134, -0.1523) |
|  | (-4.6030, -8.0915, -1.0500) | (0.4302, -0.8132, -0.3919) |
|  | (-9.5312, -4.0857, -7.3044) | (-0.5073, -0.3479, -0.7884) |
| Hydrophobe | (-7.6710, -4.5475, -8.2231) | (0.0, 0.0, 0.0) |
|  | (-11.6032, -2.0268, -6.4003) | (0.0, 0.0, 0.0) |
|  | (-4.6030, -8.0915, -1.0500) | (0.0, 0.0, 0.0) |
|  | (-9.5312, -4.0857, -7.3044) | (0.0, 0.0, 0.0) |
|  | (-5.6776, -0.4634, -1.0174) | (0.0, 0.0, 0.0) |
|  | (-6.9299, -5.4087, -1.1730) | (0.0, 0.0, 0.0) |
|  | (-6.7990, -4.5530, -3.9980) | (0.0, 0.0, 0.0) |
|  | (-7.9010, 0.1820, -0.7890) | (0.0, 0.0, 0.0) |

Using the extracted ESP ($\bm{x}_{3}$) and pharmacophore ($\bm{x}_{4}$) interaction profiles, we used the version of ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq to generate 1000 structures from $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ via inpainting. We used $n_{3}=75$, specified $n_{4}=27$, and swept over 40 values of $n_{1}$ in the range $[50,89]$ at 25 samples each. For the valid samples, we relaxed the generated 3D molecular conformation with xTB, and realigned the relaxed structure to the natural product by optimizing our ESP similarity scoring function after resampling $\bm{x}_{3}$ for both structures at $n_{3}=400$. We used ESP-based alignment instead of pharmacophore-based alignment as ESP surfaces better capture 3D shapes, and because ESP-based alignments are more robust. Given the ESP-aligned relaxed structure, we compute the ESP and pharmacophore similarity scores. Fig. 9 visualizes four ShEPhERD-generated samples with neutral charge and SA$\leq 4.0$ that were in the top-10 by combined ESP and pharmacophore similarity scores.

Figure: Figure 9: Examples of four top-scoring ShEPhERD-generated ligands obtained via bioisosteric fragment merging. Shown are the xTB-relaxed structures of the generated 3D molecules after being realigned to the fragments’ ESP and pharmacophore interaction profiles. Alignment is performed by optimizing ESP similarity. The ligands are overlaid on the target interaction profiles. The target ESP surface has been upsampled for visualization.
Refer to caption: extracted/5947291/figures/appendix_figures/fragment_merging_samples.jpg

### A.3 Defining and analyzing validity of generated 3D molecules

Table: Table 2: Global validity reported as percentages for 1000 unconditional and 2000 conditional samples from various ShEPhERD models trained on ShEPhERD-GDB17. Here, “Validity” is the percentage of valid molecules (see definition below); “Validity xTB” is the percentage of valid molecules after xTB optimization; “Graph consis.” is the percentage of samples where the atom connectivity remains the same given that the molecules are valid both pre- and post- xTB-optimization; and “Valid neutral” is the percentage of valid samples that have molecular charge $=0$.

Defining validity. We define validity as the fraction of molecules that can converted from the atomic point cloud $(\bm{a},\bm{C})$ (i.e., element types and coordinates)(^3^33Although ShEPhERD generates formal charges $\bm{f}$ and the bonds $\bm{B}$, we ignore these outputs when extracting a molecule from $\bm{x}_{1}$. We argue that the validity of a generated 3D molecular structure should depend only on the generated 3D structure, not on a jointly generated 2D graph. Hence, we use the diffusion processes over $\bm{f}$ and $\bm{B}$ primarily as a way to help regularize the diffusion processes over $\bm{a}$ and $\bm{C}$. Moreover, by defining validity only with respect to the element types and coordinates, our validity checks may be applied to other 3D molecular generative models which do not generate the 2D molecular graph.) to an RDKit Mol object via RDKit’s internal xyz2mol functionality that satisfy valencies as well as other constraints enumerated below.
Empirically, we find that this method of extracting molecules from xyz-formatted structures is simple, convenient, and quite reliable when (1) explicit hydrogens are included in the 3D structure, (2) the 3D molecular geometry is of high quality, and (3) the overall charge of the molecule is known.

In particular, RDKit is used to (1) generate a Mol object from $(\bm{a},\bm{C})$ formatted as an .xyz file; (2) determine bonds via a rule-based algorithm based on interatomic distances; (3) sanitize the molecule; (4) ensure that there are no radical electrons (which aren’t included in our datasets); (5) check that structure is not fragmented; and (6) check that there are no more than six atomic formal charges that are non-zero (which usually indicates that a valid resonance structure of an aromatic ring could not be determined). In an effort to handle charged molecules, we run this procedure for molecular charges in this order: [0, +1, -1, +2, -2] and break the loop if extraction of a molecule is successful. Finally, an xTB single point calculation is run to obtain partial charges. The molecular charge is specifically used to determine the bonds, to sanitize the molecule, and to run the single point calculation. The sampled $\bm{x}_{1}$ is determined invalid if any of the aforementioned steps fails.

We point out that in some cases, generated molecules may fail our stringent validity checks due to having imperfect molecular geometries (e.g., a bond length is slightly out of distribution). To consider these cases and to measure a notion of distance from the training distribution, we perform an xTB local geometry optimization with the generated $(\bm{a},\bm{C})$ (formatted as an .xyz file) and repeat the validity checks. If the natively generated structure was determined invalid pre-optimization, the xTB optimization is performed with a neutral overall charge; otherwise the same charge was used as was determined upon initial molecule extraction with RDKit.

If both the pre- and post-optimized samples are valid, we check if the graphs are consistent by matching their SMILES strings. If consistent, we also compute the strain energy and heavy-atom RMSD between the relaxed and unrelaxed conformations.

Analyzing validity. Metrics of validity for unconditional and conditional across all samples from ShEPhERD-GDB17 are presented in Table 2. There is a trend of relatively high validity ($\geq 96\%$) for the joint generation of $\bm{x}_{1}$ and $\bm{x}_{2}$, slightly lower validity ($91$-$93\%$) for the joint generation of $\bm{x}_{1}$ and $\bm{x}_{3}$, and a significant decrease for the joint generation of $\bm{x}_{1}$ and $\bm{x}_{4}$. We hypothesize that the lower validity of the pharmacophore-based models is because of the discrete nature of the pharmacophore interaction profile relative to the surface or ESP-attributed surface. Notably, validity improves in the pharmacophore-conditioned sampling scheme compared to the unconditional setting. Constraining the sampling space by specifying the pharmacophores during inpainting may be an easier task as it forces the model to form certain groups/substructures. Specifying the pharmacophores based on a target profile also avoids relying on ShEPhERD’s jointly-generated $\bm{x}_{4}$ (which may have errors) to guide the denoising of $\bm{x}_{1}$.
One way to improve joint pharmacophore generation may be to introduce a noise schedule that denoises the pharmacophores earlier than $(\bm{a},\bm{C})$, as a means to implicitly condition the denoising of the 3D molecules on less noisy pharmacophore representations.

Graph consistency is high for most models which implies that (1) the inter-atomic geometry of valid structures are within the bounds for the the bond-determining rules, and (2) the generated conformations are precise enough that relaxation with xTB is stable. The fraction of valid molecules with neutral charge ranges from $84.1\%-93.2\%$ for unconditional samples and $72.6\%-84.5\%$ for the conditional samples. This is somewhat surprising because the training distribution only contains neutral species. We attribute this behavior to the presence of zwitterions (neutral molecules that contain balanced atoms/groups of opposite charge) in the training sets; the model has learned to generated charged substructures, but hasn’t been trained long enough to learn to neutralize the overall molecule every time. Encouragingly, there is a higher rate of generating neutral species for the $P(\bm{x}_{1},\bm{x}_{3})$ and $P(\bm{x}_{1},\bm{x}_{4})$ models relative to the $P(\bm{x}_{1},\bm{x}_{2})$ model, which perhaps implies that using representations that implicitly encode information about the charge improves the model’s ability to learn the correct distribution of molecular charge. In contrast, the conditional $P(\bm{x}_{1}|\bm{x}_{3})$ and $P(\bm{x}_{1}|\bm{x}_{4})$ models perform worse than $P(\bm{x}_{1}|\bm{x}_{2})$ model in terms of generating neutral molecules. This may be because the model learns to associate certain ESP or pharmacophoric criteria with the presence of charged groups (e.g., carboxylate anions yield negative ESP surfaces), and hence introduces charged groups more often while failing to neutralize the overall molecule.

Fig. 10 details how validity is impacted by the selection of the number of atoms ($n_{1}$) for unconditionally generated samples (top row) and conditionally generated samples (bottom row) from models trained on ShEPhERD-GDB17. In both unconditional and conditional settings, versions of ShEPhERD that jointly model shape or electrostatics attain high validity rates when the number of atoms is in-distribution of the training set (Fig. 21). On the other hand, the pharmacophore models generally generate more invalid structures as the number of atoms increase. We hypothesize that this is related to the difficulty of modeling discrete pharmacophores and satisfying all (jointly-generated or conditional) pharmacophore criteria as the generated molecules become larger and hence more complex. We emphasize that in all cases, ShEPhERD shows the capacity to extrapolate into unseen chemical space (e.g., $n_{1}$ exceeding the bounds of the training datasets), albeit at lower validity rates.

Table 3 and Figures 11-12 show the strain energy and heavy-atom RMSD of ShEPhERD-generated molecular structures compared to their xTB-relaxed geometries. Note that RMSDs are computed after realigning (by minimizing RMSD) the relaxed geometry to the unrelaxed geometry. Mean strain energies remain relatively low for the model trained on ShEPhERD-GDB17 (mean strain energies ${<}10$ kcal/mol). Indeed, these strain energies are typically lower than the strain energies of molecules with MMFF94-optimized geometries that have their heavy-atom RMSD and strain energies measured after relaxing the MMFF94 geometries with xTB. However, the high standard deviations in the strain energy distributions indicate that there are still many generated molecules that are highly strained (up to nearly 1000 kcal/mol). Nevertheless, almost all of the RMSD’s are ${<}1$Å which implies that the generated molecules are close to the true xTB geometry. Indeed, the unconditionally-generated samples for models trained on ShPhERD-GDB17 have average heavy-atom RMSDs below $0.08$ Å. We note that slightly distorted atomic positions can vastly increase strain energy, which makes strain energy an unforgiving (but still useful) measure of absolute geometric quality. We notice a similar trend to validity where samples from $P(\bm{x}_{1},\bm{x}_{2})$ model perform the best in terms of strain energy while samples from $P(\bm{x}_{1},\bm{x}_{4})$ model perform the worst, though the RMSD’s are quite low across all the unconditional models.

The models trained on ShEPhERD-MOSES-aq perform worse in terms of both the strain energy and RMSD metrics (Table 3), though the trends with respect to $n_{1}$ are similar (Fig. 12). We note that the molecules in ShEPhERD-MOSES-aq are generally larger than those in ShEPhERD-GDB17. For instance, the largest molecules in MOSES have 27 non-hydrogen atoms, whereas the largest molecules in GDB17 have 17 non-hydrogen atoms. ShEPhERD-MOSES-aq also contains fewer training examples than ShEPhERD-GDB17; we expect the performance of models trained on ShEPhERD-MOSES-aq will improve with further training on larger drug-like datasets.

Figure: Figure 10: Validity as a function of the number of generated atoms ($n_{1}$) for 3D molecules $\bm{x}_{1}$ generated by various ShEPhERD models trained on ShEPhERD-GDB17 in (Top Row) unconditional or (Bottom Row) conditional settings. Blue circles and orange triangles indicate the validity of generated molecular structures before and after xTB optimization, respectively.
Refer to caption: x2.png

Table: Table 3: The overall mean and standard deviation of strain energies and RMSDs for unconditional samples from versions of ShEPhERD trained on ShEPhERD-GDB17 and ShEPhERD-MOSES-aq (Fig. 28). We compute the metrics from the valid samples amongst 1000 total samples from each model. The column “MMFF94” is a baseline that subsamples 1000 molecules from each respective dataset, generates an MMFF94-level conformer for each molecule, and measures the strain energy and heavy-atom RMSD upon xTB relaxation. Molecules are relaxed with xTB in the gas phase for ShEPhERD-GDB17, and in implicit water for ShEPhERD-MOSES-aq.

Figure: Figure 11: Strain energy and heavy-atom RMSD for samples unconditionally generated by various models trained ShEPhERD-GDB17. The box plots are binned by the number of atoms. Though the strain energies can be high, the RMSDs are quite low. Summary statistics are reported in Table 3.
Refer to caption: x3.png

Figure: Figure 12: Strain energy and heavy-atom RMSD for valid samples (amongst 1000 total samples) unconditionally generated by ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq. The box plots are binned by the number of atoms. Though the strain energies can be high, the RMSDs are quite low. Summary statistics are reported in Table 3.
Refer to caption: x4.png

Figure: Figure 13: Validity as a function of the number of generated atoms for ligands generated with ShEPhERD for the natural product (NP) ligand hopping, bioactive hit diversification, and bioisosteric fragment merging (FM) experiments. For the NP analog generation and FM experiments, we report validity both before and after relaxation of $\bm{x}_{1}$ in implicit water with xTB. For the bioactive hit diversification experiments, we only check the ability to extract SMILES strings from the generated $\bm{x}_{1}$ pre-xTB-optimization. There are no clear patterns in the validity of the molecules generated by ShEPhERD in these experiments, but we attribute the overall lower validity rates to the difficulties of each task and the need for out-of-distribution sampling (Fig. 20), particularly for the NP and FM experiments. Validity rates for each experiment are summarized in Table 4.
Refer to caption: x5.png

Table: Table 4: Global validity metrics for molecules generated by ShEPhERD for the natural product ligand hopping, bioactive hit diversification, and bioisosteric fragment merging experiments.

### A.4 Unconditional and conditional generation metrics

Novelty is defined as the fraction of valid molecules that are not contained in the training set. Uniqueness is defined as the fraction of valid molecules that are only generated once across the samples from a particular model in a given experiment. Novelty and uniqueness metrics for all generated samples (in both unconditional- and conditional-generation experiments) are reported in Table 5. ShEPhERD attains nearly 100% novelty and uniqueness for all cases.

Fig. 14 also reports specific property distributions (SA Score, logP, QED, Fsp3) of unconditionally-generated molecules from versions of ShEPhERD trained on ShEPhERD-GDB17 or ShEPhERD-MOSES-aq. Property distributions for generated molecules are compared against the distributions from a random sampling of molecules from the corresponding datasets.

**Table 5: Novelty and uniqueness metrics across different versions of ShEPhERD trained on either ShEPhERD-GDB17 or ShEPhERD-MOSES-aq. The first six rows report metrics for ShEPhERD models when trained on ShEPhERD-GDB17 and applied to either unconditional ($P(\bm{x}_{1},\bm{x}_{i})$) or conditional ($P(\bm{x}_{1}|\bm{x}_{i})$) generation. The remaining rows report metrics for ShEPhERD when trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq. $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ indicates this model when applied to unconditionally generate 1000 samples, whereas $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ indicates this model when applied to conditional generation in the natural product (NP), bioactive hit diversification, and bioisosteric fragment merging (FM) experiments.**
|  | Samples | Novelty | Uniqueness |
| --- | --- | --- | --- |
|  |  | (%) | (%) |
| ShEPhERD-GDB17 | $P(\bm{x}_{1},\bm{x}_{2})$ | 96.4 | 99.5 |
|  | $P(\bm{x}_{1},\bm{x}_{3})$ | 96.7 | 99.6 |
|  | $P(\bm{x}_{1},\bm{x}_{4})$ | 96.0 | 99.9 |
|  | $P(\bm{x}_{1}|\bm{x}_{2})$ | 99.4 | 97.9 |
|  | $P(\bm{x}_{1}|\bm{x}_{3})$ | 99.3 | 97.6 |
|  | $P(\bm{x}_{1}|\bm{x}_{4})$ | 99.0 | 96.0 |
| ShEPhERD-MOSES-aq | $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| NP 1 | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 99.7 |
| NP 2 | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| NP 3 | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 1iep (lowest energy) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 1iep (pose) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 3eml (lowest energy) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 3eml (pose) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 3ny8 (lowest energy) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 3ny8 (pose) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 4rlu (lowest energy) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 99.5 | 100.0 |
| 4rlu (pose) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 99.8 | 100.0 |
| 4unn (lowest energy) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 4unn (pose) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 5mo4 (lowest energy) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 5mo4 (pose) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 7l11 (lowest energy) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| 7l11 (pose) | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |
| Bioisosteric FM | $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ | 100.0 | 100.0 |

Figure: Figure 14: Distributions of molecular properties for unconditionally-generated molecules compared to samples from the respective datasets. SA Score is the synthetic accessibility score, logP is the octanol-water partition coefficient, QED is a measure of drug-likeness, and Fsp3 is the fraction of sp3 carbons in a molecule. (Top) The distributions relevant to ShEPhERD-GDB17. The plots show the property distributions for the valid molecules amongst 1000 ShEPhERD-generated samples from $P(\bm{x}_{1},\bm{x}_{2})$, $P(\bm{x}_{1},\bm{x}_{3})$, and $P(\bm{x}_{1},\bm{x}_{4})$. “Dataset Samples” refers to 1000 randomly sampled molecules from the ShEPhERD-GDB17 dataset. (Bottom) The same distributions relevant to ShEPhERD-MOSES-aq. The plots show the property distributions for the valid molecules amongst 1000 ShEPhERD-generated samples from $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$. “Dataset Samples” refers to 1000 randomly sampled molecules from the ShEPhERD-MOSES-aq dataset.
Refer to caption: x6.png

### A.5 Calculating interaction profiles from 3D molecular structures

#### A.5.1 Shapes/surfaces

We represent the shape of a molecule as a point cloud sampled from the solvent accessible surface with probe radius $r_{\text{probe}}$. Here, we describe how to extract $\bm{x}_{2}=\bm{S}_{2}\in\mathbb{R}^{n_{2}\times 3}$ from $\bm{x}_{1}=(\bm{a},\bm{C})$ with $\bm{C}\in\mathbb{R}^{n_{1}\times 3}$.

First, we sample points on the surfaces of $n_{1}$ unit spheres using Fibonacci sampling. The number of sampled points for each sphere is $25(r_{a[k]}/1.7)^{2}$ where $r_{a[k]}$ is the van der Waals (vdW) radius of the $k$th atom in Å. The points are then scaled to have norm $r_{a[k]}+r_{\text{probe}}$, and the sampled spheres are translated to their respective atomic coordinates in $\bm{C}$. Points that fall within $r_{a[k]}+r_{\text{probe}}$ of any other atom are filtered out. Next, we use Open3D (Zhou et al., ) to generate a TriangleMesh using the ball pivoting algorithm with a ball radius of $1.2$. Finally, we sample a surface point cloud with approximately evenly spaced points using Open3D’s sample_points_poisson_disk method. This final step is stochastic. We used $r_{\text{probe}}=0.6$ for the modeling in this work, but tuning of similarity scoring function parameters (App. A.6) originally used $r_{\text{probe}}=1.2$.

#### A.5.2 Electrostatic potential surfaces

We obtain the electrostatic potential (ESP) surface by computing the Coulombic potential at each point on a surface point cloud. Here, we describe how to extract $\bm{x}_{3}=(\bm{S}_{3},\bm{v})$ from $\bm{x}_{1}=(\bm{a},\bm{C})$.

We obtain $\bm{S}_{3}\in\mathbb{R}^{n_{3}\times 3}$ using the same procedure as for $\bm{S}_{2}$. We then obtain the partial charges $\bm{q}\in\mathbb{R}^{n_{1}}$ from xTB by running a single point calculation on $\bm{x}_{1}$. Through internal testing, we found that using partial charges obtained from the cheaper MMFF94 to be of insufficient quality and detrimental to the performance of both ShEPhERD and the ESP similarity scoring function. Computing partial charges with xTB is also favorable, as single point calculations can be computed for arbitrary atomistic systems that aren’t parameterized by the MMFF94 force field. Given the partial charges of each atom, we can compute the electrostatic potential at each surface point:

$$ $\displaystyle\bm{v}[k]=\frac{1}{4\pi\epsilon_{0}}\sum^{n_{3}}_{j=1}\frac{q[k]} {\|\bm{r}[k]-\bm{r}[j]\|^{2}}$ $$

where $\epsilon_{0}$ is the vacuum permittivity with units $e^{2}(\text{eV}\cdot\text{\AA})^{-1}$.

#### A.5.3 Pharmacophores

Pharmacophores $\bm{x}_{4}=(\bm{p},\bm{P},\bm{V})$ are abstracted representations of chemical substructures that are associated with common biochemical interactions. They can be directional or non-directional.

The directional pharmacophores considered in this work are:

- •
Hydrogen bond acceptor (HBA) – typically an electronegative nitrogen or oxygen with lone pairs. Vector(s) point from the heavy atom to the direction of a lone pair based on sp, sp2, or sp3 geometries.
- •
Hydrogen bond donor (HBD) – typically an electronegative nitrogen or oxygen with an attached hydrogen. Vector(s) point from the heavy atom to the attached hydrogen(s)
- •
Aromatic ring – aromatic rings can form special interactions such as $\pi$-$\pi$ stacking or cation-$\pi$ interactions. We position an aromatic pharmacophore at the centroid of an aromatic ring. Two antiparallel vectors point in opposite directions, orthogonal to the plane of the aromatic ring.
- •
Halogen-carbon bond – a halogen (F, Cl, Br, I) that is connected to a carbon can form halogen bonds that can display HBA/HBD-like properties. The vector points from the halogen and in the anti-parallel direction to the neighboring carbon.

The non-directional pharmacophores are:

- •
Hydrophobe – a sulfur atom, a halogen atom, or a group of atoms that are hydrophobic. We cluster groups of hydrophobic atoms such that distinct hydrophobes are $>2\text{ \AA}$ apart. Note that aromatic rings are also hydrophobes.
- •
Anion – a negatively charged atom or the centroid of a common anionic group.
- •
Cation – a positively charged atom or the centroid of a common cationic group.
- •
Zinc binder – a functional group capable of coordinating with a zinc ion.

We use SMARTS patterns from Pharmer, RDKit, and Pmapper (Koes & Camacho, ; Landrum et al., ; Kutlushina et al., ) to identify pharmacophores in a molecule, with three notable changes: 1) feature definitions of positive/negative ionizable are altered to be strictly anions/cations; 2) all hydrophobes within 2 Å are clustered using Berenger & Tsuda ()’s protocol; and 3) separate pharmacophore identities for halogens were made due to their unique interaction geometry (e.g., sigma holes) that warranted a distinction from HBA/HBD pharmacophores (Lin & MacKerell Jr, ; Politzer et al., ). Most of the pharmacophore extraction utilizes altered code from RDKit (Features.FeatDirUtilsRD and Features.ShowFeats modules) (Landrum et al., ). Some pharmacophores can be associated with two or more vectors (e.g., a carbonyl group has strong HBA qualities in the directions of the two lone pairs), but we average these together for ease of modeling in this work – except for aromatic rings where we do not average but choose one randomly, instead. Zinc binders, anions, and cations are rare in our datasets which significantly decreases the probability of generating these pharmacophores with ShEPhERD. Examples of each type of pharmacophore are shown in Fig. 15.

The placement of pharmacophore positions $\bm{P}$ are either at the location of the identified atom or at the centroid of a larger substructure – depending on the SMARTS feature definition. For example, some Zn binder SMARTS feature definitions from RDKit place different weights on certain atoms (Fig. 15). Each pharmacophore is associated with a vector $\bm{V}$ and can either be a unit vector ($|\bm{V}[k]|=1$) for directional pharmacophores or the zero vector ($|\bm{V}[k]|=0$) for directionless pharmacophores. For directional pharmacophores, a unit vector in the direction of a potential interaction (e.g., in the direction of a hydrogen for an HBD or a lone pair for an HBA) is extracted. It is common for multiple pharmacophores to be located at the same position. For example, aromatic rings and halogen-carbon bonds are also hydrophobes (Fig. 15), and hydroxyl groups are both HBDs and HBAs.

Figure: Figure 15: Examples of each type of pharmacophore. (Top) A molecule is shown that contains common pharmacophores, including hydrophobes (grey), HBA (purple), HBD (light blue), aromatic rings (orange), and halogen-carbon bonds (green). The left 3D structure highlights how multiple vectors are used for HBAs when they have multiple lone pairs, while the right structure shows these vectors averaged into a single direction (not shown here but the same applies for HBDs). Additionally, the clustering scheme from (Berenger & Tsuda, ) is used to reduce the number of hydrophobes. This molecule is an slightly altered example from [https://greglandrum.github.io/rdkit-blog/posts/2023-02-24-using-feature-maps.html](https://greglandrum.github.io/rdkit-blog/posts/2023-02-24-using-feature-maps.html). (Bottom) Examples are provided for pharmacophores that are rare in our datasets: cations, anions, and zinc binders. The pyridinium structure features a cation (large dark blue sphere) on the nitrogen atom. The hydrogen sulfate example shows an anion (large red sphere), with the SMARTS pattern covering the entire group, placing the feature at the centroid. The hydroxamate example contains both a cation on the oxygen and a zinc binder (pink sphere). Although the entire hydroxamate group is classified as a Zn binder according to the SMARTS pattern, the pharmacophore weights are specifically assigned to the two oxygens capable of coordinating a zinc ion, which places the pharmacophore between the oxygens rather than at the group centroid.
Refer to caption: extracted/5947291/figures/appendix_figures/pharmacophores.jpg

### A.6 Parameterization of 3D similarity scoring functions

Recall from Sec. 3 that for two point clouds $\bm{Q}_{A}$ and $\bm{Q}_{B}$ we model each point $\bm{r}_{k}$ as an isotropic Gaussian in $\mathbb{R}^{3}$ (Grant & Pickup, ; Grant et al., ). We compute the overlap between them $O_{A,B}$ with the first order Gaussian overlap. The Tanimoto similarity function $\text{sim}^{*}(\bm{Q}_{A},\bm{Q}_{B})\in[0,1]$ is used to define a similarity of two point clouds based on their overlaps:

$$ $\displaystyle O_{A,B}$ $\displaystyle=\sum_{a\in\bm{Q}_{A}}\sum_{b\in\bm{Q}_{B}}w_{a,b}\left(\frac{\pi }{2\alpha}\right)^{\frac{3}{2}}\exp{\left(-\frac{\alpha}{2}\|\bm{r}_{a}-\bm{r} _{b}\|^{2}\right)}$ $\displaystyle\text{sim}^{*}(\bm{Q}_{A},\bm{Q}_{B})$ $\displaystyle=\frac{O_{A,B}}{O_{A,A}+O_{B,B}-O_{A,B}}$ $$

Here, $\alpha$ parameterizes the Gaussian width and $w_{a,b}$ is a weighting factor that is tuned for each representation. Note that $n_{Q_{A}}$ need not equal $n_{Q_{B}}$, generally.

Since 3D similarities are sensitive to SE(3) transformations of $\bm{Q}_{A}$ with respect to $\bm{Q}_{B}$, their optimal alignment can be defined as:

$$ $\text{sim}(\bm{Q}_{A},\bm{Q}_{B})=\max_{\bm{R,t}}{\text{sim}^{*}(\bm{RQ}_{A}^{ T}+\bm{t},\bm{Q}_{B}})$ $$

where $\bm{R}\in SO(3)$ and $\bm{t}\in T(3)$. The optimization landscape is non-convex and noisy. For global optimizations (e.g., aligning random 3D molecules from the dataset to a target natural product), we optimize the alignment of $\bm{Q}_{A}$ with respect to $\bm{Q}_{B}$ with $N$ different initializations of ($\bm{R}$, $\bm{t}$). We found that $N=50$ is sufficient for convergence. For local optimizations (e.g., aligning a 3D molecule generated by ShEPhERD to a target natural product), we use the original alignment (e.g., the pose natively generated by ShEPhERD) as the only initialization and optimize directly via gradient descent.

For initializations of global alignment optimization, $\bm{t}$ is set so that the center of mass (COM) of $\bm{Q}_{A}$ is aligned with the COM of $\bm{Q}_{B}$. One of the initializations uses the identity $\bm{R}=\bm{I}$ (to keep the initial orientation). We also use alignments to four orientations relative to the principal moments of inertia of $\bm{Q}_{B}$ (the positive and negative orientations for the two largest principal components). For all other initializations when $N>5$, we use Fibonacci sampling to obtain evenly spaced out rotations in $SO(3)$.

We use unit quaternions to represent $\bm{R}$ to decrease the number of parameters and increase optimization stability. We achieve optimal alignment through automatic differentiation and use Adam optimizer for a maximum of 200 steps with a learning rate of 0.1. By implementing in PyTorch and Jax, we run batched computations across the $N$ initializations.

#### A.6.1 Shape similarity scoring function

Recall that we follow Adams & Coley () in defining the volumetric shape similarity between two atomic point clouds $\bm{C}_{A}$ and $\bm{C}_{B}$ as $\text{sim}^{*}_{\text{vol}}(\bm{C}_{A},\bm{C}_{B})$ with $w_{a,b}=2.7$ and $\alpha=0.81$. The surface shape similarity between two surfaces $\bm{S}_{A}$ and $\bm{S}_{B}$ is defined as $\text{sim}^{*}_{\text{surf}}(\bm{S}_{A},\bm{S}_{B})$ with $w_{a,b}=1$, and $\alpha=\Psi(n_{2})$.
Here, $\Psi$ is a function fitted to $\text{sim}^{*}_{\text{vol}}$ depending on the choice of $n_{2}$.

We fit $\Psi(n_{2})$ by minimizing RMSE between $\text{sim}^{*}_{\text{vol}}$ and $\text{sim}^{*}_{\text{surf}}$ for $n_{2}\in[50,100,150,200,300,400]$ on 1K randomly sampled molecules from [https://www.kaggle.com/datasets/art3mis/chembl22](https://www.kaggle.com/datasets/art3mis/chembl22) with the number of heavy atoms ranging from 6 to 62. For parameter fitting, we used a probe radius of 1.2 Å when sampling the surfaces. This differs from the 0.6 Å probe radius that we ultimately used when modeling and scoring the similarities for $\bm{x}_{2}$ and $\bm{x}_{3}$ in our experiments, but we do not expect this discrepancy to meaningfully affect our analyses. The conformations were generated with RDKit ETKDG and optimized with MMFF94. Results and correlations from tuned $\alpha$’s are found in Fig. 16. We do not expect, nor want, perfect correlation with the volumetric similarity since we aim to decouple the shape representation from the atomic coordinates. $\Psi(n_{2})$ is defined as an quadratic interpolation function from SciPy used to fit the optimal $\alpha$’s found in Table 6. Note that for surface similarity, we require $n_{S_{A}}=n_{S_{B}}$ due to the tuned $\alpha$. We show the effects of stochasticity in Figure 17.

Figure: Figure 16: The correlations between surface similarity and volumetric similarity after tuning the Gaussian width parameter $\alpha$ to minimize RMSE between the two scoring functions.
Refer to caption: extracted/5947291/figures/appendix_figures/surf_vs_vol_sim_alpha_tuning.jpg

**Table 6: Tuned Gaussian width parameters ($\alpha$) for the surface similarity scoring function. Larger $\alpha$ corresponds to smaller width.**
| Number of surface points ($n_{2}$) | 50 | 100 | 150 | 200 | 300 | 400 |
| --- | --- | --- | --- | --- | --- | --- |
| Gaussian width ($\alpha$) | 0.6011 | 0.8668 | 1.022 | 1.118 | 1.216 | 1.258 |

Figure: Figure 17: The effect of stochastic surface sampling on the surface and ESP scores, as measured by the average self-similarity as a function of the number of atoms $n_{1}$. For each of ${\sim}900$ (xTB-relaxed) generated compounds from ShEPhERD trained on ShEPhERD-GDB17, we repeatedly extract the true shape $\bm{x}_{2}$ or ESP $\bm{x}_{3}$ interaction profile 5 times, score the similarity between pairs of extracted profiles, and compute the average similarity. We call this quantity the shape or ESP “self-similarity” for a given molecule. We plot the average self-similarity as a function of the number of atoms in the molecule $n_{1}$ when extracting $\bm{x}_{2}$ or $\bm{x}_{3}$ with either 75 or 400 surface points. As the surfaces are more densely sampled — either by reducing $n_{1}$ for fixed $n_{2}$ or $n_{3}$, or increasing $n_{2}$ or $n_{3}$ for fixed $n_{1}$ — the average self-similarity approaches 1.0. Recall that we use 75 points when evaluating the self-consistency of ShEPhERD’s jointly generated $(\bm{x}_{1},\bm{x}_{2})$ and $(\bm{x}_{1},\bm{x}_{3})$ samples, but use 400 points when scoring similarities between two molecules in all conditional experiments. Using 400 points for shape and ESP similarity scoring gives a much more precise estimate of interaction similarity.
Refer to caption: x7.png

#### A.6.2 ESP surface similarity scoring function

The similarity between two electrostatic potential (ESP) surfaces $\text{sim}_{\text{ESP}}^{*}(\bm{x}_{3,A},\bm{x}_{3,B})$ is defined with an overlap function:

$$ $O_{A,B}^{\text{ESP}}=\sum_{a\in\bm{Q}_{A}}\sum_{b\in\bm{Q}_{B}}\left(\frac{\pi }{2\alpha}\right)^{\frac{3}{2}}\exp{\left(-\frac{\alpha}{2}\|\bm{r}_{a}-\bm{r} _{b}\|^{2}\right)}\exp{\left(-\frac{\|\bm{v}_{A}[a]-\bm{v}_{B}[b]\|^{2}}{ \lambda}\right)}$ $$

where $\alpha=\Psi(n_{3})$, and $\lambda=\frac{0.3}{(4\pi\epsilon_{0})^{2}}$. $\epsilon_{0}$ is the permittivity of vacuum with units $e^{2}(\text{eV}\cdot\text{\AA})^{-1}$. We chose $\lambda$ through manual inspection of optimal alignments. Fig. 17 shows the effects of stochasticity.
Others have also considered different formulations of ESP scoring functions to characterize bioisosteres (Good et al., ; Cleves et al., ; Bolcato et al., ; Osman & Arabi, ).

#### A.6.3 Pharmacophore similarity scoring function

For pharmacophore similarity, we follow the the vector weighting formulation of PheSA (Wahl, ) and define the overlap similarity between two sets of pharmacophores $\bm{x}_{4,A}$ and $\bm{x}_{4,B}$ as:

$$ $\displaystyle O^{\text{pharm}}_{A,B;m}$ $\displaystyle=\sum_{a\in\bm{Q}_{A,m}}\sum_{b\in\bm{Q}_{B,m}}w_{a,b;m}\left( \frac{\pi}{2\alpha_{m}}\right)^{\frac{3}{2}}\exp{\left(-\frac{\alpha_{m}}{2}\| \bm{r}_{a}-\bm{r}_{b}\|^{2}\right)}$ $\displaystyle\text{sim}_{\text{pharm}}^{*}(\bm{x}_{4,A},\bm{x}_{4,B})$ $\displaystyle=\frac{\sum_{m\in\mathcal{M}}O_{A,B;m}}{\sum_{m\in M}O_{A,A;m}+O_ {B,B;m}-O_{A,B;m}}$ $$

where $\mathcal{M}$ is the set of all pharmacophore types ($|\mathcal{M}|=N_{p}$), $\alpha_{m}=\Omega(m)$ where $\Omega$ maps each pharmacophore type to a Gaussian width using the parameters from Pharao (Taminau et al., ). $w_{a,b;m}$ is defined as:

$$ $w_{a,b;m}=\begin{cases}1&\text{if }m\text{ is non-directional},\\ \frac{\bm{V}[a]_{m}^{\top}\bm{V}[b]_{m}+2}{3}&\text{if }m\text{ is directional }.\end{cases}$ $$

We take the absolute value of $\bm{V}[a]_{m}^{\top}\bm{V}[b]_{m}$ for aromatic groups as we assume their $\pi$ interaction effects are symmetric across their plane. Note that Pharao does not define the halogen-carbon bond pharmacophore so we assign $\alpha=1.0$.

### A.7 Additional details on model design, training protocols, and sampling

Feature scaling. Hoogeboom et al. () and Peng et al. () found that linearly scaling certain one-hot features in 3D molecular DDPMs can improve model performance as measured via sample validity, presumably because increasing/decreasing the magnitude of certain features causes them to become resolved earlier/later in the denoising process. We loosely follow these prior works to scale $\bm{a}$ and $\bm{f}$ by 0.25 and $\bm{B}$ by 1.0. We additionally scale $\bm{v}$ and $\bm{p}$ by 2.0. Finally, we scale the pharmacophore directional vectors $\bm{V}$ by 2.0. However, we did not rigorously tune these scaling factors. We also did not explore the use of different noise schedules for different variables, although this strategy has also been found to be helpful for 3D molecule DDPMs (Peng et al., ).

Noise schedule.
ShEPhERD uses the same noise schedule for all diffusion/denoising processes. We choose the variance schedule $\sigma^{2}_{t}$ to be a weighted combination of the linear schedule used by RFDiffusion (Watson et al., ) and a cosine schedule introduced by Nichol & Dhariwal (). Fig. 18 plots our noise schedule in terms of $\sigma_{t}$, $\alpha_{t}$, $\overline{\sigma}_{t}$, and $\overline{\alpha}_{t}$ for $t\in[1,T]$. We use $T=400$ for all ShEPhERD models.

Figure: Figure 18: Noise schedule used by ShEPhERD for all forward noising processes in terms of $\sigma_{t}$, $\alpha_{t}$, $\overline{\sigma}_{t}$, and $\overline{\alpha}_{t}$ for $t\in[1,400]$.
Refer to caption: x8.png

#### A.7.1 Further details on ShEPhERD ’s denoising modules

“E3NN-style” vs. “EGNN-style” coordinate predictions. As described in section 3.3, we predict the $l{=}1$ coordinate noises $\bm{\hat{\epsilon}}_{\bm{S}_{2}}^{(t)}$, $\bm{\hat{\epsilon}}_{\bm{S}_{3}}^{(t)}$, and $\bm{\hat{\epsilon}_{P}}^{(t)}$ from the corresponding $l{=}1$ node features $\bm{\tilde{z}}^{(t)}_{i}$ directly with equivariant feed forward networks that appear in EquiformerV2 and other E3NN-based architectures. In these “E3NN-style” coordinate predictions, the noise for each node $\bm{\epsilon}_{i}[k]$ if predicted from only the $l{=}1$ feature $\bm{\tilde{z}}_{i}[k]$ of that node. While each $\bm{\tilde{z}}_{i}[k]$ implicitly contains information about the node’s surrounding environment (e.g, neighboring nodes) due to ShEPhERD’s joint module, the actual denoising step is not defined with respect to the positions of the neighboring nodes. This notably differs from Hoogeboom et al. ()’s original formulation of equivariant DDPMs for 3D molecule generation, which used “EGNN-style” (Satorras et al., ) coordinate predictions to denoise the molecule’s coordinates $\bm{C}$. For a single-layer EGNN, the coordinate prediction step takes the form (using our notation):

$$ $\bm{C}[k]^{\prime}=\bm{C}[k]+\sum_{j\neq k}\bigg{(}\frac{\bm{C}[k]-\bm{C}[j]}{ d_{kj}+1}\bigg{)}\phi_{\bm{C}}\bigg{(}\bm{z}[k],\bm{z}[j],d_{kj}\bigg{)}\ \ ; \ \ d_{kj}=\bigg{|}\bigg{|}\bm{C}[k]-\bm{C}[j]\bigg{|}\bigg{|}_{2}$ $$

where $\phi_{\bm{C}}$ is an MLP and $\bm{C}[k]^{\prime}$ can be interpreted as an updated coordinate for the node $k$, which is equivariant to E(3) transformations of $\bm{C}$. With this updated coordinate, the noise $\bm{\hat{\epsilon}}_{i}[k]$ may be obtained by computing $\bm{\hat{\epsilon}}_{i}[k]=\bm{C}[k]^{\prime}-\bm{C}[k]$. In “EGNN-style” coordinate predictions, therefore, the predicted coordinate/noise for node $k$ is defined with respect to the other nodes in the 3D graph.

We hypothesize that E3NN-style coordinate predictions are sufficient for predicting $\bm{\hat{\epsilon}}_{\bm{S}_{2}}^{(t)}$, $\bm{\hat{\epsilon}}_{\bm{S}_{3}}^{(t)}$, and $\bm{\hat{\epsilon}_{P}}^{(t)}$ (e.g., the $l{=}1$ coordinates of $\bm{x}_{2}$, $\bm{x}_{3}$, and $\bm{x}_{4}$), as these noise predictions do not need to be so highly sensitive to the positions of neighboring nodes. However, for the coordinates $\bm{C}$ of the 3D molecule $\bm{x}_{1}$, inter-node geometries (e.g., bond lengths and angles) are quite relevant to the quality of the molecule’s conformation. Hence, we hypothesize that it may be favorable to use EGNN-style coordinate predictions to obtain $\bm{\hat{\epsilon}}_{\bm{C}}$. However, EGNN-style coordinate predictions are less expressive than E3NN-style coordinate predictions, as the difference $\bm{C}[k]^{\prime}-\bm{C}[k]$ is just a linear combination $\sum_{j\neq k}w_{j}(\bm{C}[k]-\bm{C}[j])$ with (learnt) weights $w_{j}$. To denoise $\bm{C}$, we thus choose to use both E3NN- and EGNN-style coordinate predictions. We first predict a set of updated coordinates $\bm{C}^{\prime(t)}$ from $\bm{C}^{(t)}$ using E3NN-style coordinate predictions, and then refine $\bm{C}^{\prime(t)}$ with EGNN-style coordinate predictions to obtain $\bm{C}^{\prime\prime(t)}$. The noise $\bm{\hat{\epsilon}}_{\bm{C}}$ is then computed as $\bm{\hat{\epsilon}}_{\bm{C}}=\bm{C}^{\prime\prime(t)}-\bm{C}^{(t)}$. Algorithm 1 outlines this process within the context of the entire denoising module for $\bm{x}_{1}$.

Forward-noising and denoising the covalent bond graph $\bm{B}$.
Multiple works have found that jointly diffusing the 2D covalent bond graph with the 3D atomic coordinates helps improve sample quality for 3D molecule DDPMs (Peng et al., ; Vignac et al., ). While using discrete diffusion for bond/atom types has been found to further improve sample quality compared to diffusing/denoising continuous one-hot representations of bond/atom types, we choose to use continuous representations of bond/atom types to emphasize overall modeling simplicity, consistency, and extensibility. In this section, we clarify details regarding the forward and reverse processes for $\bm{B}$.

We represent $\bm{B}\in\mathbb{R}^{n_{1}\times n_{1}\times 5}$ as an $n_{1}\times n_{1}$ adjency matrix attributed with one-hot features of the covalent bond order between each pair of atoms in $\bm{x}_{1}$. The bond types include single, double, triple, and aromatic bonds, as well as an extra bond type denoting the absence of a covalent bond. It is important to note that we treat bonds between pairs of atoms as undirected edges by symmetrizing $\bm{B}^{(t)}$ in both the forward noising and reverse denoising processes. Namely, we only forward-noise and denoise the upper triangle of $\bm{B}^{(t)}$, and copy the upper triangle to the lower triangle. We mask-out the diagonal of $\bm{B}^{(t)}$ so that self-loops are not modeled.

To maintain permutation invariance when predicting the noise $\bm{\hat{\epsilon}}_{B}[k,j]$ for the bond between atoms $k$ and $j$ from their $l{=}0$ node representations $\bm{z}_{1}[k]$ and $\bm{z}_{1}[j]$, we use a permutation-invariant MLP:

$$ $\bm{\hat{\epsilon}}_{B}[k,j]=\frac{1}{2}\bigg{(}\bm{b}_{kj}+\bm{b}_{jk}\bigg{)}$ $$

$$ $\bm{b}_{kj}=\texttt{MLP}_{\bm{B}}\bigg{(}\bigg{[}\bm{B}^{(t)}[k,j],f_{RBF}(d_{ kj}),\bm{z}_{1}[k],\bm{z}_{1}[j]\bigg{]}\bigg{)}$ $$

$$ $\bm{b}_{jk}=\texttt{MLP}_{\bm{B}}\bigg{(}\bigg{[}\bm{B}^{(t)}[j,k],f_{RBF}(d_{ jk}),\bm{z}_{1}[j],\bm{z}_{1}[k]\bigg{]}\bigg{)}$ $$

where
$d_{kj}=d_{jk}=||\bm{C}^{(t)}[k]-\bm{C}^{(t)}[j]||$, $f_{RBF}$ is a radial basis function expansion using Gaussian basis functions, and $[.,.]$ denotes concatenation in the feature dimension.

#### A.7.2 Training and sampling procedures

Algorithms 2, 3, 4, 5 outline the forward pass of ShEPhERD’s denoising network and detail how we train and sample from ShEPhERD. Algorithm 1 details the forward pass of ShePhERD’s denoising module for $\bm{x}_{1}$. The denoising modules for $\bm{x}_{2}$, $\bm{x}_{3}$, and $\bm{x}_{4}$ are simpler and described in section 3.

Training details.
We train ShEPhERD with the Adam optimizer using a constant learning rate of 3e-4 and an effective batch size ranging from 40 to 48. We clip gradients that have norm exceeding 5.0. App. A.11 describes overall training times, and App. A.7.3 lists all hyperparameters.

Each training example requires sampling a time step $t\in[1,T]$ (we use $T=400$). Rather than uniformly sampling $t$ within that range, we upsample $t$ in the range $[50,250]$ to accelerate model convergence, since this interval captures the most important part of the denoising trajectory. Specifically, we uniformly sample $t\in[50,250]$ $75\%$ of the time, $t\in[0,50)$ $7.5\%$ of the time, and $t\in(250,400]$ $17.5\%$ of the time.

Figure: Algorithm 1 Denoising Module for $\bm{x}_{1}$

Figure: Algorithm 2 Forward Pass of ShEPhERD’s Denoising Network $\bm{\eta}$

Figure: Algorithm 3 Training Algorithm

Figure: Algorithm 4 Sampling Algorithm for Unconditional Generation

Figure: Algorithm 5 Sampling Algorithm for Conditional Generation with Inpainting

#### A.7.3 Model hyperparameters

Table 7 lists hyperparameters relevant to training ShEPhERD. Table 8 lists the total number of learnable parameters for each version of ShEPhERD analyzed in the main text. Table 9 lists hyperparameters related to ShEPhERD’s denoising network architecture. Note that model hyperparameters were selected early during model development based on the memory limits of our GPUs, and were not specially tuned to optimize model performance. A notable exception is the probe radius for extracting $\bm{x}_{2}$ and $\bm{x}_{3}$, which was manually tuned to $0.6$ Å in order to decouple the shape/surface representation from the 2D molecular graph or exact atomic coordinates while still yielding well-defined shapes. Note that a probe radius of 0.0 corresponds to the van der Waals surface, whereas a very large probe radius would cause all molecules to have uninformative spherical shapes. If one requires the shape/surface representation to be more sensitive to the exact atomic coordinates, one can easily retrain ShEPhERD with a smaller probe radius.

**Table 7: Hyperparameters used for training ShEPhERD**
| Training Parameters |  |
| --- | --- |
| Parameter | Value |
| Effective batch size | 48 |
| Learning rate | 0.0003 |
| Gradient clipping value | 5.0 |
| $T$ | 400 |
| $n_{2}$ | 75 |
| $n_{3}$ | 75 |
| surface probe radius | 0.6 |
| $\bm{a}$ and $\bm{f}$ scaling factors | 0.25 |
| $\bm{B}$ scaling factor | 1.0 |
| $\bm{v}$ scaling factor | 2.0 |
| $\bm{p}$ scaling factor | 2.0 |
| $\bm{V}$ scaling factor | 2.0 |

**Table 8: Total number of learnable parameters for each ShEPhERD model**
| Total Number of Learnable Parameters |  |
| --- | --- |
| Model | Number of Parameters |
| ShEPhERD-GDB17 $P(\bm{x}_{1},\bm{x}_{2})$ | 4,375,054 |
| ShEPhERD-GDB17 $P(\bm{x}_{1},\bm{x}_{3})$ | 4,387,407 |
| ShEPhERD-GDB17 $P(\bm{x}_{1},\bm{x}_{4})$ | 4,407,321 |
| ShEPhERD-MOSES-aq $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ | 6,010,427 |

**Table 9: Hyperparameters for ShEPhERD’s denoising network $\bm{\eta}$**
| Denoising Network Hyperparameters |  |
| --- | --- |
| Parameter | Value |
| Default EquiformerV2 Parameters |  |
| num_node_channels | 64 |
| lmax_list | [1] |
| mmax_list | [1] |
| ffn_hidden_channels | 32 |
| grid_resolution | 16 |
| num_sphere_samples | 128 |
| edge_channels | 128 |
| activation_function | silu |
| norm_type | layer_norm_sh |
| use_sep_s2_act | True |
| use_grid_mlp | True |
| use_gate_act | False |
| use_attn_renorm | True |
| use_s2_act_attn | False |
| Joint Module Parameters |  |
| num_EquiformerV2_layers | 2 |
| attention_channels | 24 |
| num_attention_heads | 2 |
| radius_graph_cutoff | 5.0 |
| RBF_cutoff | 5.0 |
| $\bm{x}_{1}$ Embedding Module Parameters |  |
| num_EquiformerV2_layers | 4 |
| attention_channels | 32 |
| num_attention_heads | 4 |
| ffn_hidden_channels | 64 |
| radius_graph_cutoff | $\infty$ (fully connected) |
| RBF_cutoff | 5.0 |
| $\bm{x}_{2}$, $\bm{x}_{3}$, $\bm{x}_{4}$ Embedding Module Parameters |  |
| num_EquiformerV2_layers | 2 |
| attention_channels | 24 |
| num_attention_heads | 2 |
| ffn_hidden_channels | 32 |
| radius_graph_cutoff | 5.0 |
| RBF_cutoff | 5.0 |
| x3_scalar_RBF_expansion_min | -10.0 |
| x3_scalar_RBF_expansion_max | 10.0 |
| $\bm{x}_{1}$ Denoising Module Parameters |  |
| MLP_hidden_dim | 64 |
| num_MLP_hidden_layers | 2 |
| e3nn_ffn_hidden_channels | 32 |
| egnn_normalize_vectors | True |
| egnn_distance_expansion_dim | 32 |
| $\bm{x}_{2}$, $\bm{x}_{3}$, $\bm{x}_{4}$ Denoising Module Parameters |  |
| MLP_hidden_dim | 64 |
| num_MLP_hidden_layers | 2 |
| e3nn_ffn_hidden_channels | 32 |

### A.8 Symmetry breaking in unconditional generation

We find that during unaltered unconditional generation, ShEPhERD tends to generate spherical molecules regardless of the choice of $n_{1}$ (Fig. 19). We attribute this behavior to ShEPhERD being unable to reliably break spherical symmetry when denoising the Euclidean coordinate components of $\bm{X}^{(t)}$. This phenomenon is particularly evident for models trained on ShEPhERD-GDB17, which contains a number of spherical cage-like molecules in the true data distribution. But, this also occurs for models trained on ShEPhERD-MOSES-aq, which does not contain such molecules. Intuitively, we can understand this phenomenon by considering that we denoise molecular structures and their interaction profiles starting from isotropic Gaussian noise. At early time steps in the denoising trajectory (e.g., $t$ close to $T$), the model learns to collapse the isotropic noise into a near point-mass; this behavior minimizes the L2 denoising losses at early time steps as predicting clean structures from uninformative noisy input states is impossible. Due to ShEPhERD’s strict SE(3)-equivariance, though, the denoising network cannot natively break the spherical symmetry of this point mass; this ultimately leads the network to generate spherically-shaped (yet often still valid) molecules.

Fortunately, we can easily solve this problem by manually adding extra noise to the coordinate components of $\bm{X}^{(t-1)}$ during the most informative time steps. With our choice of noise schedule, the states $\bm{X}^{(t)}$ undergo the most significant evolution in terms of their spatial coordinates from $t=130$ through $t=80$. For $t\in[80,130]$, then, we add extra (isotropic) noise to $\bm{C}^{(t-1)},\bm{S}_{2}^{(t-1)},\bm{S}_{3}^{(t-1)},\bm{P}^{(t-1)}$ so that (as an example):
$\bm{C}^{(t-1)}=\frac{1}{\alpha_{t}}\bm{C}^{(t)}-\frac{\sigma_{t}^{2}}{\alpha_{
t}\overline{\sigma}_{t}}\bm{\hat{\epsilon}}^{(t)}+\frac{\sigma_{t}\overline{
\sigma}_{t-1}}{\overline{\sigma}_{t}}\bm{\epsilon}^{\prime}+\xi\bm{\epsilon}_{
\text{symmetry-breaking}}$
where $\bm{\epsilon}_{\text{symmetry-breaking}}\sim N(\bm{0},\bm{1})$ and $\xi$ scales the noise (although we just use $\xi=1$). For $\bm{C}^{(t-1)}$, we still remove the COM from $\bm{\epsilon}_{\text{symmetry-breaking}}$.

Naively adding extra noise to the coordinates of the denoised states could quickly cause the model to go out-of-distribution, leading to low-quality or invalid samples. To remedy this, we use harmonization/resampling (Lugmayr et al., ) to re-sample the state $\bm{X}^{(t=80+\Delta t)}$ given the symmetry-broken $\bm{X}^{(t=80)}$ (we use $\Delta t=20$). Empirically, we find this adequate to maintain sample quality, as evidenced by the high validity rates and low RMSDs upon xTB-relaxation for unconditional samples from ShEPhERD when trained on ShEPhERD-GDB17. Nevertheless, we admit there may be less intrusive methods of breaking symmetry; e.g., adding non-isotropic noise with smaller $\xi$. However, we find our symmetry-breaking procedure to be simple and effective while only marginally increasing sampling time.

Finally, we emphasize that this symmetry breaking procedure is only performed during unconditional generation. When inpainting during conditional generation, symmetry is naturally broken by the conditioning information, so we do not add any extra noise (although doing so could increase the diversity of conditionally-sampled structures).

Figure: Figure 19: Without symmetry breaking during unconditional generation, ShEPhERD tends to generate spherical molecules, particularly for models trained on ShEPhERD-GDB17. Shown are unconditional samples from ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{2})$ on ShEPhERD-GDB17, generated either with or without symmetry breaking.
Refer to caption: extracted/5947291/figures/appendix_figures/symmetry_breaking.jpg

### A.9 Characterizing the need for out-of-distribution performance

Multiple of our experiments, particularly those related to ligand-based drug design, require ShEPhERD to generate larger molecules with more pharmacophores than the molecules that ShEPhERD was trained on. For instance, the largest molecules in ShEPhERD-MOSES-aq contain 27 heavy atoms, whereas the best-scoring ligands generated in our natural product ligand hopping experiments contain more than that, presumably because these larger molecules fit the shapes of the very large natural products better than molecules with fewer heavy atoms. Moreover, in our fragment merging experiment, we conditioned ShEPhERD on a pharmacophore profile containing 27 pharmacophores, which is many more than those contained by the molecules in ShEPhERD-MOSES-aq. To quantitatively illustrate the need for out-of-distribution performance and ShEPhERD’s extrapolation ability, Fig. 20 shows the distributions of $n_{1}$ vs. $n_{4}$ for molecules from ShEPhERD-MOSES-aq compared against the molecules conditionally generated by ShEPhERD across our bioisosteric drug design experiments. Also included are the top-scoring ligands for each experiment, to highlight that ShEPhERD finds high-scoring compounds that are both in-distribution and meaningfully out-of-distribution. Fig. 21 shows similar distributions for our unconditional and conditional generation experiments on the ShPhERD-GDB17 dataset.

Figure: Figure 20: Distributions of the number of pharmacophores relative to the number of atoms (including hydrogens) for conditional samples from ShEPhERD when trained on ShEPhERD-MOSES-aq. For comparison, the kernel density plot of a randomly sampled subset of $\sim$320k molecules from the ShEPhERD-MOSES-aq dataset is shown in blue on both figures. (Left) Conditional ShEPhERD-generated samples from $P(\bm{x}_{1}|\bm{x}_{3},\bm{x}_{4})$ for our natural product (NP) ligand hopping, bioactive hit diversification, and bioisosteric fragment merging experiments. Note that the ordering of the NPs correspond to the order found in Fig. 4. Because we sweep over $n_{1}$ while specifying $n_{4}$ based on the target $\bm{x}_{4}$, the ShEPhERD-generated samples appear as horizontal lines. (Right) “Top-10” scoring samples generated by ShEPhERD for each experiment. The “Top-10” samples for the NP and fragment merging experiments are defined as the valid samples that have the highest combined ESP & pharmacophore score post xTB-relaxation and ESP-based realignment. The “Top-10” samples for the bioactive hit diversification experiments are defined by their Vina scores. We include the top samples when conditioning on both the lowest energy conformer (LE) and the bound pose of each PDB ligand. Top-scoring samples are found both in- and out-of-distribution in terms of ($n_{1}$, $n_{4}$).
Refer to caption: x9.png

Figure: Figure 21: The distributions of the number of pharmacophores relative to the number of atoms (including hydrogens) for molecules in the ShEPhERD-GDB17 dataset and molecules sampled by ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{4})$. We include the distributions for $\sim$1M randomly sampled molecules from the ShEPhERD-GDB17 dataset, unconditionally-generated molecules from $P(\bm{x}_{1},\bm{x}_{4})$, and conditionally-generated molecules from $P(\bm{x}_{1}|\bm{x}_{4})$. All distributions are visualized as kernel density estimate plots. All distributions show significant overlap.
Refer to caption: x11.png

### A.10 Additional examples of generated molecules

Figures 22, 23, and 24 show additional random examples of unconditionally-generated molecules and their jointly generated interaction profiles from versions of ShEPhERD trained on ShEPhERD-GDB17. Figures 25, 26, and 27 show additional random examples of ShEPhERD-generated molecules when conditioning on target interaction profiles via inpainting. Although we are primarily interested in using the models trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq for conditional generation tasks relevant to ligand-based drug design, we also show unconditional samples from this model in Fig. 28.

Figure: Figure 22: Random unconditional samples from ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{2})$ on ShEPhERD-GDB17. We show the natively generated molecular structures (without xTB-relaxation) and their jointly generated shapes. We show examples for both large and small $n_{1}$.
Refer to caption: extracted/5947291/figures/appendix_figures/random_unconditional_shape_samples_v1.jpg

Figure: Figure 23: Random unconditional samples from ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{3})$ on ShEPhERD-GDB17. We show the natively generated molecular structures (without xTB-relaxation) and their jointly generated ESP surfaces. We show examples for both large and small $n_{1}$.
Refer to caption: extracted/5947291/figures/appendix_figures/random_unconditional_esp_samples_v1.jpg

Figure: Figure 24: Random unconditional samples from ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{4})$ on ShEPhERD-GDB17. We show the natively generated molecular structures (without xTB-relaxation) and their jointly generated pharmacophores. We show examples for both large and small $n_{1}$. $n_{4}$ is sampled from the empirical data distribution $P(n_{4}|n_{1})$. In these images, grey spheres represent hydrophobes (the most common pharmacophore in ShEPhERD-GDB17); purple and blue arrows represent HBAs and HBDs, respectively; orange arrows represent aromatic groups; and green arrows represent halogen-carbon bonds.
Refer to caption: extracted/5947291/figures/appendix_figures/random_unconditional_pharm_samples_v1.jpg

Figure: Figure 25: Conditional samples from ShEPhERD – trained to learn $P(\bm{x}_{1},\bm{x}_{2})$ on ShEPhERD-GDB17 – obtained by inpainting the 3D molecule given target shapes extracted from molecules held-out from training. We show the natively generated molecular structures (without xTB-relaxation or realignment) overlaid on the target shapes. The shape surfaces are upsampled for visualization.
Refer to caption: extracted/5947291/figures/appendix_figures/random_conditional_shape_samples.jpg

Figure: Figure 26: Conditional samples from ShEPhERD – trained to learn $P(\bm{x}_{1},\bm{x}_{3})$ on ShEPhERD-GDB17 – obtained by inpainting the 3D molecule given target ESP surfaces extracted from molecules held-out from training. We show the natively generated molecular structures (without xTB-relaxation or realignment) overlaid on the target ESP surfaces. The surfaces are upsampled for visualization.
Refer to caption: extracted/5947291/figures/appendix_figures/random_conditional_ESP_samples.jpg

Figure: Figure 27: Conditional samples from ShEPhERD – trained to learn $P(\bm{x}_{1},\bm{x}_{4})$ on ShEPhERD-GDB17 – obtained by inpainting the 3D molecule given target pharmacophores extracted from molecules held-out from training. We show the natively generated molecular structures (without xTB-relaxation or realignment) overlaid on the target pharmacophores.
Refer to caption: extracted/5947291/figures/appendix_figures/random_conditional_pharm_samples.jpg

Figure: Figure 28: Random unconditional samples from ShEPhERD trained to learn $P(\bm{x}_{1},\bm{x}_{3},\bm{x}_{4})$ on ShEPhERD-MOSES-aq. We show the natively generated molecular structures (without xTB-relaxation or realignment) and their jointly generated ESP surfaces and pharmacophores. Note that a number of these structures would not pass our validity criteria.
Refer to caption: extracted/5947291/figures/appendix_figures/random_unconditional_mosesaq_samples_v1.jpg

### A.11 Training and inference resources

Training. We train all models with V100 GPUs (32 GB memory). For both the ShEPhERD-GDB17 and ShEPhERD-MOSES-aq datasets, we train each of the P($\bm{x}_{1},\bm{x}_{2}$), P($\bm{x}_{1},\bm{x}_{3}$), and P($\bm{x}_{1},\bm{x}_{3},\bm{x}_{4}$) models on 2 GPUs for approximately 2 weeks. We train the P($\bm{x}_{1},\bm{x}_{4}$) model on 1 GPU for approximately 2 weeks. We configure the number of GPUs, nominal batch sizes, and the number of gradient accumulation steps to attain effective batch sizes of approximately 48 (molecules).

Inference. For each of the P($\bm{x}_{1},\bm{x}_{2}$), P($\bm{x}_{1},\bm{x}_{3}$), and P($\bm{x}_{1},\bm{x}_{3},\bm{x}_{4}$) models, generating a batch of 10 independent samples (either unconditionally or via inpainting) takes approximately 3-4 minutes on a V100 GPU, when using iterative denoising with $T=400$ steps. Sampling times are reduced by $\sim 50\%$ for the P($\bm{x}_{1},\bm{x}_{4}$) models. Using harmonization when sampling increases the sampling time proportionally to the number of extra denoising steps. On a single CPU, inference requires 5-10 minutes per sample depending on system hardware and the choice of model.

Long sampling times are somewhat intrinsic to the DDPM paradigm due to needing numerous forward passes per sample. Our sampling is also slowed down by our use of memory-intensive E3NN modules throughout ShEPhERD’s denoising network. In this work, we did not attempt to accelerate sampling (e.g., by using fewer denoising steps at inference, reducing $T$, using less expressive networks, reducing $n_{2}$ and $n_{3}$, etc.). We leave such performance engineering to future work.

## Peer Reviews

confidence: 5
contribution: 4
initial_rating: 6
presentation: 4
questions: - To what extent does conditioning with the ShEPhERD method limit diversity of new samples? Could be that previous DDPMs produce on average worse samples but allow for the exploration of more novel chemotypes? (My guess is this method is still very useful but I would like to see this discussed/measured)
- Have you tried an ablation of the joint noising regime? I believe in DiffSBDD its benefits was found to be marginal.
- Could your pharmacophore similarity metric be seen as related to interaction fingerprint similarity? 
- Is it possible to use ShEPhERD completely de novo? I.e. for a new target without example molecules in the PDB? Could pharmacophore features be constructed from hotspot maps/manual feature engineering?
- The fragment merging protocol is nice, can a semi-automated pipeline be provided in the code?
- Section A.9 correctly points out the disparity between the size of ligand seen during training and those seen during evaluation. Why not include larger ligands in the training set?
- Given the large training set size, it would be nice to see some scaling law experiments in terms of dataset and model size. (Given that I am not aware of the resources available to the authors this will not impact my review score - I am just interested if its possible).
soundness: 2
strengths: - The paper is comprehensive and well-written, with excellent figures and molecular visualizations.
- Although conditioning on pharmacophoric features is not new in the field, this is the most thoughtful implementation of the approach to date and a very important direction more works should look into. The authors obviously have a very good understanding of drug discovery and SBDD that is not open seen at ML conferences.
- In my opinion, the field is overly obsessed with only focusing on de novo design models that cannot leverage any target specific data. Therefore, I really like the fact they were able to use experimental fragments screens to get better designed molecules.
- The inclusion of comparisons to REINVENT and random docking screens as baselines is a valuable addition, as these are seldom seen in DDPM literature. (However, I believe there is still room for improvement in the comparisons and evaluations—see below).
summary: This paper introduces ShEPhERD, a generative SE(3)-equivariant diffusion model for molecular drug design that incorporates 3D structural features inspired by medicinal chemistry, namely shape, electrostatics, and pharmacophoric properties. The authors demonstrate ShEPhERD's utility for ligand-based drug design tasks, including ligand hopping, bioactive hit diversification, and bioisosteric fragment merging.

I appreciate the approach of leveraging prior target and chemical information to enhance generative models and commend the authors for advancing in this direction. While the case studies on specific targets are excellent, the general evaluation and particularly the comparative analysis with other methods could be significantly strengthened. I am inclined to give a strong score if these concerns are addressed.
weaknesses: **Main limitations:**

- While the method is strong, main limitation of the paper as it currently stands is the extreme lack of comparisons to related works in the experiments conducted. I kindly request the authors to benchmark against previous works in the relevant tasks. I am not looking for SOTA on all metrics, but better benchmarking is essential to understand the strengths and limitations of the work and I cannot recommend a strong accept at ICLR without more comparisons.
    - There are multiple works in a similar flavour that propose inpainting and/or conditioning on pharmacophore features of some kind. E.g.  DiffSBDD+inpainting [1], MolSnapper [2] and Diffint [3].
    - Compare to DiffHopp [4]/TurboHopp [5] in case of hit diversification. Obviously these models are only in the structure-based context.
    - Compare to ChemProjector [6]/SynFormer [7] incase of ligand hopping.
- Furthermore, the metrics could be greatly improved:
    - Heavy reliance on only Vina scores and shape/tanimoto similarity. 
    - The work cites the PoseCheck [8] paper but does not use any of the metrics. For example, clashes are not measured. In the target conditioned cases, does ShEPhERD provide poses without clashes before relaxation? Is the model able to design hydrogen bonding networks better than previous SBDD models?
    - Strain energy is calculated in A.3 and the model appeared to perform well. However, the impact of this is greatly diminished by not comparing to other works. I would be grateful if the authors could include this to give context. Do the authors use the implentation in PoseCheck or their own? Is the strain energy calculated for poses bound inside the protein pocket?
    - PoseBusters [9] validity could also be used to access the quality of the structures.
- An expanded related work section covering models that perform ligand hopping, scaffold hopping and fragment merging would base the work in the literature much better.
- Would be nice to see anonymous code in the review stage.
- In the hit diversiation section, the authors claim that ShephERD produces enriched sampled “Despite having no knowledge of the protein targets”. While it is technically true the model has never directly see the receptor coordinates, a large amount of structural information is leaked/provided based on the complimentarity of the 7 experimental ligands. This is not a weakness of the method but I believe the language to be misleading. A simple rewording should do.
- The joint model/inpainting method used to condition on the target profile is not novel, as this was done in DiffSBDD 2 years ago [1]. A citation to that work and the original RePaint paper reproducing the ‘replacement’ method of inpainting in the main text would suffice.
- In the interaction-conditioned generation section you only generate 20 samples per target and filter based on interaction similarity. How many samples are left after filtering? How can you ensure statistical significance?
- It is not clear to me if any attempt was made to reduce data leakage between train and test sets (except for size) as the splits seem random. How can you be sure that the target interaction profile on the test samples are not the same/similar to those seem during training?
- Comparison to REINVENT.
    - Was REINVENT trained to convergence? 10,000 oracle calls seems quite low in my experience? (Happy to be proved wrong)
    - Can the authors clarify exactly which version of REINVENT they use? They seem to use an old version from the PMO benchmark. This is the 2017 implementation I believe? I would encourage them to use the largest REINVENT4, this is easy to use and a strong baseline in my experience (https://github.com/MolecularAI/REINVENT4).
    - Little measure of synthetic accessibility in the work. REINVENT does surprisingly well on SA score in my experience. Can this be benchmarks quantitively?

**Minor points:**

- Some background on various inpainting/conditioning strategies for DDPMs would be nice (i.e. replacement method [what is used here], classifier guidance, reconstruction guidance, Doobs-H transform)
- The generated samples in Appendix A.10 are quite small. It would be nice to include some larger samples as well to see how the model handles larger structures.
- The neural network architecture is described in great detail. In my opinion taking too much space. I would prefer to see this in the appendix with more discussion on experiments in the main text.
- What was the intuition for scaling the pharmacophore feature vector by 2.0?
- Can the target profile guidance be more or less ‘weakly’ enforced?
- It would be nice to see the fragment merging work replicated on more than one target. I would recommend the Mpro screen (on Fragalysis already) for example. You can then compare for novelty against the molecules derived from than screen in the COVID moonshot- To what extent does conditioning with the ShEPhERD method limit diversity of new samples? Could be that previous DDPMs produce on average worse samples but allow for the exploration of novel chemotypes? (My guess is this method is still very useful but I would like to see this discussed/measured)
- Have you tried an ablation of the joint noising regime? I believe in DiffSBDD its benefits was found to be marginal. [1]
- Could your pharmacophore similarity metric be seen as related to interaction fingerprint similarity? 
- Is it possible to use ShEPhERD completely de novo? I.e. for a new target without example molecules in the PDB? Could pharmacophore features be constructed from hotspot maps/manual feature engineering?
- The fragment merging protocol is nice, can a semi-automated pipeline be provided in the code?
- Section A.9 correctly points out the disparity between the size of ligand seen during training and those seen during evaluation. Why not include larger ligands in the training set?
- Given the large training set size, it would be nice to see some scaling law experiments in terms of dataset and model size. (Given that I am not aware of the resources available to the authors this will not impact my review score - I am just interested if its possible)..
    - Can you provide a list of the 7 PDBs used to perform these experiments?
- Figure 1 + 2:
    - Are both very good: I would consider merging to safe space and add more analysis/discussion.
- Figure 3:
    - Can ‘dataset samples’ be added to the histograms in the right of the figure?
- Figure 4:
    - The PyMol in this paper is generally very good but the overlapping molecules in this figure is highly confusing. A key/color map would be very helpful. (Especially for hit diversification section).
    - REINVENT molecules seem to consistently outperform ShEPhERD on SA scores. I would be good to include an example structure instead of just listing the metrics.
    - The PDB targets are quite well chosen (Kinases/COVID MPro etc), adding target names in brackets next to PDBs would be helpful.


1. Schneuing, Arne, et al. "Structure-based drug design with equivariant diffusion models." arXiv preprint arXiv:2210.13695 (2022).
2. Ziv, Yael, Brian Marsden, and Charlotte M. Deane. "Molsnapper: Conditioning diffusion for structure based drug design." bioRxiv (2024): 2024-03.
3. Sako, Masami, Nobuaki Yasuo, and Masakazu Sekijima. "DiffInt: A Pharmacophore-Aware Diffusion Model for Structure-Based Drug Design with Explicit Hydrogen Bond Interaction Guidance." (2024).
4. Torge, Jos, et al. "Diffhopp: A graph diffusion model for novel drug design via scaffold hopping." arXiv preprint arXiv:2308.07416 (2023).
5. Yoo, Kiwoong, et al. "TurboHopp: Accelerated Molecule Scaffold Hopping with Consistency Models." arXiv preprint arXiv:2410.20660 (2024).
6. Luo, Shitong, et al. "Projecting Molecules into Synthesizable Chemical Spaces." arXiv preprint arXiv:2406.04628 (2024).
7. Gao, Wenhao, Shitong Luo, and Connor W. Coley. "Generative Artificial Intelligence for Navigating Synthesizable Chemical Space." arXiv preprint arXiv:2410.03494 (2024).
8. Harris, Charles, et al. "Posecheck: Generative models for 3d structure-based drug design produce unrealistic poses." NeurIPS 2023 Generative AI and Biology (GenBio) Workshop. 2023.
9. Buttenschoen, Martin, Garrett M. Morris, and Charlotte M. Deane. "PoseBusters: AI-based docking methods fail to generate physically valid poses or generalise to novel sequences." Chemical Science 15.9 (2024): 3130-3139.

---

confidence: 5
contribution: 4
initial_rating: 10
presentation: 4
questions: Although it is not unreasonable, what is the logic for fitting the surfaces onto the volumes in order to obtain the parameters required for the method?  Also, why only fitting the particular set of discrete numbers of points (50, 100, 150, 200, 300, 400) and not, say, a function of the number of atoms, or of the total surface size?  Would there be a problem with scaling the method to a large number of points, or is it simply deemed unnecessary?  As a minor related question, what is the purpose of a separate set of surface points for x2 compared to x3?  Wouldn't the coordinates of x3 suffice for both tasks?

Other than following prior work, do the authors have a high level rationale for selecting the particular subset of targets from the PDB in Zheng et al (2024)?
soundness: 4
strengths: The paper presents a substantial advance compared to previous related methods; it has a strong logical foundation with clear exposition of the method and its parameters, and a diverse set of demonstrations.  The conditional generation of molecules based on ligand-only properties is important in applied drug discovery and the method offers both a practical way for thinking about some of the concepts and a tool that can be readily applied.  In addition, it introduces useful and practical datasets, and it defines a simple but clear metric for comparing surface shapes.  I think that the work opens up a logical framework for building practical tools in drug discovery and will influence followup research.
summary: This paper presents a new method for the generation of small molecules conditioned on shape, electrostatic potential surface, or pharmacophores.  This development provides a useful and practical tool that can complement multiple conventional tasks during various stages in drug design projects.  The method is a joint diffusion model that predicts molecules in terms of graphs plus coordinates, surfaces sampled on a modest number of points (with or without an electrostatic potential), and pharmacophores (types and coordinates).  The conditional generation uses inpainting to sample the molecule graphs and coordinates conditioned on one or more other properties of the molecule, which can mimic a known target molecule.
weaknesses: Most of the weaknesses of the paper are somewhat superficial and do not limit the importance of this contribution.

All the conditional generation requires providing n1 a priori for each sample, and the pharmacophore application, which is particularly interesting and useful, requires providing both n1 and n4.  In addition to sampling n4 based on the chosen n1 and the training data distribution, have the authors considered exploring ways to automate or predict values for parameters like n1 and n4?  It would be very interesting to explore the regime where compute limitations are of no concern.

The demonstration of conditional generation in Fig 3 is performed on a random subsample of the training molecules.  I understand that there is no perfect split of molecules, however, the jump from a random subset, which work rather well, to natural products, which are completely out of distribution and don't appear to work well (albeit better than reinvent), appears a bit too dramatic.  To enhance the robustness of the evaluation and to gain additional insights, the authors could create additional, harder, splits by molecular properties such as numbers of atoms or rings, splits by scaffold, or splits by similarity clustering of the training data. Because the authors provide two new training datasets, it is perhaps more important than it would otherwise have been to provide and characterize a number of distinct, increasing difficulty splits.

The surfaces were fit on modest numbers of points and the graph embedding networks are surprisingly small, which are all probably a result of trying to keep the graph neural nets from collapsing.  It is possible that these shortcuts (and the usual problems of graph neural nets and scaling to larger systems) lead to the size-dependent effects demonstrated in the supplemental materials; these limitations might also underpin the examples of natural products.  It may be worth adding a brief note to discuss the scaling of the method with the size of the small molecules, even if the detailed analysis is reserved for the supplemental materials due to space limitations.

---

confidence: 4
contribution: 3
initial_rating: 6
presentation: 3
questions: My questions are the same as mentioned in Weaknesses.
soundness: 3
strengths: 1. The method provides a comprehensive summary of related work. To my knowledge, it is currently the only 3D generative model that simultaneously incorporates shape and pharmacophore constraints.
2. The authors achieved credible results in three specific tasks, partially demonstrating the method’s effectiveness.
summary: The paper introduces the ShEPhERD method for generating 3D molecules that satisfy specific shape, electrostatic surface, and pharmacophore constraints. The method represents molecular shapes, electrostatic surfaces, and pharmacophores as point clouds and jointly models these properties with the 3D molecular structure. During sampling, inpainting techniques are used to generate molecules that meet the specified properties. This method enriches ligand-based drug design approaches and shows significant innovation in its application areas. However, it lacks substantial innovation in the basic structure of the generative algorithm.
weaknesses: 1. The method represents shapes and electrostatic surfaces as point clouds and generates them together with the molecular structure using 3D graph neural networks, which may significantly impact computational efficiency. The generative model based on DDPM may also have efficiency drawbacks. The authors should provide relevant statistical results.
2. The method lacks comparisons with other methods, casting doubt on the superiority of the results.
    (1). For shape-based molecular generation, methods like ShapeMol and SQUID can achieve similar effects and should be compared. (2) For pharmacophore-based molecular generation, both pharmacophore type matching and structural matching should be considered and evaluated separately. If combined comparison is necessary, reference to 2D pharmacophore works such as the Match Score (Nature Communications | (2023)14:6234) is suggested. (3) The authors should provide the similarity of generated molecules to reference molecules for each specific task to ensure that the generated molecules are not simple variations of the reference molecules.

Meta-review: In this work, authors introduce ShEPhERD, an SE(3)-equivariant diffusion model for 3D molecular generation that jointly models molecular graphs, shapes, electrostatic surfaces, and pharmacophores. The model enables conditional generation of molecules with desired 3D interaction properties through inpainting, facilitating tasks like natural product ligand hopping and bioisosteric fragment merging.

The main scientific contribution is demonstrating that a joint diffusion model over multiple molecular representations enables generating novel molecules that preserve key 3D interaction patterns of reference compounds. The authors validate this through experiments showing ShEPhERD can: 1) generate molecules with similar shapes to references while maintaining low graph similarity, outperforming SQUID as shown in additional experiments, 2) generate bioisosteric analogs of natural products that maintain interaction patterns while having lower synthetic accessibility scores, and 3) merge fragments based on their pharmacophoric interactions.

The reviewers identified several key strengths. Reviewer HW8X highlighted the thoughtful implementation of pharmacophore conditioning and practical utility for drug discovery applications. Reviewer pbc2 noted the strong logical foundation and clear exposition. The extensive comparisons to baselines like REINVENT and response to reviewer requests for additional comparisons to methods like SQUID, DiffSBDD, and SynFormer demonstrate robust performance.

Some limitations were raised, particularly around computational efficiency and scaling to larger molecules, as noted by Reviewer JQjD. However, the authors acknowledge these engineering challenges as future work, appropriately focusing first on establishing the method's utility. Questions about dataset splits and evaluation protocols were addressed through additional experiments showing good performance on challenging out-of-distribution tasks.

I recommend acceptance for several reasons: First, the paper makes a clear technical contribution by introducing a novel joint modeling approach for 3D molecular properties. Second, the extensive experimental validation, including new comparisons added during discussion, demonstrates practical utility across multiple drug design tasks. Third, the method enables important capabilities not possible with existing approaches, particularly around preserving 3D interactions in ligand-based design. While there are opportunities for improving efficiency and scaling, these do not diminish the core methodological advance. The thorough author response and additional experiments effectively addressed the key reviewer concerns about comparisons and evaluation protocols.

In summary, this work represents an important contribution to molecular generation methods that will benefit both the machine learning and drug discovery communities. The clear exposition, comprehensive evaluation, and demonstrated practical utility make it worthy of presentation at ICLR.

Decision: Accept (Oral)