ShEPhERD: Diffusing shape, electrostatics, and pharmacophores for bioisosteric drug design
Authors: Keir Adams
Kento Abeywardane
Jenna Fromer
Connor W. Coley
Venue: ICLR 2025 Oral (2025)
Keywords: 3D molecular generation
drug design
molecules
Primary area: applications to physical sciences (physics, chemistry, biology, etc.)
TL;DR: We design a diffusion model that jointly generates 3D molecules and explicit representations of their 3D shapes, electrostatics, and pharmacophores and demonstrate its utility in bioisosteric drug design

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, 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.

Published as a conference paper at ICLR 2025
ShEPhERD: DIFFUSING SHAPE, ELECTROSTATICS, AND PHARMACOPHORES FOR BIOISOSTERIC DRUG DESIGN
Keir Adams", Kento Abeywardane Jenna Fromer, & Connor W. Coley Massachusetts Institute of Technology, Cambridge, MA 02139, USA (keir, kento, j Fromer, ccol ey emii fmmit. edu
ABSTRACT
Engineering molecules to exhibit precise 3D intermolecular interactions with their environment forms the basis of chemical design. In ligand-based drug de- sign, bioisosteric analogues of known bioactive hits are often identified by vir- tually 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 inter- action profiles may facilitate 3D interaction-aware chemical design. We specif- ically design ShEPhERD an -Euiivaquv diffusion model which jointly diffuses/denoises 3D molecular graphs and representations of their shapes, elec- trostatic potential surfaces, and (directional) pharmacophores to/from Gaussian noise. Inspired by traditional ligand discovery, we compose 3D similarity scoring functions to assess ShEPhERD's s ability conditionally generate novel molecules with desired interaction profiles. We demonstrate ShEPhERD's potential for im- pact via exemplary drug design tasks including natural product ligand hopping, protein-blind bioactive hit diversification, and bioisosteric fragment merging.
INTRODUCTION
Designing new molecules to attain specific functions via physicochemical interactions with their environment is foundational task across the chemical sciences. For instance, early-stage drug dis- covery often involves tuning the 3D shape, electrostatic potential (ESP) surface, and non-covan non-covalent interactions of small-molecule ligands promote selective binding
protein target (Bissantz) jet al. 2010 Huggins et al.] 2012). Homogeneous catalyst design requires developing organometal- lic, organic, or even peptidic catalysts that stabilize reactive transition states via specific noncovalent interactions (Raynal et al.
et al. 2020 Knowles & JJcobsen] 2010 Wheeler et al., 2016 Toste et al.) 2017 Metrano Miller 2018). Supramolecular chemistry similarly optimizes host-guest interactions for applications across photoresponsive materials design, biomedicine, and structure-directt zeolite synthesis (Qu et al. 2015) Stoffelen Huskens 2016/Corma al., (004).
This essential challenge of designing molecules with targeted 3D interactions manifests across myr- iad tasks in ligand- drug design. In medicinal chemistry, bioisosteric seaffold hopping aims to swap-out substructures within larger molecule while preserving bioactivity (Langdon et al.) (000)). Often, the swapped scaffolds share biochemically-reeewee features such as electrostatics or pharma- cophores, which describe both non-directional (hydrophobic, ionic) and directional (hydrogen bond acceptor/donor, aromatic π--,, halogen bonding) non-covalent interactions. When scaffold hopping extends entire ligands, "ligand hopping" can help identify synthesizable analogues of complex natural products that mimic their 3D biochemical interactions (Grisoni et al. 20188). In hit expan- sion, ligand hopping is also used to diversify known bioactive hits by proposing alternative actives, ranging from topologically-similal "me-too" compounds to distinctly new chemotypes (Schneider) al. 2006 Wermuth| 2006). Notably, ligand hopping does not require knowledge of the protein target. Lastly, bioisosteric scaffold hopping extends beyond altering individual molecules; Wills al. (202)) used "bioisosteric fragment merging" to replace a ser of fragments that independently bind protein with one ligand that captures the fragments' aggregate 3D binding interactions.
"These authors contributed equally
ShEPhERD: Shape, Electrostatics, and Pharmacophore Explicit Representation Diffusion

---

Published as a conference paper at ICLR 2025
P(
ShEPhERD
P(X
Applications in bioisosteric ligand design
natural product ligand hopping
bioactive hit diversiffiation
bioisosteric fragment merging
molessiar
shape
pharmucopheee Ietegaake
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 nat- ural product ligand hopping, bioactive hit diversification, and bioisosteric fragment merging.
Figure/ 1:
Here, we consider three archetypal tasks in bioisosteric drug design: 3D similarity-constrae lig- and hopping, protein-blind bioactive hit diversification, and bioisosteric fragment merging (Fig. [ 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, tradi- tional 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 ref- erence molecule (Oprea & Matter 2004, Rush et al., 2005, Zavodszky et al. 2209) [anddrss et al.) virtual screening has acute drawbacks: it cannot explore beyond
2012). However, similarity- known chemical libraries, is inefficient by virtue of being undirected, and it can be quite slow when multiple geometries of conformationnally-eeiii flexiblle molecules must be scored.
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 in- teractions. (2) Harris et al. (2023) have found that 3D generative models for structure-based drug design struggle to generate ligands that form biochemical interactions (c.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 set- tings are often data-restricted, necessitating zero/few-shot generative approaches. To address these challenges, we introduce ShEPhERD, 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 phar- macophores that are amenable to
3E(3)-equivariant diffusion modeling.
We formulate a joint denoising diffusion probabilistic model (DDPM) that learns the joint distribu- tion 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 elec- trostatics, natively model the directionality of pharmacophores by diffusing/denoisingg vectors on the unit sphere. We sample from specific interaction coniatiioonconiiionnee distributions via inpainting.
Inspired by virtual screening, we craft shape. electrostatic, and pharmacophore similarity func- tions to score (1) the self- consissnist 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- -iirraaion-smmilaiity to target profiles, even upon geometrie relaxation with semi-empirical DFT.
After training on drug-like datasets, we demonstrate in vilico that ShEPhERD facilely design small-molecule mimics of natural products via ligand hopping, diversify bioactive hits while pre- serving protein-binni 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 cam- paigns that require the de novo design of new molecules with precise 3D interactions. However, we stress ShEPhERD's 's general applicability to other areas of interaction-awarr 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

---

Published as a conference paper at ICLR 2025
2 RELATED WORK
3D similarity scoring for ligand-based drug design. Ligand- 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. 2019 Rush et al.) 2005 Jack-
son et al. 2022). Shape similarity functions typically use dtttecentered Gaussians to compute the
volumetric overlap between active and query molecules (Grant Pickup 1995:
Grant et al. 1996).
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., 1992. Bolcato et al.
2022
OpenFye, Cadence Molecular aar Sciences (024) or pharmacophore similarity (Sprague) 1995 Dixon ct 20066] Taminau et al. 2008, Wahl 2024). To better capture intermolecular interactions, other methods quantify the Coulombic pharmacophoric "potential" felt by chemical probe at points near the molecule's 's surface (Cheeseright et al. 2006 Vainio et al.; 2009) Cleves et al. 2019). We similarly use surface representations of molecular arr 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 r-effects (Wahll 2024]. Notably, we de- velop our own 3D similarity scoring functions which natively operate on our chosen representations.
Symmeetry-reeeerving generation of molecules in 3D. Many generative models for small molecules have been developed to sample chemical space (Gómez-Bombarelli et al., 2018 Segllr al. 2018 Jin et al. 2018 Vignac al. [022] Bilodeau et al. [022] Anstine & Isayev (2023). 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 gen-
erating structures with E(3)- or SE(3) (in3)-inna internal coordinates (Gebauer et al.] [202] Luo & Ji. 2022 Roney et al.] 2022). or more recently, by generating Euclidean coordinates with equivariant networks (Vignac al. 2023 Peng et al. 2023 Irwin al. 2024). We generate molecular struc- tures with an equivariant DDPM (Ho et al. 2020 Hoogeboom et al. 2022), a strategy that has seen use in nooein-connitioned ligand and linker design (Igashov et al.)
2024 [Schneuing et al. 2022).
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)-equivarian Euclidean Neural Networks (E3NNs) (Geiger & Smidt, (2022) encode/decode coordinates, vectors, and scalar features across our heterogeneous representations.
3D interaction- aware molecular generation. Prior works have partially explored shape- and pharmacophore-uun aware molecular generative design. Multiple methods generate either (ID) SMILES or (2D) molecular graphs conditioned on 3D representations of target shapes and/or pharmacophores (kaalic al., 2019 Imrie et al. 2021 Zhu et al. 2023, Xie et al. 2024). or use 3D similarity scores to fine-tune unconditional generative models (Papadopoulos et al. 2021 Neeser et al. 2023). which we include as baseline (App. [..]. Since these generative models do not directly predict 3D struc-
tures, they require conformer generation as a post-processing step. On the other hand, Adams & Co- ley (2022) found that generating structures natively in 3D yields more chemically diverse molecules with higher 3D shape-similarity compared to competing shape- canee-conditioned ID approach. Mul- tiple other methods have since been developed for chape-conditioned 3D generation (Chen et al.]
2023 Lin et al. 2024 Le al., 2024). Regarding pharmacophores, Ziv et al. (2024) applied 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 т-- interactions, hydrophobic effects, and halo- gen bonding. Electrostatics-awarr molecular generation has been considered only by [olcato et al. (2022), via exchanging eneeenmmeeteed chemical fragments with similar electrostatics. But, they do not conditionally generate 3D molecules given global ESP surface. Our work unifies and extends prior work on interaction-aware 3D generative design by comprehensively modeling shape, electro- statics, and arbitrary pharmacophores (including their directionality) a single general framework.
3D generative structure-based drug design (SBDD). Tangential to our work are models that gen- erate ligands inside a protein pocket by training on potein-ligand, complexes with (Lee et al. 2024 Sako et al. [2224] Zhung et al. 2024 Wang et al. 2022. Huang et al. 2024) or without (Peng al. 2022 Schneuing et al.) 2022 Guan et al. 2023) explicit encodings of their interactions. In contrast, 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

---

Published as a conference paper at ICLR 2025
ShEPkERD 3D Similarity Scoring Functions
Overlap weighted by electrostatic potential
Point Similarity
Isotropic Gaussian Overlap in I
Son Sen-dheetoood
phonauopione
Marched chargey relatn everlop
Mintsatched charges
decreace overlap
Overlap weighted by pharmacophore types
High overlap
overlup
A
Directionless relain eveelop
Mismatched types have so overlop
Miseligned dicections dierease overlap
ShEPhERD Denoising Network
(--1) i
KINN
i-1
interact node features
ator types focmel charges bond types coordinates
coondinates
2
coondinatees
aD grapa
eleetrostatic potencial
Cg
pharm. types
4
directions
Input States
Embedding Modules
Joint Module
Denoising Modules
Conditional Sampling via Inpainting
Replace jointly denoised Interaction prodile(s) with nolsed targes profie((s
noising target interaction profle(s)
ShOPkERD
SUPPLERD
m
m
exploxed
polly sarget grotie
nolay
ExEget peodle
dicctbuation
Figure 2: (Top) Visualization of our 3D point-cloud similarity scoring functions, which evaluate sur- face/shape, electrostatic, and pharmacophore similarity via weighted Gaussian overlaps. (Middle) ShEPhERD's 's denoising network architecture, which uses SE(3)-equivarian neural networks to (I) 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 molecules that exhibit desired interaction profiles.
METHODOLOGGY
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 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-conditionee distributions over 3D chemical space via impainting. Below, formally define our chosen representations for 3D molecules and their interaction profiles, and also define the 3D similarity scoring functions used our evaluations. We then detail our joint DDPM and sampling protocols. Fig. 2 visualizes our overall methodology.
3.1 DEFINING REPRESENTATION OF MOLECULES AND THEIR INTERACTION PROFILES
3D molecular graph. We define organic molecule with mi atoms in specific conformation as a 3D molecular graph æ; (a, f, J,C. C, B) where a (H, C,O.... "h lists the atom types (N, options); f (-2, -, 2 specifies the formal charge of each atom; C' [ specifies the atomic coordinates; and B {0, 2,3, 1.5]")
is the adjacency matrix specifying the covalent bond order between pairs of atoms. For modeling purposes, the categorical variables a, f. B are represented as one-hot continuous features so that a R x f E R=1x6 and B E Reixn1xb
Shape. We define the shape 2e S2 E X3233 the solvent-ace accces-leeeesiblee surface of 21 n2 is fixed regardless of n, We use surface points to better decouple the shape representation from the molecule's 's exact coordinates or 2D graph A.5).
point cloud with n2 points sampled from
4

---

Published as a conference paper at ICLR 2025
Electrostatic potential (ESP) surface. We represent the ESP surface 23 = (S3, v) by the Coulom-
bic potential at each point in a surface point cloud S3 E R=*x3, with v E RRa,, We compute v from the per-atom partial charges of 11, which are computed via semi-empirical DFT (App. A.5
Pharmacophores. The set of n4 pharmacophores are represented as 24 = (p, P.V) where TR is matrix of one-hot encodings of the Np pharmacophore types; P RR443 is the pharnacophores' coordinates; and V { 0) contains the relative unit or zero vectors speci- fying the directionality of each pharmacophore. The directional pharmacophores include hydrogen bond acceptors (HBA) and donors (HBD), aromatic rings, and talogen-carbon bonds. Direction- less pharmacophores (V(k = 0) include hydrophobes, anions, cations, and zinc binders. App. AA.5 details how we extract pharmacophores from molecule by adapting known SMARTS patterns.
3.2 SHAPE, ELECTROSTATIC, AND PHARMACOPHORE SIMILARITY SCORING FUNCTIONS
To formulate scoring functions, first define point cloud Q R x3 where each point rи
is assigned an isotropic Gaussian in IR3 (Grant & Pickup 1995 Grant et al. 1996). We measure the
Tanimoto similarity sim E [0, 1] between two point clouds QA and QB B using first order Gaus- OA,B EnEQu EsEQa wa.b( exp (
sian overlap OA."
1-221 arrra r)ll3) and sim" (QA. QB)
OA.A a O where is a Gaussian width and wab is a weighting factor. Note that 3D similar- ities are sensitive to SE(3) transformations of Q with respect to Qg. We characterize the similarity at their optimal alignment by sim A QB) max sim" (RQ) + t. Qx) where R E SO(3) and t E T(3). We align using automatic differentiation. Note that niq, need not equal nQn.
Shape scoring. The volumetric similarity between two atomic point clouds CA and C3 is defined as sim (CA, Cr) with Wa,b = 2.7 and a = 0.81 (Adams & Coley| 2022). We newly define the surface similarity between two surfaces S^ and Sa sim (S SB) with W 1, and = (n2). Here, V is a function fitted to sim depending on the choice of n2 (App. A.6)
Electrostatic scoring. We define the similarity between two electrostatic potential surfaces 33,4 and ø3,B ( B3 with Wa,s exp (- [| ), (na), and À
0.3
[Irc) y
Eo is the permittivity of vacuum with units e²v - A)-1. Inspired by Hodgkin & Richards (1987) and Good (1992), we use the difference between electrostatic potentials increase sensitivity to their respective magnitudes (i.e., rather than simply comparing signs).
Pharmacophore scoring. We define the similarity between two sets of pharmacophores æ4,4 and Emem OA,Bim
F4,8 with phamm ( (I FA.B) Emere A.. 8:m -OA.a.. Here, M is the set of all phar- macophore types (IM) Np). Wa,b = 1 if m is non-directional, or scaling of vector cosine sim- ilarity Wa.b.m = Vial +2 if m is directional (Wahl, 2024). am = 9(m) where S2 maps each pharmacophore type to Gaussian width (App. A.6) We take the absolute value of V [a] V(b).. for aromatic groups as we assume their ı 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. 2020 Sohl-Dickstein et al.
2015) to decom-
pose the joint distribution over the tuple X (x1, 22, 23, 2.4) as Pamal (X) = P(X
(TX(T) II= (1X( x()), where Pama is the data distribution, P P(X (T) is (roughly) a Gaus-
sian prior, and (t-1) X (t) are Markov transition distributions learnt by a neural network. This
network is trained to reverse forward-noising process P(X x
(o, I
which gradually corrupts data X into Gaussian noise (T) according to variance preserving noise schedule given by of and a; vi - of for t 1, .... T. See Ho et al. (020) for preliminaries on DDPMs. Here, we describe the forward and reverse processes ShEPhERD's joint DDPM.
Forward noising processes. We follow Hoogeboom et al. 2002) to forward-noise the 3D molecule I1 = (a. f. C)B For a E "N,, we use Gaussian noising where (t) = (t-1) + dre for (N, 1). The processes for f and B are similar, but we symmetrize the upper/lower triangles of (()), We apply isotropic noise to C R x3, but center the noise at 0 to ensure translational
invariance: ( og ((-1) oe(E 1 S {(k)) for <(k) ~ N(0, Iạ) and e x.
5

---

Published as a conference paper at ICLR 2025
For the molecular shape 22 S2 E R=x we also forward-noise with isotropic Gaussian noise: 5 o + diC for ~ ((, We do not subtract the noise's center of mass (COM), though; this ensures that the model can learn to denoise 22 such that it is centered with respect to 11- For the ESP surface 33-- = (S3, v), we forward-noise the surface Sa E Raax3 x3 in the same manner as S2. We forward-noise byy R in the typical way: v(4) = u. + orE for 6 ~ NO0, I).
For the pharmacophores F4 (p, P,V), we forward-noise their types p just like the atom types a, and their positions P just like the shape S2. Diffussng the vector directions V E { 0}"* is complicated since some pharmacophores are directionless. To unify their treatment, we interpret the pharmacophore vectors as Euclidean points in R3 (e.g., only noiseless vectors have norm 1.0 or 0.0). We then forward-noise the vectors like any point cloud: V = ar V(t + orE for E ~ N( 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 (*) = a a (0) T,E for e ~ N(0, I). where a; = II'=1 a, and o = VI-a?
Reverse denoising process, Starting from any X (t) (but typically pure noise X (2') X(') ~ N NO, f)), the DDPM iteratively denoises X() by stochastically interpolating towards a predicted clean structure X00 se x0) ~ Patta (x (t)) resembling true samples from the data distribution. We follow the DDPM formulation where rather than predicting X (0) directly, the network predicts the true noise g(t) ~ 6 that, when applied to data X(), yields (() = C2 xx,X (0) + ore. In this case, the single-step denoising update can be derived as: (t- a, (t) of e(t) d c', where the additional noise e' ~ N (0, I) (set to E' = for t = 1) makes each denoising step stochastic.
x(0)
a.
ai
ShEPhERD employs a single denoising network 7 that is trained to jointly predict the noises _22 3z
e) a( ee()) e
=n( (t) (t)
44 .t). where e(*) (El" el() e( (e).).))... e (E
al0) ($.) el)). and 2) (E
(()) al)) At inference, we jointly apply the denoising updates EEp
to obtain ((-11), When computing 1
a(t) (t-1). we remove the COM from & and the extra noise 6G-
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 point-cloud with one-hot, scalar, and/or vector attributes. One can also directly model specific marginal distributions (e.g., P(x1, x2),, P(æ1, #2y), or P(11, .21)) )) by simply modeling a subset of the variables ( {æ22, 32, .. Finally, ShEP- hERD 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 1 to satisfy three criteria:
The noise predictions E(t) e E E are T(3)-invariant and SO(3)- equivariant with respect to global SE(3)- transfrttooos of (t) in order to (I) efficiently preserve molecular symmetries, and (2) ensure x x22 æz and x remain aligned during denoising.
Symmetryy-ressevvingg
(2) (2) (t)
(t)
Expressive: v captures both local and global relationships between æ1 : 222 203 (t) and x!"
General: To promote applications across chemical design, r) accommodates other definitions of shapeddletrosatii//phaammccccc
and permits incorporating other structural interactions, too.
To achieve these criteria, we design 77 to have three components: (1) a set of embedding modules which equivariantly encode the heterogeneous æ; into uniform sets of latent (-0 (scalar) and 1-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 é, for each x;,
The embedding modules use SE(3) equivvriaat E3NNs (we choose to use expressive EquiformerV2 modules (Lao et al.
, _x(" and æ into latent codes (z] z(11 = (1 V E (1,2,3,4] 4]. Here, z, E R., xd are invariant scalar representations of each node (e.g., atom, point, or pharmacophore), and ž, E xnxx are equivariant vector representations of each node. To make each system x; sensitive to rela- tive translations with respect to x
) to individually encode the 3D structures of (t) (t)
we also include a virtual node that is positioned at the center of mass of zz and which remains unnoised. Prior to 3D message passing with the E3NNs, scalar

---

Published as a conference paper at ICLR 2025
atom/pointyyharmscohoor features (e.g., a f(0, v(t), ( are embedded into 1=0 node features, and vector features (e.g.. V(0) are directly assigned as 1=1 node features. The pairwise bond rep- resentations (+) for æ (4) are also embedded as l=0 edge attributes. Finally, for each Φi, we embed sinusoidal positional encodings of the time step t and add these to all the I=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 x 25 z3 and x to form a heterogeneous 3D graph where the nodes are attributed with their corresponding latent features (z x V E 21, 3, 4]. We then encode this heterogeneous 3D graph with another E3NN (EquiformerV2 V2) module fomm and resid-
ually update the nodes' latent features: (z.") 2((1) += o glocal
(Ir! 2 (0)
((() Viv [1, (1,
(2) We sum-pool the updated 1=1 node features across each sub-graph, concatenate, and then embed with equivariant feed-forward network @gglahel to obtain global 1=1 code describing the overall Sglobel
(EX= - [b]) ViE E [1, 2, 4) We then apply equivariant tensor products between =[(0 E Rdx3 and an /=0 embedding tf This yields /=0 and (=1 global latent features (z joins* joint), which are residually added the node representations (z (t) (,()).
system: joint =
joint
Cat
(1) ( For 22. 33. The denoising modules predict the noises é, from the node-level features (z,
and 1s, the scalar noises (E(") ) are predicted from the corresponding (1=0)
z(() codes using multi-layer perceptrons (MLPs), whereas the vector noises (e(s). e e(e) e( are predicted from the (1=1) z(!) codes using equivariant feed-forward networks ("E3NN- style" coordinate predic-
tions). For I Ea and e are predicted from (t) with simple MLPs.
e is predicted from pairs (=("]/]], (1) [i]) using permutation-invariant MLP (App. A.7.1). e is predicted from (z(0) x(9)
using E3NN-style and EGNN-style (Satorras et al. 2021 coordinate predictions (App. A.7.I
Training. We train the denoising network 77 with unweighted L2 regression losses 4, 1| e, 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(=.). ), P(x11 222). P(x1, 233) P(x1 x24, and P(z1, D3, za)- Note that since 23 defines an (attributed) surface S3, jointly modeling (I2g, æ23) is redundant; 23 implicitly models Ø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 (T) from isotropic Gaussian noise, and then denoise for T steps to sample ( We then argmax a (0), f()), B()), and p()) to obtain discrete atom/bond/pharm. types, and round each V [k] to have norm 1.0 or 0.0. To help break the spherical symmetry of (t) at early time steps, we strategically add extra noise to the reverse process (App. (A.)). For conditional generation, we use inpainting (Lugmayr et al.) [022] Schneuing al. 2022) to sample _ conditioned on target interaction profiles. Namely, we first simulate the forward -nissng of the target profiles (at), æ3,, xi) + (az 2q. x1)(t) t E [1. Then, during the reverse denoising process, we replace the ShEPhERD-denoise (az , ,z3
(.(()) (() with the noisy
target (x)) zg (f) Like other molecule DDPMs, we must specify n. and n4 for each sample.
4 EXPERIMENTS
We train and evaluate ShEPhERD using two new datasets. Our first dataset (ShEPHERD-GDBI7 contains 2.8M molecules sampled from medicinally-relevant subsets of GDB17 (Ruddigkeit et al., 2012 Awale et al.) 2019 Bühlmann & Reymond, 2020). Each molecule contains 17 non-hydrogen atoms with element types in (H, N, O. S. F, CI, Br, 1}. and includes one conformation opti- mized with GFN2-xTB in the gas phase (Bannwarth et al. 2019). Our second dataset (ShEPhERD- MOSES-aq) contains 1.6M drug-like molecules from MOSES (Polykovskiy et al., 2020) with up to 27 hyd-hedr atoms. Each molecule contains one conformation optimized with GFN2-xTB in implicit water. In all experiments, hydrogens are treated explicitly. Whereas we use ShEP. hERD-GGDI7 evaluate ShEPhERD in straightforward unconditional and conditional generation
7

---

Published as a conference paper at ICLR 2025
Unconditional Generation self consistency of generated shapes
Conditional Generation
Diaritztion of signed dnilatty screvs
Diurhuztins sllyed comple
camples
529
Buating Pancie
Shes Smilltt
SUPERR
self-consistency of penerated ESP surfaces
shape-condirtoned samples
SMBESERD
Pe.200 Shart
0.90
0.87
samples
Rngeerdii maniaa Sirlvitty
self consistency generated plarmacoghores
SREPEERD vensicn
skartano sapple
reacouataccl jua rpherecttteed
0.84
cinmaauohorr-ccelttiion samples
0.61
0.73
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 ShE.PhERD-geneeee dataset-sampled 3D molecules (post-relaxation and realignment) and 100 target molecules from ShEPhERD-GDBIT 17, 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-ate RMSD upon xTB-relaxation. Also shown are top- rop-scoring samples overlaid on their target profiles (surfaces are upsampled for visualization), labeled with 3D similarity scores.
or
settings, we use ShEPhERD-MOSER aq to challenge ShEPhERD to design drug-like analogues of natural products, to diversify bioactive hits, and to merge fragments from a fragment screen.
Unconditional joint generation. We first evaluate ShEPhERD's s ability to jointly generate 3D molecules and their interaction profiles in self-consistent way. Namely, 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. reports distributions of the 3D similarity between gen- erated and true interaction profiles across 1000 samples (with n4 [11, [1], n44 Pauta (n4 (n)) from models trained on ShEPhERD-GDBII to learn P(x x22, P(x1,.3)) or P(x1, x. When compare against the (optimally aligned) similarities simsarss simEsP, and simpharm between the true profiles of the generated molecules and those of random molecules from the dataset, we con- firm that ShEPhERD'' generated profiles have substantially enriched similarities to the profiles in all cases. Interestingly, ShEPhERD is more self-consistent when generating shapes or ESP sur- faces than when generating pharmacophores. We partially attribute this performance disparity to the discrete nature of the pharmacophore representations and the requirement of specifying ny; ShEPhERD occasionally generates molecules that have more true pharmacophores than generated pharmacophores, as demonstrated by the samples shown in Fig. (i.e., some have missing HBAs).
Interaction-condditonne generation. We now evaluate ShEPhERD's ability sample from interaction interaction-condiiionee chemical space. To do so, we reuse the same P(æ1, Z2), P(æ1, 23), and P(æ1, .:4) models trained on ShEPhERD-GDBIT but use impainting to sample from the interaction- conditioned distributions P(æ1 (z2)). P(z, lar3), and (a4) Specifically, we extract the true in- teraction 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
but
We 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 A for all unconditional models trained ShEPhERD-GDBI7. App. A.4 reports such unconditional metries (validity, novelty, etc.).
8

---

Published as a conference paper at ICLR 2025
Natural product ligand hopping
Bioactive hit diversification
Natural product
ShEPhERD analogues (SAU) ESP sim Phann.
Socree tscatrres
Xo deiigg
033
37, 0422,
Protein Terget
Bot 033
PDB Liqand
Top analogues
Rocuark Lebray
Vias
Vie 1084
Vaar-111
Mine J0B
7133
SA:
Vila.
Top
docked poses
-0.3
-1o
160
SAx 6.6
Bost
Smauk
701
71L1
Pragment screen
Bioisosteric fragment merging
Overiaid SREPİERD poses (rckosed & realignedl)
(13 everlapinng Iragnenss)
SNEPAERD
pharmacophoers
XA 0.42
0.57,
0.30
Figure 4: (Left) Examples of ShEPhERD-generane analogues of natural product targets, labeled by SA score, ESP similarity, and pharmacophore similarity to the target. Similarities are computed after XTlrrlaxationn and ESP-optimal realignment. (Right) Distributions of Vina scores for <000 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-ag For 5mo4 and 7111, we show top-scoring ShEPhERD -egnrrated lig- ands (conditioned on low-energy poses), and overlay selection from their top- 10 docked poses on the PDB ligands. (Bottom) ShEPhERD's s bioisosteric fragment merging workflow. We extract ESP surface and pharmacophores of 13 fragments from fragment screen, and show ShEPhERD- generated ligands with low SA score and high 3D similarities to the fragments' interaction profiles.
2D graph similarities >0.3 to the target molecule, we relax the generated structures with xTB and compute their optimally-realignedd 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 (294% 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
Natural product ligand hopping- Numerous clipicall--aprrovedd 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-access accesssible small- molecule analogues of natural products that mimic their precise 3D interactions is preeminent task in scaffold/ligand hopping. To imitate this task, we select three complex natural products from CO- CONUT (Sorokina al. 2021), including two large macrocycles and a fused ring system with stereocenters. We then apply ShEPhERD (trained to learn P(z1. Zy Ix) on ShEPhERD-MOSES aq) to generate drug- drug-like molecules conditioned jointly on the ESP surface and pharmacophores of the lowest-energy conformer each natural product, again via inpainting. We emphasize that these natural products citt-f-diitrrbution compared to the drug-like molecules contained in ShEP.

---

Published as a conference paper at ICLR 2025
hERD HERD-MOSES-aG, Upon generating 2500 samples from the conditional distribution P(z1 2.)) for each natural product, ShEPhERD identifies mollcull 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, 2009). 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. (020)) (App.
Bioactive hit diversification. Whereas stucture-based drug design aims to design high-affinity lig- ands 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 anteraction-prseevvvii
scaffold hopping, often as a means to reduce toxicity, increase synthesizability, or evade patent restrictions. We simulate bioac- tive 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 2010 Eberhardt et al.) (02T) to dock the generated ligands to their respective proteins, treating the weak surrogate for bioactivity. To best imitate ligand- -bseed design, condition ShEPhERD on the ESP surface and pharmacophores of the lowest-enery 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( [ag ita), and dock the valid samples. Fig. 4 shows the distributions of Vina scores for ShEPhERD generrted analogues, compared against a docking screen of 10K random compounds from ShEPhERD-MOSES
Vina docking scores
Despite having no explicit knowledge of the protein targets, ShEPhERD enriches Vina scores in multiple cases. For the top- scoring generated ligands for 5mo4 and 7111, ShEPhERD generates substantial scaffold hops that yield 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, Willss et al. (2024) showed that merging (not linking) the fragments to form a bioisosteric ligand (which may not contain the exact fragment hits) can diversify ligand hypothe- ses 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 bind the antiviral target EV-D68 3C protease (Lithgo et al.. 2024] Wills al. 2024). We extract n2 = 27 pharmacophores by clustering common motifs and selecting interactions iden- tified by Fragalysis (Diamond Light Source 2024) (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 (n1 [50, 89)) from P(æ |æ3, 224) via inpainting. Fig. showw samples with SA < 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 n1 and n are significantly cut-of-distribution from SEEPhERD-MOSER aq (App. AA.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-ssmiilrrty 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-awart chemical design such as structure-base drug design and organocatalyst design.
target structures, as assessed via custom 3D
10

---

Published as a conference paper at ICLR 2025
REPRODUCIBILITY STATEMENT
Our main text and appendices provide all critical details necessary 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/ wtwpo//oolygouo///oooooooooo shepherd and http:s////ithuu.cooo////// Shepherd-scom score
ACKNOWLEDGMENTS
The authors would like thank Wenhao Gao for support with SynFormer. This research was sup- ported by the Office of Naval Research under grant number ONR 400014-21-1-2195. This material 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.
REFERENCES
Keir Adams and Connor W Coley. Equivariant shape-conditioned generation of 3d molecules for ligand-based drug design. arXiv preprint 10. 2022.
Dylan M Anstine and Olexandr Isayev. Generative models as an emerging paradigm in the chemical sciences. Journal of the American Chemical Society, 145(16):8736-875) 2023.
Mahendra Awale, Finton Sirockin, Nikolaus Stiefl, and Lean-Louis Reymond. Medicinal chemistry aware database gdbmedchem. Molecular informatics, 38(8-9):1900031 2019.
Christoph Bannwarth, Sebastian Ehlert, and Stefan Grimme. Gfn2-xtb an accurate and broadly parametrized self- consisennt tight-binding quantum chemical method with multipole electrostatics and Density-deeendeen dispersion contributions. Journal of chemical theory and computation, 15 (3): 1652-1111,,,,, 2019.
Francois Berenger and Koji Tsuda. 3d-sensitive encoding of pharmacophore features. Journal of Chemical Information and Modeling, 63(8):2360-239, 2023.
Helen M Berman, John Westbrook, Zukang Feng, Gary Gilliland, Talapady N Bhat, Helge Weissig, Πlya N Shindyalov, and Philip Bourne. The protein data bank. Nucleic acids research, 28(1): 235-242, 2000.
Camille Bilodeau, Wengong Jin, Tommi Jaakkola, Regina Barzilay, and Klavs F Jensen. Gener- ative models for molecular discovery: Recent advances and challenges, Wiley Interdisciplinary Reviews: Computational Molecular Science, 12(5):e 2022.
Caterina Bissantz, Bernd Kuhn, and Martin Stahl. A medicinal chemist's guide to molecular inter- 2010.
actions. Journal of medicinal chemistry,
Thomas Blaschke, Josep Arús-Pous, Hongming Chen, Christian Margreitter, Christian Tyrchan, Ola Engkvist, Kostas Papadopoulos, and Atanas Patronov. REINVENT 2.0: An AI Tool for de Novo Drug Design. Journal of Chemical Information and Modeling. 60(12):5918-59)) December 2020. ISSN 1549-9596, 1549-960X
Giovanni Bolcato, Esther Heid, and Jonas Boström. On the value of using 3d shape and electrostatic similarities in deep generative methods. Journal of Chemical Information and Modeling, 62(6): 1388-1398, 2022.
Sven Bühlmann and Lean-Louis Reymond. Chembl-likeness score and database gdbchembl. Fron- fiers in chemistry, 8:46, 2020.
Tim Cheeseright, Mark Mackey, Sally Rose, and Andy Vinter. Molecular field extrema descriptors of biological activity: definition and validation. Journal of chemical information and modeling, 46(2):665-676, 2006.
11

---

Published as a conference paper at ICLR 2025
Ziqi Chen, Bo Peng, Srinivasan Parthasarathy, and Xia Ning. Shape-conditioned 3d molecule gen- eration via equivariant diffusion models. arXiv preprint arXiv:2308. /1890, 2023.
Ann E Cleves, Stephen R Johnson, and Ajay N Jain. Electrostatic-fiel and surface-shape shape similarity for virtual screening and pose prediction. Journal of computer- aided molecular design, 33(10); 865-886, 2019.
Avelino Corma, Fernando Rey, Jordi Rius, Maria J Sabater, and Susana Valencia. Supramolecular self-assembled molecules as organic directing agent for synthesis of zeolites, Nature, 47006)))); 287-290, 2004.
W. Cowan-Jacob. Abll kinase (1334i d382n) in complex with asciminib and nilotinib. Protein Data Bank, 2016. PDB ID: 5mo4.
M.G. Deshmukh, J.A. Ippolito, E.A. Stone, W.L. Jorgensen, and K.S. Anderson. Crystal structure of the sars-cov-2(2019-ncoo) main protease in complex with compound 5. Protein Data Bank, 2020. PDB ID: 7111.
Diamond Light Source. Fragalysis, 2024. URL https://ffagalystttt http:/rrggalylss..iian.ooo.... www.diaadooo.oomommmm ac. uk/ Accessed September 15, 2024.
Steven I. Dixon, Alexander M Smondyrev, Eric Knoll, Shashidhar Rao, David E Shaw, and Richard A Friesner. Phase: new engine for pharmacophore perception, 3d qsar model develop- ment, and 3d database screening: 1. methodology and preliminary results. Journal of computer- aided molecular design, 20:647-671, 2006.
Yu Dong, Xiaodi Qiu, Neil Shaw, Yueyang Xu, Yuna Sun, Xuemei Li, Jun Li, and Zihe Rao. Molecu- lar basis for the inhibition of B-hydroxyacyl-acp dehydratase hadab complex from mycobacterium tuberculosis by flavonoid inhibitors. Protein & cell, 6(7):504 517,45171550 2015.
Jerome Eberhardt, Diogo Martos-Martins, Andreas Tillack, and Stefano Forli. Autodock vina 1.2. O: New docking methods, expanded force field, and python bindings. Journal of chemical information and modeling, 61(8)388113889) 2021.
Peter Ertl and Ansgar Schuffenhauer. Estimation of synthetic accessibility score of drug-like molecules based on molecular complexity and fragment contributions. Journal of cheminfor- matics, 1:1-11, 2009.
Alexander Fanourakis, Philip Docherty, Padon Chuentragool, and Robert Phipps. Recent devel-
opments in enantioselective transition metal catalysis featuring attractive noncovalent interactions between ligand and substrate. ACS catalysis, 10612)11077210771 2020.
Lorna R Fiedler, Kathryn Chapman, Min Xie, Evic Maifoshie, Micacla Jenkins, Pelin Arabacilar Golforoush, Mohamed Bellahcene, Michela Noseda, Dörte Faust, Ashley Jarvis, et al. Map4k4 inhibition promotes survival of human stem cell-derived cardiomyocytes and reduces infarct size 2019.
in vivo. Cell Stem Cell, 24(4):579-59)) 2019.
forlilab. Meeko, n.d. URL httppss//gituubb.cmffff////
Available at
https//giinbb.co/frrribb//eeee
Wenhao Gao, Tianfan Fu, Jimeng Sun, and Connor W. Coley. Sample Efficiency Matters: A Bench- mark for Practical Molecular Optimization. In Thirty-sixth Conference on Neural Information Processing Systems Datasets and Benchmarks Track, June 2022.
Wenhao Gao, Shitong Luo, and Connor W Coley. Generative artificial intelligence for navigating synthesizable chemical space. arXiv preprint arXiv:2410.034944, 2024.
Niklas WA Gebauer, Michael Gastegger, Stefaan SP Hessmann, Hlaus-Robert Müller, and Kristof T Schütt. Inverse design of 3d molecular structures with conditional generative neural networks. Nature communications, 13(1):933, 2022.
Mario Geiger and Tess Smidt. e3nn: Euclidean neural networks. arXiv preprint arXi::2207.09553, 2022.
12

---

Published as a conference paper at ICLR 2025
Rafael Jomeez-oomaarelli Jennifer N Wei, David Duvenaud, José Miguel Hernández-Lobato, Timothy D Hirzel,
Benjamín Baanchez-Lengllngg Dennis Sheberla, Jorge Ryan P Adams, and Alán Aspuru-Guzik. Automatic chemical design using data-driven contin- uous representation of molecules. ACS central science, 4(2):268-276, 2018.
AC Good. The calculation of molecular similarity: Alternative formulas, data manipulation and graphical display. Journal of Molecular Graphics, 10(3):144-1511 1992.
Andrew C Good, Edward E Hodgkin, and W Graham Richards. Utilization of gaussian functions for the rapid evaluation of molecular similarity. Journal of Chemical Information and Computer Sciences, 32(3):188 111-1188819, 1992.
J Andrew Grant and BT Pickup. A gaussian description of molecular shape. The Journal of Physical Chemistry, 99(1 1):3533-3333 3510, 1995.
J Andrew Grant, Maria A Gallardo, and Barry T Pickup. A fast method of molecular shape compari- son: A simple application of gaussian description of molecular shape. Journal of computational chemistry, (7(14):65331166, 1996.
Francesca Grisoni, Daniel Merk, Viviana Consonni, Jan A Hiss, Sara Giani Tagliabue, Roberto Todeschini, and Gisbert Schneider. Scaffold hopping from natural products to synthetic mimetics by holistic molecular similarity. Communications Chemistry, 1(1):44, 2018.
Jiaqi Guan, Wesley Wei Qian, Xingang Peng, Yufeng Su, Jian Peng, and Jianzhu Ma. 3d equiv- ariant diffusion for target-aware molecule generation and affinity prediction. arXiv preprint arXiv:2303.03543, 2023.
Charles Harris, Kieran Didi, Arian R Jamasb, Chaitanya K Joshi, Simon V Mathis, Pietro Lio, and Tom Blundell. Benchmarking generated poses: How rational is structure- based drug design with generative models? arXiv preprint arXiv:2308. 2.0.0:230..04113, 2023.
Jonathan Ho, Ajay Jain, and Pieter Abbeel. Denoising diffusion probabilistic models. Advances in neural information processing systems, 33:6840-6851 2020.
Edward E Hodgkin and W Graham Richards. Molecular similarity based on electrostatic potential and electric field. International Journal of Quantum Chemistry, 32($14):105-110, 1987.
Emiel Hoogeboom, Victor Gareia Satorras, Clément Vignac, and Max Welling. Equivariant diffu- sion for molecule generation in 3d. In International conference machine learning, pp. 8867- 8887. PMLR, 2022.
Kexin Huang, Tianfan Fu, Wenhao Gao. Yue Zhao, Yusuf Roohani, Jure Leskovec, Connor W Co- ley, Cao Xiao, Jimeng Sun, and Marinka Zitnik. Therapeuties data commons: Machine learning datasets and tasks for drug discovery and development. Proceedings of Neural Information Pro- cessing Systems, NeurlPS Datasets and Benchmarks, 2021.
Zhilin Huang, Ling Yang, Xiangxin Zhou, Zhilong Zhang, Wentao Zhang, Xiawu Zheng, Jie Chen, Yu Wang, CUI Bin, and Wenming Yang. Protein-ligand interaction prior for binding-aware 3d molecule diffusion models. In The Twelfth International Conference on Learning Representa- tions, 2024.
David J Huggins, Woody Sherman, and Bruce Tidor. Rational approaches to improving selectivity in drug design. Journal of medicinal chemistry, 55(4):1424-1444 2012.
Ilia Igashov, Hannes Stärk, Clément Vignac, Amne Schneuing, Victor Garcia Satorras, Pascal Frossard, Max Welling, Michael Bronstein, and Bruno Correia. Equivariant 3d-conditional diffu- sion model for molecular linker design. Nature Machine Intelligence, pp. 1-11, 2024.
Fergus Imric, Thomas E Hadfield, Anthony R Bradley, and Charlotte M Deane. Deep generative design with 3d pharmacophoric constraints. Chemical science, 12(44))4457711458, 2021.
Ross Irwin, Alessandro Tibo, Jon-Paul Janet, and Simon Olsson. Efficient 3d molecular generation with flow matching and scale optimal transport. arXiv preprint arXiv:2406.07266, 2024.
13

---

Published as a conference paper at ICLR 2025
V.-P. Jaakola, M.T. Griffith, M.A. Hanson, V. Cherezov, E.T.T. Chien, J.R. Lane, A.P. Ijzerman, and R.C. Stevens. The 2.6 crystal structure of human a2a adenosine receptor bound to zm241385. Protein Data Bank, 2008a. PDB ID: 3eml.
Veli-Pekka Jaakola, Mark T Griffith, Michael A Hanson, Vadim Cherezov, Ellen YT Chien, Robert Lane, Adriaan P Ijzerman, and Raymond C Stevens. The 2.6 angstrom crystal structure of human a2a adenosine receptor bound to an antagonist. Science, (22(5905)
2008b.
Victoria Jackson, Linda Jordan, Ryan N Burgin, Oliver JS McGaw, Calum W Muir, and Victor Cehan. Application of molecular-modeling scaffold-hopping, and bioisosteric approaches to the discovery of heterocyclic picolinamides. Journal of Agricudural and Food Chemistry, 70 136):11031-1104,, 2022.
Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Junetion tree variational autoencoder for molecular graph generation. In International conference on machine learning, pp. 2323-2332. PMLR, 2018.
Robert R Knowles and Eric N Jacobsen. Attractive noncovalent interactions in asymmetric catalysis: Links between enzymes and small molecule catalysts. Proceedings of the National Academy of Sciences, 107(48))20678 20685, 2010.
David Ryan Koes and Carlos Camacho. Pharmer: efficient and exact pharmacophore search. Journal of chemical information and modeling. 5166:1307-1314, 2011.
Alina Kutlushina, Aigul Khakimova, Timur Madzhidov, and Pavel Polishchuk. Ligand-based phar- macophore modeling using novel 3d pharmacophore signatures. Molecules, 23(12):3094, 2018.
Greg Landrum al. Rdkit: Open-source cheminformatics, 2006a.
Gregory A Landrum, Julie E Penzotti, and Santosh Putta. Feature-map vectors: a new class of informative descriptors for computational drug discovery. Journal of computer-aided molecular design, 20:751-776 2006b.
Sarah R Langdon, Peter Ertl, and Nathan Brown. Bioisosteric replacement and scaffold hopping in lead generation and optimization. Molecular informatics, 29(5):366-3855 2010.
Tuan Le, Julian Cremer, Djork-Aré Clevert, and Kristof T Schütt.. Latent-guided equivariant diffu- sion for controlled structure-based de novo ligand generation. In ICML '24 Workshop fo Life and Material Science: From Theory to Industry Applications, 2024.
Joongwon Lee, Wonho Zhung, and Woo Youn Kim. Ncidiff: Non-covalent interaction interaction-generativv diffusion model for improving reliability molecule generation inside protein pocket. arXiv preprint arXiv:2405. 16861, 2024.
Li, Dong, and Z.H. Rao. Crystal structure of (3r--hydroxyycyllacp dehydratase hadab hetero- dimer from mycobacterium tuberculosis complexed with 2' 2"444"-trihydroxxx -trihydroxychallooe.. Protein Data Bank, 2014. PDB 4rlu.
Yi-Lun Liao, Brandon Wood, Abhishek Das, and Tess Smidt. Equiformerv2: Improved equivariant transformer for scaling to higher-degree representations. arXiv preprint arXiv:2306, /2059, 2023.
Fang-Yt Yu Lin and Alexander D MacKerell Jr. Do halogen-hydrogen bond donor interactions dom- inate the favorable contribution halogens to ligand-protein binding? The Journal of Physical Chemistry B, 121(28):6813 1--26811,11 2017.
Jie Lin, Mingyuan Xu, and Hongming Chen. Diff-shape: A novel constrained diffusion model for shape based de novo drug design. ChemRxiv, 2024. doi: 10.26434/chemrxiv -024444----22-------11 km0h1.
Ryan M Lithgo, Charlie WE Tomlinson, Michael Fairhead, Max Winokan, Warren Thompson, Conor Wild, Jasmin Aschenbrenner. Blake Balcomb, Peter G Marples, Anu V Chandran, et al. Crystallographic fragment screen of enterovirus d68 3c protease and iterative design of lead-like compounds using structure-guided expansions. bioRxiv, pp. 2024-04, 2024.
14

---

Published as a conference paper at ICLR 2025
Andreas Lugmayr, Martin Danelljan, Andres Romero, Fisher Yu, Radu Timofte, and Luc Van Gool. Repaint: Inpainting using denoising diffusion probabilistic models. In Proceedings of the IEEE/CVF conference computer vision and pattern recognition, pp. 11461-11471, 2022.
Youzhi Luo and Shuiwang Ji. An autoregressive flow model for 3d molecular geometry generation from scratch. In International conference learning representations (ICLR), 2022.
Elaine C Meng, Thomas D Goddard, Eric F Pettersen, Greg Couch, Zach Pearson, John H Morris, and Thomas E Ferrin. UCSF ChimeraX: Tools for structure building and analysis. Protein Science, 32(1 I 1):422 2023.
Anthony J Metrano and Scott J Miller. Peptide-based catalysts reach the outer sphere through remote desymmetrization and atroposelectivity. Accounts of Chemical Research, 52(1):199-215, 2018.
B. Nagar, W. Bornmann, T. Schindler, B. Clarkson, and J. Kuriyan. Crystal structure of the c-abl kinase domain in complex with sti-571. Protein Data Bank, 2001. PDB ID: liep.
Bhushan Nagar, William G Bornmann, Patricia Pellicena, Thomas Schindler, Darren R Veach, W Todd Miller, Bayard Clarkson, and John Kuriyan. Crystal structures of the kinase domain of c-abl in complex with the small molecule inhibitors pd173955 and imatinib (sti-571). Cancer research, 62(15):4236-42433 2002.
Maruti Naik, Anandkumar Raichurkar, Balachandra S Bandodkar, Begur Varun, Shantika Bhat, Rajesh Kalkhambkar, Kannan Murugan, Rani Menon, Jyothi Bhat, Beena Paul, et al. Struc- ture guided lead generation for m. tuberculosis thymidylate kinase (mtb tmk): discovery of 3- cyanopyridone and 1, 6-naphthyridin-2-om as potent inhibitors. Journal of medicinal chemistry, 58(2):753-766, 2015.
Rebecca M Neeser, Mehmet Akdel, Daniel Kovtun, and Luca Naef. Reinforcement learning-drived linker design via fast attention-ban based point cloud alignment. arXiv preprint arXiv:2306.08 /8/66, 2023.
Alexander Quinn Nichol and Prafulla Dhariwal. Improved denoising diffusion probabilistic models. In International conference machine learning, pp. 8162-8171, PMLR, 2021.
OpenEye, Cadence Molecular Sciences. EON 3.0.0.0. http:://ww.eeessopee com 2024. Santa Fe, NM.
Tudor Oprea and Hans Matter. Integrating virtual screening lead discovery. Current opinion in chemical biology, 8(4):349-358, 2004.
Alaa MA Osman and Alya A Arabi. Average electron density: A quantitative tool for evaluating non-classical bioisosteres of amides. ACS omega, 9(1))13172-13182 2024.
Kostas Papadopoulos, Kathryn A Giblin, Jon Paul Janet, Atanas Patronov, and Ola Engkvist. De novo design with deep generative models based on 3d similarity scoring. Bioorganic & Medicinal Chemistry, 16308, 2021.
Xingang Peng, Shitong Luo. Jiaqi Guan, Qi Xie, Jian Peng, and Jianzhu Ma. Pocket2mol: Effi- cient molecular sampling based on 3d protein pockets. In International Conference on Machine Learning. pp.
PMLR, 2022.
Xingang Peng, Jiaqi Guan, Qiang Liu, and Jianzhu Ma. Moldiff: Addressing the atom-bond incon- sistency problem in 3d molecule diffusion generation. arXiv preprint 275.:2330.007508, 2023.
Peter Politzer, Jane S Murray, and Timothy Clark. Halogen bonding and other o-hole interactions: 2013.
A perspective. Physical Chemistry Chemical Physics,
Daniil Polykovskiy, Alexander Zhebrak, Benjamin Sanchez-Lengeling. Sergey Golovanov, Oktai Tatanov, Stanislav Belyaev, Rauf Kurbanov, Aleksey Artamonov, Vladimir Aladinskiy, Mark Veselov, et al. Molecular sets (moses): a benchmarking platform for molecular generation models.
Frontiers in pharmacology, 1 11:566644, 2020.
15

---

Published as a conference paper at ICLR 2025
Da-Hui Qu, ChaoCChn Wang. Qi-Wei Zhang, Xiang Ma, and He Tian. Photoresponsive host- functional systems. Chemical reviews, 115(15):5543-788) 2015.
Matthieu Raynal, Pablo Ballester, Anton Vidal-Ferran, and Piet WNM van Leeuwen. Supramolecu- catalysis. part I: non-covalent interactions as tool for building and modifying homogeneous catalysts. Chemicul Society Reviews, 43(5):1660-1777 2014.
J.A. Read, S. Hussein, H. Gingell, and J. Tucker. Mtb tmk in complex with compound 8, Protein Data Bank, 2014. PDB ID: 4unn.
James P Roney, Paul Maragakis, Peter Skopp, and David E Shaw. Generating realistic 3d molecules with an equivariant conditional likelihood model. 2022.
Lars Ruddigkeit, Ruud Van Deursen, Lorenz C Blum, and Jean-Louis Reymond. Enumeration of 166 billion organic small molecules in the chemical universe database gdb-17. Journal of chemical information and modeling, 52(1 0:2864-2877, 2012.
Thomas Rush, J Andrew Grant, Lidia Mosyak, and Anthony Nicholls, A shape-based 3-d scaffold hopping method and its application to a bacterial protein- protein interaction. Journal of medicinal chemistry, 48(5): 1489 1495, 2005.
Masami Sako, NOBUAKI YASUO, and Masakazu Sekijima. Diffint: A pharmacophore-awan diffusion model for srrueture-bssdd drug design with explicit hydrogen bond interaction guidance. ChemRxiv, 2024. doi: 1..26444cchemmxii-202---23
Marijn PA Sanders, Arménio JM Barbosa, Barbara Zarzycka, Gerry AF Nicolaes, Jan PG Klomp, Jacob De Vlieg, and Alberto Del Rio. Comparative analysis of pharmacophore screening tools. Journal of chemical information and modeling, 12(6):1607-1620, 2012.
Victor Gareia Satorras, Emiel Hoogeboom, and Max Welling- E (n) equivariant graph neural net- works. In International conference on machine learning, pp. 9323-9332. PMLR, 2021.
Gisbert Schneider, Petra Schneider, and Steffen Renner. Scaffold Scaffold-hopping: how far can you jump?
QSAR & Combinatorial Science, 25(12):1162-1177) 2006.
Arne Schneuing. Yuanqi Du, Charles Harris, Arian Jamasb, Ilia Igashov, Weitao Du, Tom Blun- dell, Pietro Lió, Carla Gomes, Max Welling, et al. ttructure-based drug design with equivariant diffusion models. arXiv preprint arXiv:22 13695, 2022.
Marwin HS Segler, Thierry Kogej, Christian Tyrchan, and Mark P Waller. Generating focused molecule libraries for drug discovery with recurrent neural networks. ACS central science, 4(1): 120-131, 2018
Miha Skalic, José Jiménez, Davide Sabbadin, and Gianni De Fabritüs. Shape-based generative modeling for de novo drug design. Journal of chemical information and modeling, 59(3): 1205- 1214, 2019.
Jascha Sohl-Dickstein, Eric Weiss, Niru Maheswaranathan, and Surya Ganguli. Deep unsupervised learning using nonequilibrium thermodynamics. In International conference on machine learn- ing, pp. 2256-2265. PMLR, 2015.
Maria Sorokina, Peter Merseburger, Kohulan Rajan, Mchmet Aziz Yirik, and Christoph Steinbeck. Coconut online: collection of open natural products database. Journal of Cheminformatics, 13 (1):2, 2021.
PW Sprague. Automated chemical hypothesis generation and database searching with catalyst. Perspectives in Drug Discovery and Design, 3, 1995.
Teague Sterling and John J. Irwin. Zine 15 ligand discovery for everyone. Journal of Chemical Information and Modeling, 55(11):2324-2333 2015. doi: (0.1022/acs.jiimm5bo00
Carmen Stoffelen and Jurriaan Huskens. Soft supramolecular nanoparticles by noncovalent and host-guest interactions. Small, 12(1):96-119, 2016.
16

---

Published as a conference paper at ICLR 2025
Jonatan Taminau, Gert Thijs, and Hans De Winter. Pharao: pharmacophore alignment and optimiza- tion. Journal of Molecular Graphics and Modelling, 27(2):16111) 2008.
F Dean Toste, Matthew Sigman, and Scott Miller. Pursuit of noncovalent interactions for strategic site-selective catalysis. Accounts of chemical research, 50(3):609-61 615, 2017.
Oleg Trott and Arthur Olson. Autodock vina: improving the speed and accuracy of docking with new scoring function, efficient optimization, and multithreading. Journal of computational chemistry, 31(2):55-461, 2010.
Mikko Vainio, J Santeri Puranen, and Mark S Johnson. Shaep: molecular overlay based on shape and electrostatic potential, 2009.
Clement Vignac, Igor Krawczuk, Antoine Siraudin, Bohan Wang, Volkan Cevher, and Pas- cal Frossard. Digress: Discrete denoising diffusion for graph generation. arXiv preprint arXiv:2209. 14734, 2022.
Clement Vignac, Nagham Osman, Laura Toni, and Pascal Frossard. Midi: Mixed graph and 3d denoising diffusion for molecule generation. In Joint European Conference on Machine Learning and Knowledge Discovery in Databases, pp. 560-576. Springer, 2023.
D. Wacker, G. Fenalti, M.A. Brown, Katritch, R. Abagyan, V. Cherezov, and R.C. Stevens. Crystal structure of the human beta2 adrenergic receptor in complex with the inverse agonist ici 118,551. Protein Data Bank, 2010a. PDB ID: 3ny8.
Daniel Wacker, Gustavo Fenalti, Monica A Brown, Vsevolod Katritch, Ruben Abagyan, Vadim Cherezov, and Raymond C Stevens. Conserved binding mode of human 32 adrenergic recep- tor inverse agonists and antagonist revealed by x-ray crystallography, Journal of the American Chemical Society, 132(33):1 1443 11445, 2010b.
Joel Wahl. Phesa: An open- soencs tool for pharmacophore-cnhaneedd shape alignment. Journal of Chemical Information and Modeling, 2024.
Mingyang Wang, Chang-Ya Yu Hsich, Jike Wang, Dong Wang, Gaoqi Weng, Chao Shen, Xiaojun Yao, Zhitong Bing, Honglin Li, Dongsheng Cao, et al. Relation: A deep generative model for structure-based de novo drug design. Journal of Medicinal Chemistry, 65(13):9478 9492, 2022.
Joseph I. Watson, David Juergens, Nathaniel R Bennett, Brian L. Trippe, Jason Yim, Helen E Eise- nach, Woody Ahern, Andrew J Borst, Robert J Ragotte, Lukas F Milles, al. De novo design of 2023.
protein structure and function with rfdiffusion. Nature, 62097976;:008--11
Camille G Wermuth. Similarity in drugs: reflections on analogue design. Drug Discovery Today, 11(7-8):348-33) 2006.
Steven E Wheeler, Trevor Seguin, Yanfei Guan, and Analise C Doney, Noncovalent interactions in organocatalysis and the prospect of computational catalyst design. Accounts of Chemical Re- search, 1061-1699,, 2016.
Stephanie Wills, Ruben Lanchez-Garcia Stephen D Roughley, Andy Merritt, Roderick E Hubbard, Frank von Delft, and Charlotte M Deane. Expanding the scope of a catalogue search to bioisos- teric fragment merges using graph database approach. bioRxiv, pp. 2024-08, 2024.
Andrew A Wylie, Joseph Schoepfer, Wolfgang Jahnke, Sandra W Cowan-Jacob, Alice Loo, Pascal Furet, Andreas I. Marzinzik, Xavier Pelle, Jerry Donovan, Wenjing Zhu, et al. The allosteric inhibitor ab1001 enables dual targeting ber-abl1. Nature, 543(7647):733-737 2017.
Weixin Xie, Jianhang Zhang, Qin Xie, Chaojun Gong, Youjun Xu, Luhua Lai, and Jianfeng Pei. Accelerating discovery of novel and bioactive ligands with pharmacophore-infonee informmedd generative models. arXiv preprint (rXiv:2401.01059 2024.
Maria Zavodszky, Anjali Rohatgi, Jeffrey R Van Voorst, Honggao Yan, and Leslie A Kuhn. Scoring ligand similarity in structure- based virtual screening. Journal of Molecular Recognition, 22(4):
280-292, 2009.
17

---

Published as a conference paper at ICLR 2025
Chun-Hui Zhang, Elizabeth A Stone, Maya Deshmukh, Joseph A Ippolito, Mohammad M Ghahre- manpour, Julian firado-Rives, Krasimir A Spasov, Shuo Zhang, Yuka Takeo, Shalley N Kudalkar, et al. Potent noncovalent inhibitors of the main protease of sars-cov-2 from molecular sculpting of the drug perampanel guided by free energy perturbation calculations. ACS central science, (3):467-4755 2021.
Kangyu Zheng, Yingzhou Lu, Zaixi Zhang, Zhongwei Wan, Yao Ma, Marinka Zitnik, and Tianfan Fu. Structure-based drug design benchmark: Do 3d methods really dominate? arXiv preprint arXiv:24066033403, 2024.
Qian-Yr Zhou, Jaesik Park, and Vladlen Koltun. Open3D: modern library for 3D data processing. arXiv:1801.09847, 2018.
Huimin Zhu, Renyi Zhou, Dongsheng Cao, Jing Tang, and Min Li. A Sharmaccophoregguice deep learning approach for bioactive molecular generation. Nature Communications, 14(1):6234, 2023.
Wonho Zhung, Hyeongwoo Kim, and Woo Youn Kim. 3d molecular generative framework for interaction-guided drug design. Nature Communications, 15(1):2688, 2024.
Yael Ziv, Brian Marsden, and Charlotte Deane. Molsnapper: Conditioning diffusion for structure based drug design. bioRxiv, pp. 2024-03, 2024.
18

---

Published as a conference paper at ICLR 2025
A
APPENDIX
CONTENTS
Introduction
Related Work
3
Methodology
3.1
Defining representations of molecules and their interaction profiles
3.2
Shape, electrostatic, and pharmacophore similarity scoring functions
3.3
Joint diffusion of molecules and their interaction profiles with ShEPhERD
Experiments
Conclusion
10
Appendix
19
A.1
Comparing ShEPhERD to REINVENT and virtual screening for natural product
ligand hopping
20
A.1.I
REINVENT convergence
21
A.2
Additional details on experiments
23
A.2.1
Unconditional joint generation CEEEPERRDGGDDD
23
A.2.2
Conditional generation on ShEPHERD-GDBI7
23
A.2.3
Natural product ligand hopping
23
A.2.4
Bioactive hit diversification)
24
A.2.5
Bioisosteric fragment merging
28
A.3
Defining and analyzing validity of generated 3D molecules)
30
Unconditional and conditional generation metrics
36
A.5
Calculating interaction profiles from 3D molecular structures
................................
38
A.5.I
Shapes/surfaces
38
A.5.2
Electrostatic potential surfaces
38
A.5.3
Pharmacophores
38
A.6
Parameterization of 3D similartyy scoring functions
39
A.6.I
Shape similarity scoring function
41
A.6.2
ESP surface similarity scoring function
42
A.6.3
Pharmacophore similarity scoring function
42
A.7
Additional details on model design, training protocols, and sampling
43
A.7.I
Further details on ShEPhERD's denoising modules
43
A.7.2
Training and sampling procedures
45
A.7.3
Model hyperparameters
50
19

---

Published as a conference paper at ICLR 2025
A.8
Symmetry breaking in unconditional generation
52
A.9
Characterizing the need for out-of-distribuui performance
53
A. 10 Additional examples generated molecules
55
A.11 Training and inference resources
58
A. 12 Additional Experiments and Comparisons to Related Work
59
A.1 Comparison to SQUID for shape conpe-cond generation
59
A. 12.2
Comparisons against tiructure-base drug design models that explicitly en- code the protein pocket
59
A. 12.3
Comparisons against inpainting with DiftSBDD
61
A. 12.4 Comparisons against SynFormer on natural product ligand hopping
61
A.1 COMPARING ShEPhERD TO REINVENT AND VIRTUAL SCREENING FOR NATURAL PRODUCT LIGAND HOPPING
REINVENT, state-of-the-at baseline for generative molecular design and optimization, applies reinforcement learning policy to iteratively update SMILES recursive neural network (RNN) with a provided reward function (Blaschke al.] 2020). We applied REINVENT the three natural product ligand hopping tasks defined in section 4 using a combination of our ESP and pharma- cophore 3D similarity scoring functions (section 3.2 as REINVENT's reward function. We fol-
lowed the REINVENT implementation in the Practical Molecular Optimization (PMO) benchmark (Gao et al., 2022). REINVENT was pretrained on the ZINC database (Sterling & Irwin 2015) and was deployed with a batch size of 64, 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 x 10-4. Any SMILES which failed during pharmacophore scoring were assigned score of 0.
Since REINVENT generates molecules in a ID SMILES representation, we generate up to 5 con- formers for each SMILES in order apply 3D similarity scoring functions as the reward func- tion. The procedure to compute the reward for a single generated SMILES is as follows: 1) embed 5 conformers with RDKit's EmbedMultipleeorr each conformer with MMFF94 for a maximum of 200 steps; 3) cluster the conformers with Butina clustering using an RMSD threshold of 0.1 A; 4) relax each remaining conformer with xTB in im- plicit water; 5) extract the ESP surface 23 with 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 con- formers; 8) add the ESP and pharmacophore similarity scores to obtain one combined score per conformer; and 9) take the maximum score across the conformers the reward.
function, which uses ETKDG; 2) optimize
= 400 and pharmacophore profile a4 of each
Fig. 5) compares the distributions ESP and pharmacophore similarity for samples obtained by (I) using ShEPhERD (trained ShEPhERD-MOSES-are to sample 2500 molecules from P(x1 Lys x44) ) via inpainting; (2) virtually screening (VS) 2500 random 3D molecules from ShEP- hERD- aq;; and (3) optimizing REINVENT with oracle/sampling budget of 10,000. Note that REINVENT was pretrained on ZINC, and ShEPhERD's training set (MOSES-ag) is a small sub- 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-MOSEL (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 ~500 samples from ShEPhERD against ~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.
20

---

Published as a conference paper at ICLR 2025
Natural product (I)
Natural product (2)
Natural product (3)
Pharmacophore Similaris
Pharmacophork Similarity
RDWOKT
Similarity
2:30
SPEPPERD 2.50
stem tRD (hx2.99)
9:20
0.33
0.20
01
Pharmacophh
8:35
0.00
830
0.02
0.03
0.4
25
0n
0.
2t ESP Surface Similarity
0.3
04
26
0.5
04
0.
07
ESP Surface Similarity
ESP
Surface Similarity
Natural product
ShEPhERD samples
REINVENT samples
a
(3)
eeaa
Figure 5: Distributions of 3D ESP and pharmacophore similarity to each of three natural product tar- gets for small-molecules sampled via ShEPHERD, REINVENT, or virtual screening (VS). Similarity scores 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 included. Also visualized examples of top-scoring samples generated by ShEPhERD and REINVENT.
We also emphasize that unlike REINVENT, ShEPhERD is not trained to directly optimize 3D sim- ilarity 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, 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 genetic algorithm which iteratively mutates and evolves the best-scoring samples by partial forward-noising and subsequent denoising.
A.1.1 REINVENT CONVERGENCE
We plot the average score of each REINVENT batch vs. the iteration number for each natural product target in Fig. 6 For the first natural product, REINVENT appears converged. For the second, REINVENT appears close to converged. For the third natural product, REINVENT could likely find better scoring molecules if trained for longer. We emphasize that could also run ShEPhERD for longer, too. Although REINVENT may not necessarily be converged for each natural product, we follow the PMO benchmark (Gao al. 2022) and limit the number of oracle calls to 10,000. In practice, even allowing for just 10,000 REINVENT samples required multiple days of computation per target because of the expensive conformer sampling, xTB optimization, and alignment with our scoring functions. We also note that in contrast to REINVENT, ShEPhERD is not specifically trained to optimize similarity for any natural product (or for any particular molecule). Fig. also plots SA score distributions for both REINVENT- and ThEPhERDD-geneaale molecules.
21

---

Published as a conference paper at ICLR 2025
Natural product (1)
Natural product (2)
Natural product (3)
0.70
0.66
Averege betch score
betch
0.65
0.62
0.58
0.55
50
100
150
50
100
150
50
100
150
Iteration
Iteration
Iteration
4.0
4.0
3.5
3.5
3.0
3.0
2.5
2.0
1.5
REINVENT
SNEPHERD
REINVENT
ShEPPERD
REINVENT
ShEPhERD
Figure 6: (Top) The average score (ESP+pharmacophore 3D similarity) across batch for each iter- ation of the optimization for REINVENT. Only natural product (1) seems to have fully converged. All optimizations were concluded 10,000 oracle calls. (Bottom) The distributions of SA score for REINVENT and SLEPhERD for molecules with SA score < 4.5.
22

---

Published as a conference paper at ICLR 2025
A.2 ADDITIONAL DETAILS ON EXPERIMENTS
A.2.1 UNCONDITIONAL JOINT GENERATION ON SHEPHERD-GDBI7
We define self-consistency as the 3D similarity between the ShEPhERD's generated interaction pro- file (i.e., Z2, 23. or 24) and the true interaction extracted from the generated molecule æ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 preserve the natively generated pose of 21 as much as possible. Since the P(x1. x2) and P(11.20 233) models employ n22, R2 = 75, all selfs-oon similarity evaluations also use 75 surface points when extracting surface and ESP interaction profiles from the (xTB- relaxed and RMSID-aligned) generated structures. Note that this differs from our other experiments involving conditional generation, which employ n22 fig 400 for similarity scoring.
For each valid ShEPhERD -eenerrtd molecule, the lower bound of self-consistency is calculated by (globally) aligning sampled molecule from ShEPhERD-GDBIT 17 the true interaction profile. Since surface sampling is stochastic (App. A.5), the self-similarity between two surfaces extracted from the same 31 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 expec- tation of the self-similarity between repeated extractions of the same profile from the same 21. We estimate these upper bounds by resampling the true profiles for each of the xT'B-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-GIBI
We use our 3D similarity scoring functions to evaluate the conditional ShEPhERD models. For each SeEPPhRR-generatte sample that 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 f22 fi3 = 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. 19). which should be considered when comparing any two similarity scores. As baseline, we also sample 20 molecules from the ShEPhERD-GIE 1DB1 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. by manually searching the COCONUT database pittps://cooe
na tura products s net/) / for natural product with complex structures that are meaningfully out of distribution compared to SEEPHERD-MOSESS 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 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 A.5
A.5
For each natural product, sampled 2500 structures from P(æ: z. xz) by inpainting, where (æg, æt) are the ESP and pharmacophore interaction profiles of the natural product. We use the version of ShEPhERD trained learn P(æ1, 33, 224) on ShEPhERD-MOSES-e specified based on the target æ;. and swept over 25 values of nı 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 ni. 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 33 for both structures at n23 = 400. We used ESP-based alignment instead of pharmacophore- capture the 3D shapes of the molecules, and because ESP- based alignments are more robust. Given the all-nnd 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.I) A.I).A.I)
We used na = 75,
alignment ESP surfaces better
23

---

Published as a conference paper at ICLR 2025
A.2.4 BIOACTIVE HIT DIVERSIFICATION
We used the same PDB co-crystal structures as used in (Zheng et al. 2024) and use the prepared
receptor files. . pdbbt AutoDock Vina files provided by Therapeutic Data Commons (TDC) (Huang
ct al. 2021). These targets liep (Nagar et al.) 2001] 2002). 3eml (Jaakola et al. 22008a b). 3ny8
[Wacker et al. 2010a b). 4rlu (Li et al., 2014; Dong et al. 2015). 4unn (Read et al.j 2014
Naik
et al.; 2015), 5mo4 (Cowan-Jacob) 2016; Wylie al.] 2017), and 7111 (Deshmukh al. 2020
Zhang et al. 2021). We use the corresponding co-crystal ligands: STI, ZMA, JRZ, HCC,
QZZ,
Nilotinib (NIL), and XFI.
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. 2000) to also extract the interaction profiles and obtain the "PDB pose" conformation. We added hydrogens to these PDB poses using RDKit's 's AddHs function with addCoo TrreTrue, 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 single point xTB calculation in implicit water.
For the two poses of each experimental ligand, we generate 500 samples from P(r) (( an)) by inpainting, where (arz, mz) are the ESP and pharmacophore interaction profiles of the ligand. We use the ShEPhERD model that trained to Iearn P(x1 33 x2) on SHEPhERD-MOSES-
In similar manner to the natural products, we used n3 = 75 and specified n, based on the target vi. We swept over 25 values of 71 in the range [ni 12, n + 12] where ni is the total number of atoms in the PDB ligand. We generate 20 structures for each choice of n1. 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, implement adapted version of the Vina smi les 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 n.d.); 4) dock the molecule with Autodock Vina v1.2.5 (Trott & Olson 2010 Eberhardt et al.) 2021 using& random seed of 987654321 and exhaustiveness of 32; and 5) obtain a Vina score in kcal/mol, which we use as (poor) proxy for binding affinity.
Fig. showcases examples for the 5mo4 and 7111 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- scoring poses. We then visualize the top scoring pose of three top-scoring ShEPhEERD-eeneranee ligands, overlaid on the best docking pose of the PDB ligand. Although the top-scoring docked pose of the ThEPhERD-generatee 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 CHEEPERD--genraatee ligands that closely aligns with the top pose of the PDB ligand.
Fig. shows the distributions of the top- 10 Vina docking scores for ShE:hERRD-generaae ligands compared to the top-10 docking scores obtained by virtually screening 10K molecules from ShEP- hERD-MOSES-we 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 's bound pose. Note that conditioning on the PDB ligand's s bound pose indirectly leaks in- formation 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 MMEPhERD--ennraatee conditioning the PDB ligand's bound pose, the Vina docking score of the PDB ligand is rela- tively 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.
ligands when
To investigate the potential of using ShEPhERD to optimize docking score, we also condition ShEP. hERD the lowest energy conformer and docked pose of the best scoring compound from the 10K
24

---

Published as a conference paper at ICLR 2025
5mo4
PDB docked pose
Top-10 dacked PDB- ligand pose ensemble
scoring poses of ShEPLERD samples
7111
PDB pose best docked pose
Top-10 docked PDB- ligand pase
scoring ShEPhERD samples
Figure 7: (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 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 Fig. (green-grey), overlaid on the top-scoring re-docked pose for the experimental ligand (light blue). The top and bottom show results for 5mo4's s and 7111 's experimental co-crystal ligands, respectively. Although the top-1 docked pose for the ShEEPERRD- of the experimental ligand, the ensembles of their docked poses do closely match, indicating that the experimental ligands and the ShEPhERD-diversiffe 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.
ligands may not perfectly align with the top scoring pose
docking screen. Results are shown Fig. 9 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 experi- mental PDB ligand, shifting the mean of the top- 10 scores from -9.96 to 11.70 kcal/mol (Fig. [) This PDB ligand scores poorly (the worst among all PDB ligands). Another interesting phenomenon that neither the 10K screened molecules nor ShEPhERD -eenerated molecules (conditioned on the best from the screen) manage to surpass the docking score of 7111's experimental PDB ligand. However, ShEPhERD does improve this docking score when directly conditioning on the 7111 PDB ligand (Fig. 8 highlighting ShEPhERD's flexibility in structure-baat based drug design.
25

---

Published as a conference paper at ICLR 2025
8.5
-9.0
Vina Score (kcal/mol)
100.0
100.5
11..5
SibrfeeD - bound pose Snffheri loweet
-12.0
Igand
3eml
3ny8
4riu
5mo4
7i11
Target
Figure 8: The distributions of the top- 10 scoring ShEPhERD SreeereRD-generaeed 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 valid SMILES string from the generated 3D conformation.
All compounds
Cemerated coeking sererrated- eners 10k docking screen Redocied FCB ligand Cur Carditioning matocule
Vina Score (kcalimol)
-2
-10
-12
liep
3eml
3ryy
4umn
Smoo4
7111
Protein Target
Top-10
-8
Vina Score (kcalmmoll
-10
-11
Geraneed loved 10k docking
-13
Redocked P36 igand Condittoring melecute
liep 3eml
3ry8
4rlu
Aunn
Smo4
7111
Protein Target
Figure 9: Docking results for ShEPhERD-generaee molecules when conditioning on the best scor- ing 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 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- scoring compounds.
26

---

Published as a conference paper at ICLR 2025
5mo4
7111
Figure 10: The individually visualized docking poses for each sampled molecule for PDB targets 5mo4 and 7111 from Fig.
27

---

Published as a conference paper at ICLR 2025
A.2.5 BIOISOSTERIC FRAGMENT MERGING
We demonstrate ShEPhERD's native capacity for bioisosteric fragment merging on the antiviral tar- get EV-D68 3C protease. Starting from the 25 fragments hits in the catalytic pocket identified by Wills et al. (2024). 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 2024 and ChimeraX (Meng et al. 2023). Specifically, we removed fragments that did not have common interactions, that significantly deviated from the fragment ensemble's ss aggregate shape, added extra infor- mation. The final list of fragments were: x0147 _OA, x1071 1071_0A, x1083 0A, x1084_0A, xl 14, x1163 OA, x1498 OA,, x1537 0A, x1594 OA, x1919 _0A, x2021 _OA, x2099 0A, and x2135_0A 0A.
For each fragment, we added hydrogens using RDKit's AddHs function with SedCoords-Trus We then generated an aggregate shape representation by sampling from 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 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 threshold distance. These threshold distances for hydrophobes. These were somewhat arbitrary cutoffs that yielded reasonable pharmacophore hypotheses upon visual inspection. Table lists the final 27 selected pharmacophores.
1.0 Á for HBAs, 2.0 Á for HBDs, 1.5 Á for aromatic groups, and 2.0 Á
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 HBA
Coordinates (P)
Vector (V)
(-7.5155, -3.7035 -..065) -9.0635)
(0.1900, 0.7041, -0.6842)
(-6.7175,
(0.8205, -0.2287, -0.5240)
-,.3300-7.33
(0.7719, 0.4600, -0.4388)
(-6.5197, 43333,
(0.5106, 0.7063, -0.4905)
(-6.9450, -.9966,,5.,
(-0.1783, -0.4378, -0.8812)
(33.7780.77.217,, 2.2860)
60.7287.0062777 -0.2755)
-..7250, -0.2220
-0.631,, 0.6316)
(-1.9580, -8.8820, 3.3910)
(0.7956, -0.5564, -0.2396)
(-5.0980, -1.5370
(0.4851, 0.6853, -0.5432)
15.8830,-111999 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.217, 2.2860)
(0.7506, 0.2245, 0.6214)
(-1.9580, -8.880,,,
(-0.2298, 0.1607, 0.9599)
-.00000,-70774
(-0.3293, 0.3106, -0.8917)
Aromatic
(-7.6710, -.54775,-8.2231
(-0.5639, -0.4342, -0.7025)
(-11.6032, 2,006,,,
(0.3775, 0.9134, -5.15
-.000000-..09115-11000000
(0.4302, -0.8132, -0.3919)
(-9.5312, 4.3847,
(-0.5073, -0.3479, -0.7884)
Hydrophobe
(-7.6710, 4.5475.-8..231)
(0.0, 0.0, 0.0)
(-11.6032, -2.0268, -6.4003)
(0.0, 0.0, 0.0)
(-4.6030, -.0500
(0.0, 0.0, 0.0)
(-9.5312, -4.0044,,77
(0.0, 0.0, 0.0)
(-5.6776, -0.4630,
(0.0, 0.0, 0.0)
(-6.9299, -..430
(0.0, 0.0, 0.0)
(-6.7990, -455530, -3.9980)
(0.0, 0.0, 0.0)
(-7.9010, 0.1820, -0.7890)
(0.0, 0.0, 0.0)
28

---

Published as a conference paper at ICLR 2025
Using the extracted ESP (æx) and pharmacophore (xt) interaction profiles, used the version of ShEPhERD trained to learn P(11, 33. T4) on ShEPhERD BOEEERRD-OOES-AR aq to generate 1000 structures from p(x, [2, Ø44 via inpainting. We used ng 75, specified n4 27, and swept over 40 values of n1 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 the natural product by optimizing our ESP similarity scoring function after resampling 23 for both structures at n2 400. We used ESP-based alignment instead of pharmacophore-bais 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. [T] visualizes four ShEPhERDD-geneatte samples with neutral charge and SA< 4.0 that were in the top-10 by combined ESP and pharmacophore similarity scores.
alignment ESP surfaces better
Figure 11: Examples of four top-scoring ShEPhERD-gennratte ligands obtained via bioisosteric fragment merging. Shown are the xTB-relaxed structures of the generated 3D molecules after being realigned the fragments' ESP and pharmacophore interaction profiles. Aligument is performed by optimizing ESP similarity. The ligands are overlaid on the target interaction profiles. The target ESP surface has been upsampled for visualization.
29

---

Published as a conference paper at ICLR 2025
A.3 DEFINING AND ANALYZING VALIDITY OF GENERATED 3D MOLECULES
Table 2: Global validity reported percentages for 1000 unconditional and 2000 conditional samples from various ShEPhERD models trained on ShEPhERD-GDBIT, 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 samples where the atom con- nectivity remains the same given that the molecules are valid both pre- and post- oTtttiiimi and "Valid neutral" is the percentage of valid samples that have molecular charge = 0.
Unconditional Generation
Samples
P(æ (2)
P zzg)
P(æ Z44
Validity (%)
96.2
92.7
73.7
Validity xTB (%)
96.3
93.0
73.7
Graph consis. (%)
99.1
99.5
98.4
Valid neutral (%)
84.1
93.2
89.7
Conditional Generation
Samples
P(æ (22)
P(æ |2g)
P(æ 224))
Validity (%)
96.0
91.9
80.7
Validity xTB (%)
96.2
92.9
79.9
Graph consis. (%)
97.4
92.8
93.8
Valid neutral
(%)
84.2
82.5
72.6
Defining validity. We define validity as the fraction of molecules that can converted from the atomic point cloud (a, C) (i.e., element types and coordinates to RDKit Mol object via RDKit's in- ternal 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 sim- ple, convenient, and quite reliable when (1) explicit hydrogens are included in the 3D structure, (2) the 3D molecular geometry of high quality, and (3) the overall charge of the molecule known.
In particular, RDKit is used to (1) generate a Mol object from (a, C)) formatted as an xyz file; (2) determine bonds via rule-based algorithm based on interatomic distances; (3) sanitize the molecule; (4) ensure that there are 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 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 the single point calculation. The sampled æ 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 (a, C') (formatted as an .yzz file) and repeat the validity checks. the natively generated structure was determined invalid pre-optimization, the xTB optimization is performed with 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.
AAthoough ShEPhERD generates formal charges and the bonds B, ignore these outputs when extract- ing molecule from T1. We argue that validity of 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 I and B primarily as way to help regularize the diffusion processes over a and C. Moreover, by defining validity only with respect to the element types and coordinates, our validity checks may he applied to other 3D molecular generative models which do not generate the 2D molecular graph.
30

---

Published as a conference paper at ICLR 2025
Analyzing validity. Metries of validity for unconditional and conditional across all samples from ShEPhERD-GDBIT are presented in T There is a trend of relatively high validity (2 96%) for the joint generation of 2 and 22. slightly lower validity (91-93%) for the joint generation of Φ and 23,, and a significant decrease for the joint generation of æ1 and 244- We hypothesize that the lower models is hecause of the discrete nature of the pharmacophore interaction profile relative to the surface or ESP-attributed surface. Notably, validity improves in the pharmacophore-conditoen
validity of the pharmacophore-basee
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/substrnetures. Specifying the pharmacophores based on a target profile also avoids relying on ShEPhERD's ' jointly-generated 24 (which may have errors) to guide the denoising of æ+. One way to improve joint pharmacophore generation may be to introduce a noise schedule that denoises the pharmacophores earlier than (a, 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 confor- mations 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 Iearned to generated charged substructures, but hasn't 't been trained long enough to learn to neutralize the over- all molecule every time. Encouragingly, there is a higher rate of generating neutral species for the P(æ1, aDss and P(æ1 X24) models relative to the P(z1.22) 222 model, which perhaps implies that us- ing 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(æ (æa) and P(æi (a:4) models perform worse than P(æ) x22) 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 (c.g., carboxylate anions yield negative ESP surfaces), and hence introduces charged groups more often while failing neutralize the overall molecule.
Fig. 12 details how validity is impacted by the selection of the number of atoms (nz) for uncondi- tionally generated samples (top row) and conditionally generated samples (bottom row) from models trained on ShEPhERD-GDDE
In both unconditional and conditional settings, versions of ShEP- hERD that jointly model shape or electrostatics attain high validity rates when the number of atoms is in-distribution of the training set (Fig. 23). On the other hand, the pharmacophore models gen- erally 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 goinele-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., n1 exceeding the bounds of the training datasets), albeit at lower validity rates.
Table 3 and Figures 13 14 show the strain energy and heavy-atom RMSD of ShEPhERD-generaate 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-GDE 17 (mean strain energies <10 kcal/mol). Indeed, these strain energies are typically lower than the strain energies of molecules with MMFF94-optimizedd 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- -gnerated samples for models trained on ShPhERD-GDBI7 have average heavy- atom RMSDs below 0.08 A. We note that slightly distorted atomic positions can vastly increase strain energy, which makes strain energy an unforgiving (but still useful) measure of absolute geo- metric quality. We notice a similar trend to validity where samples from P(æ1, model perform the best in terms of strain energy while samples from P( 24) model perform the worst, though the RMSD's 's are quite low across all the unconditional models.
The models trained on ShEPhERD-MOSESS
perform worse in terms of both the strain energy and RMSD metrics (Table 3) though the trends with respect to n, are similar (Fig. [4). We note that
31

---

Published as a conference paper at ICLR 2025
the molecules in ShEPhERD-MOSEN aq are generally larger than those in ShEPhERD-GDE For instance, the largest molecules in MOSES have non-hydrogen atoms, whereas the largest molecules in GDB17 have 17 non-hydrogen atoms. SMEPhERD-MOSES, training examples than SHEPHERD-GDBIT we expect the performance models trained on ShEP- hERID-MOSES-ag will improve with further training larger drug-like datasets.
also contains fewer
P X2)
P(X1 X3)
P(X1, X4)
100
Validity (%)
40
Before aptimization
Atter xTB optimization
20
60
40 Number of Atoms
40
Number of Atoms
Number of Acoms
P(X₂ (x2)
P(x1 (2)
P(Xy (x4]
80
Validity
60
20
Before XTB optimization
After XTB optimization
20
40
60
20 Number Atoms
40 60
20
40 60
Number of Atoms
Number of Acoms
Figure 12: Validity as function of the number of generated atoms (n1) for 3D molecules æ1 gen- erated by various ShEPhERD models trained on ShEPHERD-GDBIT 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.
Table 3: The overall mean and standard deviation of strain energies and RMSDs for unconditional samples from versions of ShEPhERD trained on ShEPhERD-GiDB17 and ShEPhERD-MOSES (Fig. 30). We compute the metrics from the valid samples amongst 1000 total samples from each model. The column "MMFF94" is baseline that subsamples 1000 molecules from each respective dataset, generates
MMFF94-Ieve 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-GDBI 17, and in implicit water for
ShEPHERD-GDBI7
P(æ,, 222
P(a1,
P 1 244
MMFF94
Strain Energy
2.81
5.80
5.87
7.34
(kcal/mol)
+5.32
+73.17
+47.36
+9920
RMSD
0.084
0.081
0.089
0.193
(A)
+0.133
#0.144
+0.155
#0.184
ShEPHERD-MOSES-an
P( 1..1.8g.3.. 344
MMFF94
Strain Energy
17.43
7.99
(kcal/mol)
+124.02
13.26
RMSD
0.273
0.337
(A)
+0.292
+0.274
32

---

Published as a conference paper at ICLR 2025
P(X1, X2)
100
Strain energy (kcal/mol)
10a
13000
100 0.00
RMSD (A)
1am0
01000
1000
10-1
100
10-1
10-2
P(x1, X3)
100
0OOO
Strain energy
400 Φ
(kcal/moll)
102
RMSD (A)
- CBD
1000
1000 0
10-1
420000
40000
100
10-2
P(X1, Xa)
Strain energy
8
100
10
(kcal/mol)
102
RMSD (Å)
10
100
10
8.58
64
Number of atoms
Number of atoms
Figure 13: Strain energy and heavy-atom RMSD for samples unconditionally generated by various models trained ShEPhERD-GDBI7 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
P(X1, X3 Xa)
poo
Strain energy
10a
100
(/cll/mol)
RMSD (A)
102
DODO
1000
400
10-1
10
100
10-2
Number of atoms
Number of atoms
Figure 14: Strain energy and heavy-atom RMSD for valid samples (amongst 1000 total samples) unconditionally generated by SłEPhERD trained to learn P 2' 2.) on ShEPhERD-MOSES- aq. The box plots are binned by the number of atoms. Though the strain energies can be high, the RMSDs quite low. Summary statistics are reported in Table 3
33

---

Published as a conference paper at ICLR 2025
NP
NP 2
NP 3
100
80
60
Validity
50
70
80
40
50
70
80
40
80
Number of Atoms
Number of Atoms
Number of Atoms
1iep (LE)
3eml (LE)
3ny8 (LE)
100
80
Validity (%)
60
40
50
50
60
50
Number of Atoms
Number of Atoms
Number of Atoms
4rlu (LE)
4unn (LE)
5mo4 (LE)
100
Validity
50
20
30
60
70 Number of Atoms
Number of Atoms
Number of Atoms
7111 (LE)
Bioisoteric FM
100
Validity (8)
Legend
Reference
optimization
After xTB optimizaaiinn
SMILES extractable
40
50
60
80
90
Number of Atoms
Number of Atoms
Figure 15: Validity as function of the number generated atoms for ligands generated with ShEP- hERD 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 æ; in implicit water with xTB. For the bioactive hit di- versification experiments, we only check the ability to extract SMILES strings from the generated 2h xTooooTiiipttitiinnn There no clear patterns in the validity of the molecules generated by ShEPLERD 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. 22 particularly for the NP and FM experiments. Validity rates for each experiment are summarized in Table
34

---

Published as a conference paper at ICLR 2025
Table 4: Global validity metries for molecules generated by ShEPhERD for the natural product
ligand hopping, bioactive hit diversification, and bioisosteric fragment merging experiments.
Natural product ligand hopping
NP 1
NP 2
NP 3
Validity (%)
65.6
56.5
55.6
Validity xTB (%)
67.6
62.4
56.0
Graph consis. (%)
78.5
88.6
68.2
Valid neutral (%)
50.5
72.3
47.0
Bioactive hit diversification
Tiep
3eml
3ny8
4rlu
4unn
5mo4
7111
Validity (%)
32.7
73.3
52.3
88.5
66.5
71.5
53.5
Valid neutral
47.1
60.0
75.7
78.8
55.2
71.4
62.6
Bioisosteric fragment merging
Bioisosteric FM
Validity (%)
45.0
Validity xTB (%)
37.5
Graph consis. (%)
76.6
Valid neutral (%)
46.0
35

---

Published as a conference paper at ICLR 2025
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 the fraction of valid molecules that 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-eeneratio experiments) are reported in ShEPhERD attains nearly 100% novelty and uniqueness for all cases.
Table 5
Fig. 16 also reports specific property distributions (SA Score, logP. QED, Fsp3) of unconditionally- generated molecules from versions of ShEPhERD trained on ShEPhERD-GDDE
or ShEPhERD- MOSES-aq. Property distributions for generated molecules are compared against the distributions from random sampling of molecules from the corresponding datasets.
Table 5: Novelty and uniqueness metrics across different versions of ShEPhERD trained either ShEPhERD-GDBI7 or ShEPHERD-MOSESS-R The first six rows report metries for ShEPhERD models when trained on ShEPhERD-GIDE and applied to either unconditional (P(æ1. ax)) or con- ditional (P(x (x:)) generation. The remaining rows report metrics for ShEPhERD when trained to 1 33. x.1) indicates this model when applied to unconditionally generate 1000 samples, whereas P(x [ag 224 indicates this model when applied to conditional generation in the natural product (NP), bioactive hit diversification, and bioisosteric fragment merging (FM) experiments.
learn 1 33, 224) ShEPhERD-MOSES.
Samples
Novelty
Uniqueness
(%)
(%)
ShEPhERD-GDBI7
P(æ1,
96.4
99.5
P(x1, 233
96.7
99.6
P(æ1, 344
96.0
99.9
P( 222
99.4
97.9
P(æ | (C3
99.3
97.6
P(æ (24)
99.0
96.0
SHEPHERD-MOSES
P( 323 304)
100.0
100.0
NP 1
P x1
D33 4AA)
100.0
99.7
NP 2
P
D33 (0
100.0
100.0
3P
P 2
233 07.4
100.0
100.0
Tiep (lowest energy)
Pl x
323 (74
100.0
100.0
liep (pose)
P 22
033 (F4)
100.0
100.0
3eml (lowest energy)
P(æ x
033
374)
100.0
100.0
3eml (pose)
P(x x
E3g
244)
100.0
100.0
3ny8 (lowest energy)
P(x
Ø33
100.0
100.0
3ny8 (pose)
P(x x
Ø33 324)
100.0
100.0
4rlu (lowest energy)
P x1
23, 424)
99.5
100.0
4rlu (pose)
P x
E3, ..a)
99.8
100.0
4unn (lowest energy)
P x
033 274)
100.0
100.0
4unn (pose)
P(æ
23, (7.4)
100.0
100.0
5mo4 (lowest energy)
P
Ø3,
100.0
100.0
5mo4 (pose)
P(x
100.0
100.0
7111 (lowest energy)
P x
Ø3g 444)
100.0
100.0
7111 (pose)
P x1
E3, (.1)
100.0
100.0
Bioisosteric FM
P(x
033, 27,))
100.0
100.0
36

---

Published as a conference paper at ICLR 2025
ShEPHERD -GDB17
38
loge
8.4
loge
0.5-
3.4
6.3
3.
3.0
Figure 16: Distributions of molecular properties for uncondittionally-gneratdd molecules compared samples from the respective datasets. SA Score is the synthetic accessibility score, logP the octanol-water partition coefficient, QED is measure of drug-likeness, and Fsp3 is the fraction of sp3 carbons in a molecule. (Top) The distributions relevant to ShEPhERD-GDER show the property distributions for the valid molecules amongst 1000 SEEEPhERD-generann samples from P(æ1, 222), ((x) 223), and P(x1, X24). "Dataset Samples" refers to 1000 randomly sampled molecules from the ShEPhERD-GDI GDB17 dataset. (Bottom) The same distributions relevant to ShEP- hERD-MOSES-ag The plots show the property distributions for the valid molecules amongst 1000 ShEPhERD- samples from ( 23, 2:4). "Dataset Samples" refers to 1000 randomly sampled molecules from the ShEPHERD-MOSEN -aq dataset.
The plots
37

---

Published as a conference paper at ICLR 2025
A.5 CALCULATING INTERACTION PROFILES FROM 3D MOLECULAR STRUCTURES
A.5.1 SHAPES/SURFACES
We represent the shape of molecule as point cloud sampled from the solvent accessible surface
with probe radius "probe- Here, we describe how extract 22 S E 1 x3 from æ = (a, C) with E _
First, we sample points on the surfaces of nı unit spheres using Fibonacci sampling. The number of sampled points for each sphere is 25(ra (araaje//L77 /1.72 where ra(k] the van der Waals (vdW) radius of the kth atom in A. The points are then scaled to have norm [a(k] Tprobe and the sampled spheres are translated to their respective atomic coordinates in C. Points that fall within ra[A]] t of any other atom are filtered out. Next, we use Open3D (Zhou et al. 2018) to generate a Triangl aeesh eMesh using the ball pivoting algorithm with ball radius of 1.2. Finally, we sample surface point cloud with approximately evenly spaced points using Open3D's s ammpleppontt s.pontssppoii. sapoisson_dis.
method. This final step is stochastic. We used Ppoobe 0.6 for the modeling in this work, but tuning of similarity scoring function parameters (App. A.6 originally used pprobe 1.2.
A.5.2 ELECTROSTATIC POTENTIAL. SURFACES
We obtain the electrostatic potential (ESP) surface by computing the Coulombic potential at each point surface point cloud. Here, describe how to extract 33 (S,, v) from æ1 (a, C)..
We obtain S3 E x3 using the same procedure as for S2. We then obtain the partial charges q E Th1 from xTB by running single point caleulation ø1- Through internal testing, found that using partial charges obtained from the cheaper MMFF94 to be of insufficient quality and detrimental 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:
the performance of both ShEPhERD and the ESP similarity scoring function.
ng
[[k]
v u[k]
4T
[k
r +bll
where 6o is the vacuum permittivity with units (eV A)-
A.5.3 PHARMACOPHORESS
Pharmacophores (p. P. 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 क- stacking or cation-n interactions. We position aromatic pharmacophore the centroid of aro- matic ring. Two antiparallel vectors point in opposite directions, orthogonal to the plane of the aromatic ring.
Halogen-carbon bond halogen (F, CL. Br, I) that is connected to a carbon can form halogen bonds that can display HBA/HBD-like propertics. The vector points from the halogen and in the aati-parallel direction the neighboring carbon.
The non-directional pharmacophores are:
Hydrophobe- sulfur atom, halogen atom, group of atoms that are hydrophobic. We cluster groups of hydrophobic atoms such that distinct hydrophobes are > 2 A apart. Note that aromatic rings are also hydrophobes.
38

---

Published as a conference paper at ICLR 2025
Anion negatively charged atom the centroid of common anionic group.
Cation 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 2011 Landrum
et al. 2006a Kutlushina et al.] 2018) 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 hydrophohes within 2 Á are clustered using Berenger & Tsuda (2023) 's protocol; and 3) sepa- rate pharmacophore identities for halogens were made due their unique interaction geometry (c.g., sigma holes) that warranted a distinction from HBA/HBD pharmacophores (Lin & MacKerell Jr 2017 Politzer et al. 2013). Most of the pharmacophore extraction utilizes altered code from RD- Kit (Features Feat MettttlssAA 1l and Features. ShowF eats modules) (Landrum et al.) 2006b). 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 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.
P are either at the location of the identified atom or at the centroid of 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.
The placement of p
17 Each pharmacophore is associated with vector V and either be a unit vector () (V[[ 1) For directional pharmacophores or the zeTo vector () [V[] = 0) for directionless pharmacophores. For directional pharmacophores, a unit vector the direction of potential interaction (e.g., in the direction of hydrogen for an HBD or a lone pair for an HBA) is extracted. is common for multiple pharmacophores to be located the same position. For example, aromatic rings and halogen-carbol bonds are also hydrophobes (Fig. 17). and hydroxyl groups are both HBDs and HBAs.
A.6 PARAMETERIZAATION OF 3D SIMILARITY SCORING FUNCTIONS
Recall from Sec. that for two point clouds Q. and Q z we model each point T&e as an isotropic
Gaussian in R33 Grant & Pickup 1995 Grant et al. 1996). We compute the overlap between them OA.B with the first order Gaussian overlap. The Tanimoto similarity function sim' (Q Q 8) e
[0, is used define a similarity of point clouds based on their overlaps:
OA.B
Σ
Σ BEQ
Wo,b
exp
Ira - r r6613)
20
a
A
OA,B
sim" QB) O + O OA,B
Here, a parameterizes the Gaussian width and a. is a weighting factor that is tuned for each representation. Note that no, need not equal nis+ generally.
Since 3D similarities are sensitive to SE(3) transformations of QA with respect QB, their optimal alignment can be defined as:
sim Q QB) max sim' (RQ) + QB)
R,t
where R E SO(3) and E 733). The optimization landscape is non-convex and noisy. For global optimizations (e.g.. aligning random 3D molecules from the dataset to target natural product), we optimize the alignment of QA with respect to QB with N different initializations of (R, t). We found that N 50 is sufficient for convergence. For local optimizations (e.g., aligning 3D molecule generated by ShEPhERD to 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, is set that the center of mass (COM) of QA is aligned with the COM of Q.. One of the initializations uses the identity R = I (to keep the initial orientation). We also use alignments to four orientations relative to the principal moments
39

---

Zn binder Published as a conference paper at ICLR 2025
hydrophobe clustering Mutiple vectors
Hydrophobe clustering Averaged vectors
Pyridinium
Hydrogen sulfate
Hydroxymate
Examples of each type of pharmacophore. (Top) 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 struc- ture highlights how multiple vectors are used for HBAs when they have multiple lone pairs, while the right structure shows these vectors averaged into single direction (not shown here but the same applies for HBDs). Additionally, the clustering scheme from (Berenger Tsuda, (023) is used to reduce the number of hydrophobes. This molecule an slightly altered example from httpp::/ggregaandre
Figure 17:
github io/ wawii-kkioooooo
posts/2023 02 24 using- feature maps mamasshtmll (Bottom) Examples are provided 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 zine binder (pink sphere). Although the entire hydroxamate group is classified according to the SMARTS pattern, the pharmacophore weights are specifically assigned to the two oxygens capable of coordinating zinc ion, which places the pharmacophore between the oxygens rather than at the group centroid.
40

---

Published as a conference paper at ICLR 2025
of inertia of QB (the positive and negative orientations for the two largest principal components). For all other initializations when N > 5, we use Fibonacei sampling to obtain evenly spaced out rotations in S0(3).
We use unit quaternions to represent R to decrease the number of parameters and increase opti- mization 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 2022 in defining the volumetric shape similarity between two atomic point clouds CA and C3 as vol
CA. Cg) with Wa,b = 2.7 and = 0.81. The surface shape similarity between two surfaces Sa and Sn is defined as sim (S AA, Sn) with Wo.b 1, and (12)) Here, V is function fitted to som depending on the choice of V12-
We fit vp (m2) by minimizing RMSE between sim sim for n2 E 50, 100, 150, 200, 300, 400 on IK randomly sampled molecules from https: /www. kaggle.gg....................... com/ datasets/ with the number of heavy atoms ranging from to 62. For parameter fit- ting, used 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 I2 and Vy in experiments, but do not expect this discrepancy to meaningfully affect our analyses. The conformations were generated with RDKit ETKDG and optimized with MMFF94, Results and cor- relations from tuned a's 's found in Fig. 18 We do not expect, nor want, perfect correlation with the volumetric similarity since we aim to decouple the shape representation from the atomic coordi- nates. V/ (12) is defined an quadratic interpolation function from SciPy used fit the optimal a's found in Table 6 Note that for surface similarity, we require n = nss,, due to the tuned a. We show the effects of stochasticity Figure
400 surface points
300 surface points
200 surface points
Heeaey Aloit
Tanimoto Scores
Tonimoto Scores RMSE 0.0391
Tontmoto RMSE 0,0406
RMSE 0.0393
0.
Surface Similarity
0.25
150 surface points
100 surface points
50 surface points
1.0
Tanimoto Scores
Tarimate Scores
Tanimoto Scoms R² - 0.683 0.
0.3
A- 0.8
0,0424
6? 0.775
0.0519
0.
0.2
0.4
68
0.0
0.2
0.4
0.
0.6
Volumetric Similarity
Figure 18: The correlations between surface similarity and volumetric similarity after tuning the Gaussian width parameter a to minimize RMSE between the two scoring functions.
Table 6: Tuned Gaussian width parameters (o) for the surface similarity scoring function. Larger corresponds to smaller width.
Number of surface points (nz)
50
100
150
200
300
400
Gaussian width (a)
0.6011
0.8668
1.022
1.118
1.216
1.258
41

---

Published as a conference paper at ICLR 2025
LoW0
Self-Simi selfSiillarit
6.13
EgP Self Stmillaill
237
5312
2333
-400
539
Number atoms (ng)
Number
Figure 19: The effect of stochastic surface sampling on the surface and ESP scores, measured by the average self-similarity as function of the number of atoms 711 For each of ~900 (BTB-relaxed) generated compounds from ShEPhERD trained on ShEPhERD-GDBI7 17, repeatedly extract the true shape 22 or ESP 33 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 given molecule. We plot the average self simiaattyy as function of the number of atoms in the molecule n when extracting 32 or æg with either 75 or 400 surface points. As the surfaces are more densely sampled - either by reducing the average self-similarity approaches 1.0. Recall that we use 75 points when evaluating the self- consistency of ShEPhERD's jointly generated (æ1, x22) and (æ1, Ø3) samples, but use 400 points when scoring similarities between two molecules all conditional experiments. Using 400 points for shape and ESP similarity scoring gives much more precise estimate of interaction similarity.
for fixed 722 na, or increasing 72 or ni for fixed
A.6.2 ESP SURFACE SIMILARITY SCORING FUNCTION
The similarity between two electrostatic potential (ESP) surfaces simEsp (33,AA B3,B) 3) is defined with an overlap function:
|| [a] u [b]
OSSP ABB
Σ QEQ A DEQu
Σ
exp
2 ||r
ro|l3) exp
2a
is the permittivity of vacuum with units ee A)-1 We chose through manual inspection of optimal alignments. Fig. 19 shows the effects of stochas- ticity. Others have also considered different formulations of ESP scoring functions to characterize
where (n3), and X
0.3
[Amro
bioisosteres (Good al. 1992, Cleves et al. 2019 Bolcato et
2022 Osman & Arabi 2024).
A.6.3 PHARMACOPHIORS SIMILARITY SCORING FUNCTION
For pharmacophore similarity, follow the the vector weighting formulation of PheSA (Wahl 2024) and define the overlap similarity between two sets of pharmacophores 24,4 and #4, as:
A.B.m
Σ
Wa,dim
π
exp
r6l13
20m
2
a aEQA.m BEQB.m
EmEm OA,Bim
sim pharm (F4,A,
EmEM OA,A.m O B:m
-AAB.m B:m
where M the of all pharmacophore types (IM) Np). Om 3(m) where Ω maps each
pharmacophore type to a Gaussian width using the parameters from Pharao Taminau et al. 2008).
Wa.bm is defined as:
if m is non-directional,
Wam
Viai_, vibl
if m is directional.
We take the absolute value of V[a]] effects are symmetric across their plane. Note that Pharao does not define the halogen-carbon bond pharmacophore so we assign 1.0.
T ((b)... for aromatic groups as we assume their interaction
42

---

Published as a conference paper at ICLR 2025
A.7 ADDITIONAL DETAILS ON MODEL DESIGN, TRAINING PROTOCOLS, AND SAMPLING
Feature scaling. Hoogeboom et al. (0022) and Peng et al. (2023) found that linearly scaling certain one-hot features in 3D molecular DDPMs improve model performance measured via sample validity, presumably because increasing/decreaiin the magnitude of certain features causes them become resolved earlicr/later in the denoising process. We loosely follow these prior works to scale a and f by 0.25 and B by 1.0. We additionally scale and p by 2.0. Finally, we scale the pharmacophore directional vectors 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. (023).
Noise schedule. ShEPhERD uses the same noise schedule for all diffusion/denoissngg processes. We choose the variance schedule of to be a weighted combination the linear schedule used by RFDiffusion (Watson et al., 2023) and cosine schedule introduced by Nichol & Dhariwal (202I) Fig. 20 plots our noise schedule in terms oz, Ct,, oe, and a for t [1,T]. We use 400 for all ShEPhERD models.
0.30
0.2
0.4
0.20
0.5
6
15
057
3.1
000
30
300
139
200
156
300
350
400
30
200
130
200
250
300
350
400
Time step t
Time step t
4.5
63
63
Time step
Time step C
Figure 20: Noise schedule used by ShEPhERD for all forward noising processes in terms of at Ors F,, and Ti, for & E 400).
A.7.1 FURTHER DETAILS ON ShEPhERD'S DENOISING MODULES
"E3NN-style" vs. "EGNN-style" coordinate predictions. As described in section 3.3 predict the /=1 coordinate noises (()) att) and 2tt) from the corresponding /=1 node features 2 directly with equivariant feed forward networks that appear in EquiformerV2 and other E3NN-based archi- tectures. In these "E3NN-style coordinate predictions, the noise for each node [ predicted from only the [=1 feature ā [k] of that node. While each [A] implicitly contains information about the node's 's surrounding environment (e.g. neighboring nodes) due to ShEPhERD's joint module, the actual denoising step is not defined with respect the positions of the neighboring nodes. This notably differs from Hoooeboom et al. (2022)'s 's original formulation of equivariant DDPMs for 3D molecule generation, which used "EGNN-style" (Satorras al. 202I) coordinate predictions de- noise the molecule's s coordinates C. For single-layer EGNN, the coordinate prediction step takes the form (using our notation):
C(k] [[j dkj +
C(KY C(k]
Σ jzk
Φc 2[{)), 210) daj
CTN -cbi
43

---

Published as a conference paper at ICLR 2025
where &c is an MLP and C(k]' can be interpreted
updated coordinate for the node k, which
equivariant to E(3) transformations of C. With this updated coordinate, the noise [j(k] may be obtained by computing é; C[k] Clk]. In "EGNN-style" coordinate predictions, therefore, the predicted coordinate/noise for node is defined with respect to the other nodes in the 3D graph.
e()) all) We hypothesize that E3NN-style coordinate predictions are sufficient for predicting S, S and
e(()) Ep (e.g.. the [=1 coordinates of 22, 33, and æ4), these noise predictions do not need to be so highly sensitive to the positions of neighboring nodes. However, for the coordinates of the 3D molecule 21, inter-node geometries (c.g., bond lengths and angles) are guite relevant the quality of the molecule's conformation. Hence, we hypothesize that it may be favorable to use EGNN- style coordinate predictions to obtain éc. However, EGNN-style coordinate predictions expressive than E3NN-style coordinate predictions, as the difference Cīkl' C.k] is just a linear combination Ejvx Wj (C(A) - C(jl) with (leartt) weights wj. To denoise C, thus choose to both E3NN- and EGNN style coordinate predictions. We first predict set of updated coordinates (*(() from C using E3NN-style coordinate predictions, and then refine C"(t) with EGNN-style coordinate predictions to obtain C"(t), The noise è is then computed as & = C C(²)
less
Algorithm outlines this process within the context of the entire denoising module for æ1.
Forward- Barward-noising noising and denoising the covalent bond graph B. Multiple works have found that jointly diffusing the covalent bond graph with the 3D atomic coordinates helps improve sam- ple quality for 3D molecule DDPMs (Peng et al. [223] Vignac et al. (023). While using discrete diffusion for bond/atom types has been found to further improve sample quality compared to diffus- ing/denoising continuous one-hot representations of bond/atom types, choose to use continuous representations of bond/atom types to emphasize overall modeling simplicity, consistency, and ex- tensibility. In this section, clarify details regarding the forward and reverse processes for B.
We represent B E R Xmy as an n1 X n1 adjency matrix attributed with one-hot features of the covalent bond order between each pair of atoms in æ 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 note that treat bonds between pairs of atoms undirected edges by symmetrizing (t in both the forward noising and reverse denoising processes. Namely, we only forward-noise and denoise the upper triangle of (t), and copy the upper triangle to the lower triangle. We mask- out the diagonal of (4) so that self-loops are not modeled.
To maintain permutation invariance when predicting the noise & [k, j] for the bond between atoms k and j from their 7=0 node representations 21 [k] and 21 [j]. we use permutation invvariant MLP:
{k 5l 2 buj + byı
B [k]]k, [RBF(d 2 [k]. [i]]
buj MLP B
byk = MLP
B fRFF =1 [b]]
where dæj dzk- |C00 [k| c bاll,+ Lilll
FABF is radial basis function expansion using Gaussian basis functions, and [...] denotes concatenation in the feature dimension.
44

---

Published as a conference paper at ICLR 2025
A.7.2 TRAINING AND SAMPLING PROCEDURES
3 4 outline the forward pass of ShEPhERD's denoising network and detail how train and sample from ShEPhERD. Algorithm details the forward pass of ShePhERD's ' denoising module for 21. The denoising modules for I2, 33 and 24 are simpler and described in section [
Algorithms 2
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. dessriics overall training times, and App. A.P. lists all hyperparameters.
Each training example requires sampling a time step E [] (we use T 400). Rather than uniformly sampling f within that range, we upsample in the range [50, 250] to accelerate model convergence, since this interval captures the most important part of the denoising trajectory. Specif- ically, we uniformly sample E [50, 250] 75% of the time, / (0,50) 7.5% of the time, and t E (250, 400) 17.5% of the time.
Algorithm Denoising Module for #
1: Given:
2:
(21, 21)
3: C(
4: B(I)
5: Predict noise for /=0 node (atom) features:
6: (Eas εÿ) = (*1)
7: Predict noise for 1=0 edge (bond) features, using permutation-invariant MLP:
b RBF indicates radial basis function expansion with Gaussian functions
8:
bRj MLP
B k. j]. RBF([
C(* (]I)) [1k], =1
9: by: = MIP B B [ kj. BBB(BBCCC []] ( [+]) [k) 2 (1]
10:
é j = 2 (bb, + byk)
11:
k 2 (bbj bjk)
12:
Predict noise for coordinates C: ACEANN = w.3NN_eeed_forraaeeee
13:
14:
C"(1) = ( + AC FANN
15:
C"(()) FGNN(21, C*O)
2(t) = CC'() (t))
16:
(0) (t) 17: return Ca . .
45

---

Published as a conference paper at ICLR 2025
Algorithm Forward Pass ShEPHERD'ss Denoising Network
1: Given:
2: z (t) (al9) f()) C0 B()
3: 22 (6) (
4: 2y (S v(()) ((0))
5: 444 (plt), p(O), V(O)
6: t (time step)
7: Add virtual node to each æ; system with positions at the COM of C(+) (0 if centered):
8: Øi 222 2y Vi add.virtual. (0) (t) (t)
(0) (f)
22 23 x44
9: Embedding modules: embed each x into / and node features:
ž (ar (t)
EquiformerV2 module
11: 22, 22 = 22(x1) (xg
EquiformerV2 module
ø3, Ø3 d3 2z (t)
EquiformerV V2 module
13: 21, 21 = d1
4 EquiformerV? V2 module
14: Joint module:
Collate each system to form a heterogeneous 3D graph with ny + n2 + ng + 24 nodes:
16: Chetero concat ([C(t) S/
s(1) p dim = 0)
Thetero concat [1, 22, #3. 21. dim =0)
18: Shetero concat ( &2, Æ3, 21. dim = 0)
g 20: >hetero _ hetero @loou Chetero: Zhetero, Rhetero)
Encode heterogeneous graph and update node embeddings:
EquiformerV2 module
21: Zheterot hetero
22: Σhetero
ketera
23: Pool I = node embeddings across the nodes of each homogeneous sub-graph:
24: 2ppol sum Ehatero [k] Vk if k E ax
25: z sum Zhotero kk] Vk if k E x;
26: zpoot sum Rhetero [k]
Vk if k x3 (t)
27: xpoot = sum h
[k]
Vk k E (1)
28: Embed global = code for entire system:
29: zalobat
concattttcccccccccc ([=poo spool ppool pooo dim -1)
30: zalobal
o globul (zglobal)
E3NN feed-forward network
ioint
31: Apply E3NN/Equiformer's fully-conmected tensor product (FCTP):
32: Lembed = s i misuuoooii dal pos ional eew.ddicnngg
33: Zjoint, joint = squiforrer.FFrrr SitneerrFCTCCC (ltembeds 3ooint [tembed: joint
2global1,
gglobal)
34: Residually update node embeddings:
35: for each 1, do
36:
2: + jjoint
37:
2; Ejoint
38: end for
39: Denoising modules:
(ēa, Ca, &C, &B) denoise &1)
41: (Esa) denoise = &2)
42: (Esg, denoise x3(23, 2n)
43: (pp,Ep,EV) dencise (24, 24)
44: return (Ea, ε,, &C, ε εS,, εs, eor Ep, εp, εv)
46

---

Published as a conference paper at ICLR 2025
Algorithm Training Algorithm
1: Given:
2: RDKit mol object mol containing 3D conformation molecule with explicit hydrogens,
locally optimized with xTB.
3: Noise schedule oe, [te, Tt. ae for t E [1,7].
4:
Denoising network to be trained
5: Obtain unnoised input molecule and its interaction profiles:
6: mol = rrecenter.molmmme
b> Recenter molecule to have an unweighted COM of
7: a, C. get_x1 (mol)
8: hartial.charger
get
(mol)
single point
calculation with xTB in implicit solvent gas phase
9: vdw radii = get vdw radii(a)
>> list of Van der Waals radii of each atom S2 solvetl wwo.vvet.surfaceeCC vdw radii, radius 0.6, n 2 75)
11: S3 = S2
at htt.ss..aa.surfacee55..... C.aartial.chaagees 13: p. P. V get.ppharmccophhoooo (mol)
12: v
14: Linearly scale certain features:
15: a = scale a(a)
16: f scale f(f)
17: B sca scale_B(B)
v scale_v(t)
19: p scale.p(p)
scale.V(V))
21: Forward-noise input state:
22:
Sample / E [1. and obtain the corresponding ot, 0r, and
23:
Sample independent noise e ~ ~V(00, for = (a. Ja C, B, S2. S3, U, p. P, V)
24:
Subtract center of mass from EC
25:
Symmetrize the noise for EB
26:
Forward-noise each æ [a, f, C,C, S3. v,p, P, V] via: x(1) a.æ + F
27:
Perform forward-pass with ShEPhERD's denoising network:
28:
a (a ((()) (()
42 (0) (
30:
rz (S ,(()))
(6)
31:
v1 ((1), v(t)
p(o)
32:
(Ea: é &c, &B, Es. êsg éur Eps εy, x2 (1) x4
SShePhERo_forwarde
33:
Compute losses:
34:
Lg
E E
Vo S.
35:
L
L
36:
Optimize network parameters:
Take gradient step w.r.t. the denoising network's parameters to minimize
47

---

(8) Published as a conference paper at ICLR 2025
Algorithm 4 Sampling Algorithm for Unconditional Generation
1: Given:
2: Trained ShEPhERD model 7xx that learns P(X) for X C {0 22 33, 24).
3: T1
b number of atoms generate.
4: (if æ4 E X) T44
number of pharmacophores to generate.
5: Noise schedule de. at tt. a for E [1.
6: extra noise ts ordered list of time steps to add extra noise to coordinate states extra_noise.scaless list of scaling factors apply to extra coordinate noise harmonization.ts. ordered list of time steps at which to perform harmonization. harmonization.stege
7:
8:
9:
list time horizons At for each harmonization action.
10: Sample initial states at = 7 from Gaussian priors:
11:
Sample noises xTT) (0, 1) for æ [a, f. C. C. S2. u.p. P.V..
12:
Remove center of mass from ((T).
13: Symmetrize (T))
14: Denoising loop:
t= T
16: while f > 0 do
a (t) (a f())) ((t), (()).
22 (5)
18:
xg (t) (5 (.))).
20:
(p((), p()) V()).
21:
if hharmonizaaionn (s| = then
Resample states At steps in forward direction
22:
harmonizattion.ss.pppll
23:
At ---aamoniiaaiio.stepss. (op(())
z; (+++t) = forward noise( ,, At, noi se schedule) v E [1,2,3,4]
24:
25:
t= At
26:
continue
end if
28:
Predict true noise with denoising network 1)xx
(éa, Ec, ES2, ES₁, & Ep. εv,
ShEPhERO-forwarddllt 32 13 .)
(1) (1)
30:
31:
for E C. Ss. U. do
Apply reverse denoising equation
32:
c ~N,00
33:
if
C then
34:
c = c get center of mass((
35:
end if
36:
if æ E [C, S3. P and extra == then extra notraanniise
38:
else
40:
= 0
end
42:
x (t-1) 1) 1 (t)
o² Ex + (010+-1 + se E
as
a.d.
43:
(1) (t-1)
end for
45:
-=t-l
46:
end while
Final feature processing:
(a()) f()) (0)) (())) (0)) v()) undo.scaling(a (())) (()), v(0) ((()) v(0)
48:
49:
(a (0) (()) (()) v(0)
ergmax.and.round(a" f(0), B()) (()) v(oj)
48

---

Published as a conference paper at ICLR 2025
Algorithm Sampling Algorithm for Conditional Generation with Inpainting
1: Given:
2: Trained ShEPhERD model 7)x that learns P(X) for {x1, 22 23 va}
3: Y1
b number of atoms generate.
4: (if &4 x) n4
number of pharmacophores to generate
5: Noise schedule dt, at dt. a for t [1.
6: A subset the target interaction profiles I S (S2. p",
7: Simulate and store forward-noising of each target interaction profile:
8:
for
x="
S"" do
9:
2 Fnoised = 0 x
10:
(1=00 = x
noosed [0] x (t=00)
12:
for & E range(11 T) do:
Sample e N(0, 1)
14:
Fnoised [t] 0y moised
+ orE
15: end for
16:
end for
17:
Sample initial states at t = T from Gaussian priors:
18: Sample noises (xT) ~ N(0 1) for x = [a, Sa, v.p. P.
19: Remove center of mass from ((T) Symmetrize (2))
20:
21: Denoising loop:
22: t t= T
while do
24:
D Replace denoised profiles with their corresponding forward-noised target profiees
25: for (1) [S S v(t) p((), pio) [()) do
(+) x noised [t]
if @oised E I
27:
end for
z (t) (a(t), f(t)) ((t),((), B(t)
28:
2 ( (5
29:
30:
23 (5 (!))
31:
a = (t( p(t) v(t)
32:
Predict true noise with denoising network 7)x
33:
(ea, &g &c, én, E Es,, & &p. & ShEPhERD (t) x2 33 .))
(t)
34:
Apply reverse denoising equation
for x fa
C.f.C.BB.. f. U.p. do
36:
NN,N1 if C then
37:
38:
c/ = d -feeteceeteroofmassse
39:
end
40:
a (l-1) 1) (t) of (tt)
aco + 0.0 e
41:
(t) x (t-1)
42:
end for
43:
-1-1
44: end while
45: Final feature processing: (a))) f()) B (0)) p0))) v())= = undo 0) (()) B(0) v(0) ((0), V(0))
(a 00) (0) (0) (0)) v(0)-
a rodudd (0), f(o) ( (0) v
49

---

Published as a conference paper at ICLR 2025
A.7.3 MODEL HYPERPARAMETERS
Tists hyperparameters relevant training ShEPhERD. TableBJlist the total number of learn- able parameters for each version of ShEPhERD analyzed in the main text. Table hyperparam- eters related to ShEPhERD's '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 æ2 and Z3, 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 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 be more sensitive 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
722
75
n13
75
surface probe radius
0.6
and f scaling factors
0.25
B scaling factor
1.0
scaling factor
2.0
scaling factor
2.0
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(
4.375,054
ShEPhERD-GDBLI P(æ1,
4,387,407
ShriPhERD-GIDB17 P(x1, ((
4,407,321
ShiEPhRRD-MOSES-.a P(r), 33, 24
6,010,427
50

---

Published as a conference paper at ICLR 2025
Table 9: Hyperparameters for ShEPhERD's denoising network ?
Denoising Network Hyperparameters
Parameter
Value
Default EquiformerV2 Parameters
num node channels
64
Imax list
[1]
mmax_lis
32
grid resoluusotttion
16
num_sphere_samme _samples
128
edge channels
128
activation function
silu
norm_type
layer.norm_sh
use_sep ____ccca act
True
use-grid_mlp
True
use_gate_act
False
use attn renorm
True
use s2 act attn
False
Joint Module Parameters
num Equiforuurvrr 22.ymeerr layers
2
attention channels
24
num.attention.heal
radius crapluussgaphhhuofff
5.0
RBF_cutoo
5.0
TC Embedding Module Parameters
num Eeuiforuurrrrr V2 layers
4
attention channels
32
num www.aattetiionneee heads
fín hhhiddennhhannee
64
radius_grapp wrw.p...pp........
0o (fully connected)
RBF_cutoff
5.0
32 33, 24 Embedding Module Parameters
num EquiformerV2 layers
2
attention attention.channels
num aunnaatention _eads
2
fin hidden channels
32
cadius-graph.cutton
5.0
RBF cuutoff'
5.0
x3 x3_scalar CeBF_expansionnm _min
-10.0
x3 scalar RBF_exe BBBFeexansion. max
10.0
2i Denoising Module Parameters
MLP hidden dim
64
num MLP hidden layers
e3nn_ffn www..niddenncom hidden.channels
32
egnn_normaine nortoolize_vettors vectors
True
egnn wwwstanneeeppanssiiiiiee _ii
32
32g, 3z, B, Denoising Module Parameters
ML.P hidden dim
64
num MLP hidden layers
e3nn_ffn www.:fi..iiddeccchaaeee
32
51

---

Published as a conference paper at ICLR 2025
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ı (Fig. 21). We attribute this behavior to ShEPhERD being unable to reliably break spherical symmetry when denoising the Euclidean coordinate components of ( This phenomenon is particularly evident for models trained on ShEPPERDD-DDB1,, which contains number of spherical cage-like molecules in the true data distribution. But, this also occurs which does not contain such molecules. Intuitively, 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 tra- jectory (e.g.. 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 impossible. Due ShEPhERD's s strict -quivvaqurereennee,, though, the denoising network cannot natively break the spherical symmetry of this point mass; this ultimately leads the network to generate pphericall-ssapper (yet often still valid) molecules.
for models trained on ShEPhERD-MOSES--
we can easily solve this problem by manually adding extra noise to the coordinate components of (t-1)) during the most informative time steps. With our choice of noise sched- ule, the states (() undergo the most significant evolution in terms of their spatial coordinates 80. For t E 130). then, we add extra (isotropic) noise to C(t-1), S ff(t-1) S ((t11) so that (as an example): ( - c(()
Fortunately,
from 130 through
of e + -. + Oedz 100, and scales the noise (although we just use Ç 1). For ((t-1), we still remove the COM from Exymmetry breaking
& breaking where Exymmetry
1).
Naively adding extra noise to the coordinates of the denoised states could quickly cause the model to go out-of -isistiiibutinn, leading to low-quality or invalid samples. To remedy this, we use harmo- nization/resampling (Lugmayr al. 2022) to re-sample the state X ((=80+A1) given the symmetry- broken X (t=80) (we use At 20). Empirically, we find this adequate maintain sample quality, as evidenced by the high validity rates and low RMSDs upon --elrllaatoonn for unconditional sam- ples from ShEPhERD when trained on ShEPhERD-GDBIT
admit there may be less intrusive methods of breaking symmetry; e.g., adding non-isotropic noise with smaller &. How- procedure be simple and effective while only marginally
Nevertheless,
ever, we find our symmetry-breakee increasing sampling time.
Finally, we emphasize that this symmetry breaking procedure is only performed during uncondi- tional 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 aoditionalll-ssmpleed structures).
Without symmetry breaking
With symmetry breaking
Figure 21: Without symmetry breaking during unconditional generation, ShEPhERD tends gen- crate spherical molccules, particularly for models trained on ShEPHERD-GDBIZ. Shown are un- conditional samples from ShEPhERD trained to learn P P(x #2) on ShEPhERD-GDBI7, 17, generated either with or without symmetry breaking.
52

---

Published as a conference paper at ICLR 2025
A.9 CHARACTERIZING THE NEED FOR OUT-OF-DISTRIBUTIN PERFORMANCE
Multiple of our experiments, particularly those related to ligand-based drug design, require ShEP- hERD generate larger molecules with more pharmacophores than the molecules that ShEPhERD was trained on. For instance, the largest molecules in ShEPHERD-MOSES-- contain 27 heavy atoms, whereas the best-scoring ligands generated in our natural product ligand hopping experi- ments 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 pharmacophore profile containing 27 phar- macophores, which is many more than those contained by the molecules in SHEPHERD-MOSES- To quantitatively illustrate the need for out-of-distributin performance and ShEPhERD's s extrapola- tion ability, Fig. [2 shows the distributions of 1 vs. ng for molecules from ShEPhERD MEOEE----SSSAAE 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 ShE:PEERD finds high-scoring compounds that both in-distribution and meaningfully out- of-distribution. Fig. (2]ssowss similar distributions for our unconditional and conditional generation experiments on the ShPhERD-GDBI7 dataset.
Legend
Hi+ +++
25
25--
Waarmsccpphrrecooom Counrt
CHONO 300000
Phormacaphos Count
de
20.
-20
9,4
Sord
1.5
1.0-
10
morging
20
40
80
20
40
60
80
Atom Count
Atom Count
Figure 22: Distributions of the number of pharmacophores relative to the number of atoms (includ- ing hydrogens) for conditional samples from ShEPhERD when trained ShEPhERD-MOSES--e For comparison, the kernel density plot of a randomly sampled subset of ~200 molecules from the ShEPhERD-MOSES-e dataset is shown in blue on both figures. (Left) Conditional ShEPhERD-
generated samples from P(x1 |2:3, Z44) 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 n1 while specifying n4 based on the target 24, the ShEPhERD -generrtd samples appear as horizontal lines. (Right) "Top-10" scor- ing 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 both the lowest energy conformer (LE) and the bound pose of each PDB ligand. Top-scoring samples are found both in- and out-of -istrrribbttoon in terms of (n), n4)..
53

---

Published as a conference paper at ICLR 2025
20
Distributions ShEPhERD-GDSS
Pharmacophore Count
15
Unconditional generation Conditional generation
10
20 40
60
Atom Count
Figure 23: The distributions of the number of pharmacophores relative to the number of atoms (including hydrogens) for molecules the ShEPhERD-GDD GDB17 dataset and molecules sampled by ShEPhERD trained to learn P(x1, We include the distributions for ~1M randomly sam- pled molecules from the ShEPhERD-GDBI dataset, P(æ11 at4), and conditionally-generaeed molecules from P(æ 244. All distributions are visualized as kernel density estimate plots. All distributions show significant overlap.
coanditionallygennrraae
molecules from
54

---

Published as a conference paper at ICLR 2025
A.10 ADDITIONAL EXAMPLES OF GENERATED MOLECULES
Figures 24 25 and 26 show additional random examples of unconditionalll--eneratedd molecules and their jointly generated interaction profiles from versions of ShEPhERD trained on ShEP- hERD-GDBI7, Figures 27 28 and 29 show additional random examples of ShEPPERD-generane molecules when conditioning on target interaction profiles via impainting Although we are primarily interested in using the models trained to learn P(11, 33, (44) on ShEPhERD-MOSES for condi- tional generation tasks relevant ligand-bases drug design, we also show unconditional samples from this model in Fig. 30
Figure 24: Random unconditional samples from ShEPhERD trained learn P(21,02) on ShEP- hERD-GDB17. We show the natively generated molecular structures (without xTB-relaxation) and their jointly generated shapes. We show examples for both large and small n1
Figure 25: Random unconditional samples from ShEPhERD trained to learn P(æ1, xs) on ShEP.
hERD-GDBI7. show the natively generated molecular structures (without xTB-relaxation) and their jointly generated ESP surfaces. We show examples for both large and small
Figure 26: Random unconditional samples from ShEPhERD trained to learn P(æ m2A) on ShEP. hERD-GDBI7. We show the natively generated molecular structures (without xTB-relaxation) and their jointly generated pharmacophores. We show examples for both large and small nı. is sampled from the empirical data distribution P(na hydrophobes (the most common pharmacophore in ShEPhERD GDD--DDBB77 17): purple and blue arrows represent HBAs and HBDs, respectively; orange arrows represent aromatic groups; and green ar- represent halogen-carbon bonds.
these images, grey spheres represent
55

---

Published as a conference paper at ICLR 2025
Figure 27:
Conditional samples from ShEPhERD trained to learn P 22) on ShEPhERD- GDB17 - obbtaned by inpainting the 3D molecule given target shapes extracted from molecules held- out from training. We show the natively generated molecular structures (without xTl-relxxttoon or realignment) overlaid the target shapes. The shape surfaces upsampled for visualization.
Figure 28: Conditional samples from ShEPhERD trained to learn P 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 KTB-relaxation realignment) overlaid the target ESP surfaces. The surfaces are upsampled for visualization.
Figure 29: Conditional samples from ShEPhERD trained to learn P(z), x.) 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 <TB-relaxation or realignment) overlaid the target pharmacophores.
56

---

Published as a conference paper at ICLR 2025
Figure 30: Random unconditional samples from ShEPhERD trained to leamn P( D3, Fa) on ShEPHERD-MOSES-ag We show the natively generated molecular structures (without xTB- relaxation realignment) and their jointly generated ESP surfaces and pharmacophores. Note that a number of these structures would not pass our validity criteria.
57

---

Published as a conference paper at ICLR 2025
A.11 TRAINING AND INFERENCE RESOURCES
Training. We train all models with V100 GPUs (32 GB memory). For both the SAEPHERD--DBBIT and SIEPhERD-MOSSES
datasets, we train each of the , ,( and P( 33 .:4 models on 2 GPUs for approximately 2 weeks. We train the P(21, 2A4) model on 1 GPU for ap- proximately 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 22), P(æ Fz), and 1 23 244) models, generating batch of 10 independent samples (either unconditionally or via inpainting) takes approximately 3-4 minutes on a V100 GPU, when using iterative denoising with = 400 steps. Sampling times are reduced by 50% for the P( (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 system hardware and the choice of model.
Long sampling times are somewhat intrinsic to the DDPM paradigm due to needing numerous for- ward 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 acceler- ate sampling (e.g., by using fewer denoising steps inference, reducing T, using less expressive networks, reducing n2 and ng, etc.). We leave such performance engineering to future work.
58

---

Published as a conference paper at ICLR 2025
A.12 ADDITIONAL EXPERIMENTS AND COMPARISONS TO RELATED WORK
A.12.1 COMPARISON TO SQUID FOR SHAPE-CONDITIONED GENERATION
We compare ShEPhERD to two ligand- based models, SQUID (Adams & Coley 2022 and Shape- Mol (Chen et al.] 2023). for shape-conditioned molecular design. To perform this comparison, we apply the version of ShEPhERD trained on MOSES-aq to learn P(æ I 3a, (44 to generate new molecules conditioned on the surface shape (e.g., the coordinates of æa which implicitly define x²) of the 1000 3D molecules in SQUID's 's test set. Note that SQUID and ShEPhERD use the same train- ing and test splits of MOSES, and hence there was no data leakage. We condition the P( 23 21) model on the molecules' shapes (only) by inpainting just the coordinates of Ø3, and allowing ShEP- hERD to generate T1, T., and the electrostatic potentials of 233 We generate up to 10 samples from ShEPhERD and compare the average volumetric shape similarity using ROCS (the same shape scoring function that SQUID used). Since SQUID's generated molecules for each test molecule are publicly available for download via their GitHub, we re-evaluated the average shape similarity of molecules to ensure fair comparison.
Across the 1000 test-set molecules, the average shape similarity of SQUID-generen molecules is 0.70 when using X= = 1.0, and 0.74 when using À = 0.3. In contrast, the average shape similarity molecules is 0.80, substantial improvement. Moreover, the ShEPhERD. generated molecules have an average 2D graph similarity of just 0.22 compared to the reference molecule, whereas the SQUID-generated molecules have average graph similarities of 0.25 (A = 1.0) and 0.35 (A 0.3). Hence, ShEPhERD is able to generate less chemically similar molecules that have higher shape similarity to the target, on average. The full results are found in Table A.1.1
of
We emphasize that ShEPhERD greatly exceeds the performance of SQUID even with this indirect way of doing could also train 2() model directly for this task; however, we did not have time to retrain ShEPhERD from scratch within the rebuttal time period. We also note that ShEPhERD condi- tions on surface shapes, whereas have compared ShEPhERD to SQUID using volumetric shape similarity (with ROCS). We expect that we could further improve ShEPhERD's performance in shape-conditioned generation by retraining ShEPhERD with smaller surface probe radius, thereby "shrinking" the surface shape to better model volumetric shapes.
generation, i.c. via inpainting the joint P(x1, 33, 24) model. We
Table 10: A comparison of ShEPhERD to default versions of SQUID (A = (1.0, and Shape- Mol. We recompute the shape and graph similarities of SQUID using the structures provided by their GitHub. Note that the values for the evaluations of ShapeMol and ShapeMol+g were directly copied from (Chen et al., 2023) which use a different alignment protocol.
Model ShEPhERD
avgSim, (1)
avgSim, (Q)
QED (t)
0.799 + 0.058
0.223 0.065
0.723
SQUID () 1.0)
0.700 + 0.107
0.245 + 0.078 0.078
0.760
SQUID (A- = 0.3)
0.735 + 0.117
0.346 + 0.161
0.766
ShapeMol
0.689 + 0.044
0.239 + 0.042
0.748
ShapeMol+g*
0.746 0.036
0.241 + 0.050
0.749
A.12.2 COMPARISONS AGAINST MARRUCTRREBBSSEE DRUG DESIGN MODELS THAT EXPLICITLY ENCODE THE PROTEIN POCKET
While ShEPhERD is a ligand-only model and does not explicitly model protein pocket geometries, we were asked to compare to structure-based drug design (SBDD) models. We emphasize that since the methods differ in their respective applications and modeling strategies, the comparisons cannot be completely fair. Nevertheless, although ShEPhERD is not designed for structure- drug design, it can still be roughly applied to SBDD by conditioning on the bound pose of a known ligand (but still ignoring the surrounding protein pocket). This bound pose can be obtained from crystallography (e.g.. the co-crystal ligands in PDB entries) or from simulated docked poses (e.g., conditioning ShEPhERD on the top hit from small docking baseegg-based virtual screen). In order to compare to SBDD generative models, we compare ShEPhERD against the benchmarking results presented by (Zheng al. 2024), which report the average Vina docking scores amongst the top-
59

---

Published as a conference paper at ICLR 2025
10 generated compounds (amongst ~ 1000 generated samples) for PDB systems using generative SBDD methods like 3DSBDD, AutoGrow4, Pocket2Mol, PocketFlow, and ResGen. Specifically, use ShEPhERD to condition
the electrostatic potential surface (æz) and pharmacophores (z4) of the PDB ligand in its crystallographic hound pose, generating up to 1000 valid molecules per PDB target. We dock the 1000 generated molecules (starting from their SMILES strings) using the same Vina program (and version) as used by Zheng et al. (2024), and compute the average Vina scores amongst the top-10 molecules per target. These average top-10 scores are compared against the results for the SBDD models evaluated by Zheng et al. (2024) in Table 11 We also evaluate the 2D graph similarities of ShEPhERD ShenhtRD- compounds compared the PDB ligand in each case, which are generally quite low (Fig. 31
Docking scores of ShEPhERD-generee molecules and redocked PDB ligands, com- pared to ligands generated by various SBDD generative models. Scores for 3DSBDD, AutoGrow, Pocket2Mol, PocketFlow, and ResGen copied from Zheng al. (2024).
Table 11:
Model
liep
3eml
3ny8
4rlu
4unn
5mo4
7111
Redocked PDB ligand
-9.58
-9.14
-8.62
-8.70
-9.24
-9.95
-9.10
ShEPhERD
-12.25
-10.80
-10.69
-10.19
-10.60
-10.6
-9.44
= 0.18
+ 0.23
+ 0.23
= 0.23
0.27
2 0.20
+ 0.19
3DSBDD
-9.05
-1002
-10.10
-9.80
-8.23
-8.71
-8.47
: 0.38
+ 0.15
+ 0.24
$ 0.55
0.30
+ 0.45
+ 0.18
AutoGrow4
-13.23
13.03
-11.70
-11.20
-1.114
10.38
-8.84
0.11
0.09
= 0.00
# 0.00
0.12
# 0.27
= 0.33
Pocket2Mol
-10.17
12.25
-11.89
-10.57
-2.20
-10077
-9.74
0.53
0.27
+ 0.16
* 0.12
0.34
* 0.62
+ 0.38
PocketFlow
-12.49
-9.25
-8.56
-9.65
-7.90
-7.80
-8.35
0.70
0.29
+ 0.35
+ 0.25
0.78
1 0.42
+ 0.31
ResGen
-10.97
-9.25
-10.6
-11.75
-9.41
-10.44
-8.74
- 0.29
+ 0.95
+ 0.42
- 0.42
+ 0.23
# 0.39
- 0.24
Even though ShEPhERD only sees the PDB ligand and not the target protein pocket, the average top- 10 Vina scores for ShEPPERDD-eneeeeeeeeeeee -gnerrttdd molecules are very comparable to those obtained by these SBDD methods. Out of the 6 generative models (ShEPhERD and 5 SBDD methods), ShEPhERD ranks 1/6 for one target, 2/6 for one target, 3/6 for three targets, and 4/6 for two targets. These results confirm that ShEPhERD shows great promise for SBDD even though it is currently a ligand-only model.
3eml
3ry
4duu
200
330
Irocurory
Recurrtry
230
100
200
3.0-
2.0
1.0
0.0
02
0.4
0.5
0.6
3.0
3.0
0.2 5/4 0.6
1.0
Rngierpiiit Similarity
20 Rrigersf prineeiinttt..... Similartty
Frigeprint simibrity
28 Rigurpprri Similoiiyy
Aunn
5mo4
7111
300
300
100
Deswwery
200
230
200
Docking screen (random)
SnIPHERRD bound pose trandom 1k)
200
120
02
0.5
Figure 31: The graph similarity distributions of the 1000 valid SEEPPERD-ggeeeaane molecules used in App. 12.2 for each target. We also show the distributions of 1000 randomly selected samples from the 10k docking screen. Note that only 672 ShEPhERD-generated molecules were valid (and docked) for the target liep.
60

---

Published as a conference paper at ICLR 2025
1lep
3eml
3ry
4rtu
150
60
130
Freey vetcy
-00
20-
50
12-
20-
0.0
02
044
0.5
0.6
1.0
Similattty
Aunn
5mp
100
120
80-
Proqurncy
75-
mequenc
Docking screen trandom) pose
-6
50-
SNEPHERD lewest energy
0.0
0.2
2.4
0.5
0.6
1.0
0.0
02
0.4
0.5
0.6
1.0
6.2
0.4
Rngerpite Smilariy
20
Rigerpent Simbartty
20 Ennerpren Similerity
Figure 32: The graph similarity distributions of the valid CHEPhERD-ggneratee molecules used in the main text. We also show the distributions of 500 randomly selected samples from the 10k docking screen.
A.12.3 COMPARISONS AGAINST INPAINTING WITH DIFFSBDD
We also directly compare ShEPhERD to DifTSBDD, which also uses diffusion with inpainting to generate molecules, albeit in a structure-based drug design scenario where the pocket is directly encoded by the model. To compare ShEPhERD specifically against DiffSBDD, we applied ShEP-
hERD to generate analogs given the crystallographic bound pose for the ligand in the PDB entry 4tos (ligand PDB ID: 355), which DiffSBDD was evaluated in Figure of [Schneuing et al. (022). To ensure a fair comparison, we downloaded the 100 molecules generated by DiftSBDD for this PDB target, which are fortumately made publicly available by |Schneunng et al. (2022). We used the molecules generated by the all-atom version of DifTSBDD that uses inpainting, as this is the model that generally performs the best according to Schhneuing et al. (2022). We then re-docked these 100 molecules (starting from their SMILES strings) with Auto Vina (v1.1.2) to obtain their distri- bution of docking scores. We then applied ShEPhERD to generate 100 valid molecules given (X3, x4) of the bound ligand (still ignoring the actual protein pocket), and docked those molecules, again starting from their SMILES strings.
We compare the distributions of docking scores between the molecules generated by ShEPhERD and DiftSBDD in Figure 33 Even though ShEPhERD does not explicitly model the pocket struc- ture, the average Vina score for ShEPhERD -geerrtteee molecules is approximately kcal/mol lower (better) than the average Vina score for the DiffSBDD--eenerated molecules. Moreover, the ShEP- hERD-generated molecules have significantly lower (better) SA scores, as well. Mcanwhile, the ShEPhERD- -egerrtee molecules still have quite low graph similarity to the PDB ligand (vast major- ity are below 0.2). We also emphasize that this conditioning PDB ligand is clearly distrr-iistiiittionn for ShEPhERD, as the PDB ligand contains 42 non-hydrogen atoms, whereas ShEPhERD is only trained on molecules containing up to 27 non-hydrogen atoms. However, we do want to emphasize that, as previously mentioned, comparisons to SBDD methods are inherently unfair. In this case, the PDB ligand itself has a quite good docking score, which means that the ShEPhERD -eeneratdd analogs are likely to also have good docking scores, as ShEPhERD generates molecules that pre- serve 3D interactions. Nevertheless, this experiment shows ShEPhERD's 's utility in SBDD, despite being a ligand-only model.
1.14 COMPARISONS AGAINST SYNFORMER ON NATURAL PRODUCT LIGAND HOOPING
We compare against SynFormer (Gao et al. 2024 in our natural product ligand hopping experi- ments. However, note that the goal of our natural product ligand hopping experiments was to gener- ate small-molecule analogues of structurally complex natural products that (1) preserve their elec- trostatic and pharmacophore 3D interactions, and (2) have simpler chemical structures as assessed via SA score. SynFormer is non-3D generative model that can be used to generate synthesizable molellc analogues of an unsynthesizable reference molecule that have similar 2D chemical structures. We applied SynFormer using its publicly available code in its default settings generate
61

---

Published as a conference paper at ICLR 2025
All compounds
-6
ShEP1ERD bound DInSaDD
Vina Score (kcal/mol)
0.25
-B-
8
Redocked PDB ligand
Similarity
-10
SA Score
0.20
0.15
14
2D
0.10
0.05
4tos
4tos
Top-10
Vina Score lccal/mol)
Similarity
-13
0.20
5A Score
4
0.15
-15-
SNEPNERD
2-
DinSaDD
0.10
-16
Redocked PDB ligand
4tos
4Los
4tos
Figure 33: A comparison distributions for molecules generated by DifSSB and ShEPhERD for the target 4tos (PDB ligand ID 355). We show distributions of Vina Score (v1.1.2), SA score, and 2D graph similarity. (Top) The distributions all 100 valid molecules generated by both models. (Bottom) The distributions of the top-10 10 docked molecules from both models.
analogs for each of the three natural products. We took each of the generated analogs, generated up to 10 xTB- XTB-optimized conformers, and scored their 3D similarities to the reference natural products after optimal alignment.
We compare the SynFormer-geneeee 34 on the basis electrostatic and pharmacophore similarity to the reference natural products. In all cases, SynFormer underperforms ShEPhERD, which is not surprising since SynFormer was not designed to preserve 3D interaction similarity. We also note that very few (if any) of the SynFormer- generated molecules have SA scores below the threshold indicating that SynFormer's 's analogs are generally more complex, even if purportedly synthesizable.
analogs against the ShEPhERD-
analogs in Figure
enforce ShEPhERD (SA<4.5),
62

---

Published as a conference paper at ICLR 2025
Natural product (1)
Natural product (2)
Natural product (3)
Phannacophor Similarit nt
Charmacophore Similarity
Sp varl
Syslime
Syrtomer
Phamarophore
8:20
to te
6.13
2.00
04
05
07
ESP Surface Simileriy ity
ESP Surface Similarity
ESP Surface Similarity
54 MOON
SynFormer ShERLERD
SynFormer
ShERhEAD
Syntorr mer
SHERLERD
Figure 34: (Top) Plots of ESP vs. pharmacophore similarity for the molecules generated by Syn- Former and ShEPhERD in our natural product ligand hopping tasks. (Bottom) SA score distribu- tions for each model's generated molecules.
63


Reviews:
[{"review_id": "MW2n1oeP5f", "reviewer": "Reviewer_HW8X", "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.\n\nI 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.", "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)\n- Have you tried an ablation of the joint noising regime? I believe in DiffSBDD its benefits was found to be marginal.\n- Could your pharmacophore similarity metric be seen as related to interaction fingerprint similarity? \n- 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?\n- The fragment merging protocol is nice, can a semi-automated pipeline be provided in the code?\n- 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?\n- 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).", "limitations": "", "rating": 8, "confidence": 5, "soundness": 2, "presentation": 4, "contribution": 4, "strengths": "- The paper is comprehensive and well-written, with excellent figures and molecular visualizations.\n- 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.\n- 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.\n- 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).", "weaknesses": "**Main limitations:**\n\n- 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.\n    - 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].\n    - Compare to DiffHopp [4]/TurboHopp [5] in case of hit diversification. Obviously these models are only in the structure-based context.\n    - Compare to ChemProjector [6]/SynFormer [7] incase of ligand hopping.\n- Furthermore, the metrics could be greatly improved:\n    - Heavy reliance on only Vina scores and shape/tanimoto similarity. \n    - 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?\n    - 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?\n    - PoseBusters [9] validity could also be used to access the quality of the structures.\n- 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.\n- Would be nice to see anonymous code in the review stage.\n- 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.\n- 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.\n- 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?\n- 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?\n- Comparison to REINVENT.\n    - Was REINVENT trained to convergence? 10,000 oracle calls seems quite low in my experience? (Happy to be proved wrong)\n    - 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).\n    - Little measure of synthetic accessibility in the work. REINVENT does surprisingly well on SA score in my experience. Can this be benchmarks quantitively?\n\n**Minor points:**\n\n- 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)\n- 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.\n- 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.\n- What was the intuition for scaling the pharmacophore feature vector by 2.0?\n- Can the target profile guidance be more or less ‘weakly’ enforced?\n- 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)\n- Have you tried an ablation of the joint noising regime? I believe in DiffSBDD its benefits was found to be marginal. [1]\n- Could your pharmacophore similarity metric be seen as related to interaction fingerprint similarity? \n- 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?\n- The fragment merging protocol is nice, can a semi-automated pipeline be provided in the code?\n- 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?\n- 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)..\n    - Can you provide a list of the 7 PDBs used to perform these experiments?\n- Figure 1 + 2:\n    - Are both very good: I would consider merging to safe space and add more analysis/discussion.\n- Figure 3:\n    - Can ‘dataset samples’ be added to the histograms in the right of the figure?\n- Figure 4:\n    - 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).\n    - 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.\n    - The PDB targets are quite well chosen (Kinases/COVID MPro etc), adding target names in brackets next to PDBs would be helpful.\n\n\n1. Schneuing, Arne, et al. \"Structure-based drug design with equivariant diffusion models.\" arXiv preprint arXiv:2210.13695 (2022).\n2. Ziv, Yael, Brian Marsden, and Charlotte M. Deane. \"Molsnapper: Conditioning diffusion for structure based drug design.\" bioRxiv (2024): 2024-03.\n3. 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).\n4. Torge, Jos, et al. \"Diffhopp: A graph diffusion model for novel drug design via scaffold hopping.\" arXiv preprint arXiv:2308.07416 (2023).\n5. Yoo, Kiwoong, et al. \"TurboHopp: Accelerated Molecule Scaffold Hopping with Consistency Models.\" arXiv preprint arXiv:2410.20660 (2024).\n6. Luo, Shitong, et al. \"Projecting Molecules into Synthesizable Chemical Spaces.\" arXiv preprint arXiv:2406.04628 (2024).\n7. Gao, Wenhao, Shitong Luo, and Connor W. Coley. \"Generative Artificial Intelligence for Navigating Synthesizable Chemical Space.\" arXiv preprint arXiv:2410.03494 (2024).\n8. 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.\n9. 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "bkHKHrBhHE", "reviewer": "Reviewer_pbc2", "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.", "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?\n\nOther 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)?", "limitations": "", "rating": 10, "confidence": 5, "soundness": 4, "presentation": 4, "contribution": 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.", "weaknesses": "Most of the weaknesses of the paper are somewhat superficial and do not limit the importance of this contribution.\n\nAll 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.\n\nThe 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.\n\nThe 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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}, {"review_id": "vnab611gKQ", "reviewer": "Reviewer_JQjD", "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.", "questions": "My questions are the same as mentioned in Weaknesses.", "limitations": "", "rating": 6, "confidence": 4, "soundness": 3, "presentation": 3, "contribution": 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.\n2. The authors achieved credible results in three specific tasks, partially demonstrating the method’s effectiveness.", "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.\n2. The method lacks comparisons with other methods, casting doubt on the superiority of the results.\n    (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.", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "extra_scores": {}}]

Author rebuttal:

Decision: Accept (Oral)