1 Introduction↩︎

Ligand-based drug design (LBDD) uses known bioactive compounds to guide the design of new molecules with similar three-dimensional (3D) shapes and physicochemical properties [1]. In contrast to structure-based drug design (SBDD), which relies on accurate protein structures and is often more computationally demanding, LBDD remains effective when target structures are unavailable and is therefore widely used for hit discovery and hit diversification [2]. Central to LBDD is the pharmacophore, which specifies the spatial arrangement of interaction features—such as hydrogen bond donors or acceptors, charged groups, and aromatic centers—responsible for protein-ligand complementarity.

A common LBDD strategy extracts pharmacophores from known ligands and uses them in virtual screening to identify molecules with similar shapes and features [3]. Methods such as ROCS follow this approach and have shown competitive or better performance than structure-based methods [4]–[6]. However, the rapid expansion of accessible chemical space [7] renders exhaustive screening increasingly computationally prohibitive. Generative modeling has emerged as a promising alternative to brute-force search, allowing for faster inference times and the proposing of molecules outside of predefined spaces. Conditional generative models for 3D small molecule design span a range of families, from autoregressive [8] to diffusion [9]–[12], flow matching [13], [14], and variational autoencoders [15]. While performing well at capturing distributions of desired molecular motifs, these models often produce synthetically inaccessible molecules [16], [17].

Synthesis-aware generative methods for designing new molecules have typically either treated synthesizability as an additional property to optimize [18]–[20] or enforced it by generating molecules via synthetic pathways/reactions rather than as unconstrained molecular graphs [21]–[27]. However, these approaches are not conditioned on three-dimensional chemical features, which are critical for protein-ligand binding. As a result, bioactive ligand generation is typically framed either as a reward-based problem using binding-affinity proxies [23], [27] or as conditioning on known ligands via chemical space projection [28].

We bridge the gap between 3D chemical feature conditioning and synthesizable molecule generation with SynLaD(Synthesis-aware Latent Diffusion), a conditional latent diffusion model that decodes latents to molecules in both 3D and synthesizable space. An overview of SynLaD is shown in Figure 1. At its core, SynLaD employs a molecular variational autoencoder with a dual-decoder architecture: a single encoder maps a 3D molecule to the latent space, from which one decoder reconstructs the 3D atomic coordinates and corresponding atom types, while a second autoregressive decoder generates an associated synthesis pathway. Joint training of these decoders couples spatial fidelity and synthesizability in the learned latent representation.

Conditional generation is then performed in the learned latent space using a Diffusion Transformer (DiT) [29], [30]. The compressed, semantically rich latent space enables efficient sampling and improved conditional control [29], [31], while also allowing SynLaD to generate structure-informed embeddings that the synthesis head decodes into synthesis pathways. This also introduces variability in the synthesis decoder conditioning input, which we demonstrate is essential for ensuring sample diversity. Contrary to recent synthesis-constrained generative models [26], [27], [32]–[34] which rely on reaction templates, SynLaD uses a reaction prediction model trained on diverse, complementary chemical transformations to infer products. By removing the reliance on fixed templates, our approach mitigates the template bottleneck that can restrict designs to a narrow, predefined chemical space, and supports broader exploration via more general molecular-graph edits, and the possibility of extrapolation.

We evaluate SynLaD using standard molecule quality metrics, including pose validity and synthesizability, under both unconditional and pharmacophore-conditioned generation. We further assess its practical utility on bioactive analogue generation for ligands from the Lit-PCBA benchmark [35]. Across settings, SynLaD matches the molecule-quality profile of its separately trained counterparts while consistently improving the fraction of samples with plausible synthetic routes, and it outperforms existing baselines in producing synthesizable, diverse hits under pharmacophore-guided sampling.

The contributions of this work are thus: (1) a dual-constrained latent space that couples two complementary but often competing representations—3D molecular structure and synthesis plans—encouraging latents to remain simultaneously decodable into valid geometries and executable routes; (2) a pharmacophore-conditional latent generative model that enables targeted, 3D-aware design; and (3) an end-to-end evaluation protocol that jointly reports diversity, synthesizability, and target-conditioned hit discovery trade-offs—explicitly contrasting our conditional generation method with library screening and other amortized and non-amortized generative baselines. We release the code for SynLaD at https://github.com/prescient-design/synlad.

2 Background and related work↩︎

2.0.0.1 Structure- and ligand-based drug design.

Many generative models have been proposed to sample small molecules in sequence [36]–[38], 2D [39]–[41], 3D [42]–[46], and voxel spaces [47]. Conditional generation for protein binders can be broadly categorized into pocket-conditioned approaches, which condition on binding-site representations and implicitly capture protein-ligand interactions, and interaction-conditioned approaches, which condition directly on desired interaction cues such as 3D shape, pharmacophores, or electrostatic potential surfaces. The latter is most relevant to our work. Several methods generate SMILES strings or molecular graphs conditioned on 3D pharmacophore representations [15], [48]–[50]. MolSnapper [51] and, more recently, ShEPhERD [52] instead condition on richer interaction profiles and generate molecules directly in 3D space. However, MolSnapper supports a limited set of pharmacophore types, and neither method explicitly constrains generation to synthesizable chemical space.

2.0.0.2 Synthesis-aware generation.

Parallel work has focused on enhancing synthesizability by generating molecules via reaction pathways. Composing these reaction pathways into full synthesis plans is a challenging discrete search problem, and existing methods span diverse paradigms, including reinforcement learning [23], [53], search [54], surrogate-guided optimization [55], genetic algorithms [56], and GFlowNets [27], [32], [57]. One particularly effective strategy is to cast synthesis plans as sequences and learn serialized pathway representations autoregressively [25], [28], [33], [34], [58], analogous to next-token prediction in language modeling [59]–[61]. By amortizing the cost of discrete search, such models learn mappings from continuous representations to synthesis plans and can be integrated into larger frameworks (e.g., encoder-decoder architectures) for tasks such as molecular optimization or retrosynthesis (i.e., predicting a route for a given target molecule). While these approaches enforce viable synthesis routes, they do not directly control final 3D geometry. [62] recently proposed a framework to train an autoregressive synthesis pathway generator conditioned on composite properties, including pharmacophore-encoding fingerprints [63] and molecular substructures. However, the model does not condition on raw 3D pharmacophore features. Recently, CGFlow [64] and SynCoGen [65] jointly design synthesis pathways and 3D poses, but they rely on fixed reaction template libraries and work in ambient space.

2.0.0.3 Latent diffusion models.

Latent diffusion models [29], [66] perform diffusion [67], [68] in the learned latent space of an autoencoder rather than directly in the high-dimensional input space, enabling more efficient training and sampling. This paradigm has been highly effective in image, audio, and video generation, especially when combined with Diffusion Transformers (DiTs) [69], which show that standard transformer backbones scale effectively as denoisers. In molecular modeling, [70] introduced latent diffusion in the space of an equivariant autoencoder, while [30] proposed a (non-equivariant) latent diffusion model for small molecules and periodic materials. However, these approaches do not demonstrate latent generation conditioned on rich, chemically grounded features. We address this gap by enabling pharmacophore-conditioned latent generation and augmenting our autoencoder with an auxiliary synthesis-decoder head that produces explicit reaction pathways, jointly targeting native-3D conditional generation and synthesizability.

3 Methods↩︎

We present SynLaD—a latent diffusion [29] model for de novo small molecule generation conditioned on pharmacophore features (see Figure 1). SynLaD is trained in two stages: In stage 1, a variational autoencoder [71], [72] encodes molecules into a shared latent space and a two-headed decoder reconstructs both the 3D molecular structure and a synthesis pathway. In stage 2, we train a Diffusion Transformer [69] to generate new samples from the latent space, which are decoded into both 3D molecules and synthesis pathways. The two decoders are trained to reconstruct the same molecule under two different representations and at inference time the synthesis decoder is used to ensure the synthesizability of generated designs.

We hypothesize that this two-stage learning process offers two advantages. First, the shared latent space aligns the two modalities—molecular conformations and synthesis pathways—by requiring both decoders to reconstruct the same molecule from a single latent representation. Second, as shown in [30], operating in a latent space reduces the complexity of generation and enables explicit pharmacophore conditioning by decoupling discrete and continuous spaces.

3.1 Multi-headed autoencoder for multi-modal reconstruction↩︎

Our autoencoder consists of one encoder and two decoders trained using a combined reconstruction loss with a regularization term. We base our implementation of the 3D encoder and decoder on [30] and represent molecules using discrete atom types \(\boldsymbol{A}=\{a_i\}_{i=1}^{N}\in \mathbb{Z}^{1\times N}\) and continuous coordinates \(\boldsymbol{X}=\{x_i\}_{i=1}^{N}\in \mathbb{R}^{3\times N}\), where \(N\) is the number of atoms.

3.1.1 3D Encoder and Decoder↩︎

Given a molecule represented as above, the encoder \(f\) projects \(\boldsymbol{A}\) and \(\boldsymbol{X}\) into a per-atom latent representation \(\boldsymbol{Z}=f(\boldsymbol{A},\boldsymbol{X})\), and the 3D decoder \(g_{\textrm{3D}}\) reconstructs the molecule from the latent, giving \(\boldsymbol{A}',\boldsymbol{X}'=g_{\textrm{3D}}(\boldsymbol{Z})\), where \(\boldsymbol{Z}=\{z_i\}_{i=1}^{N}\in\mathbb{R}^{d\times N}\). We use a standard transformer [59] and learn molecular symmetries via random rotations; further details on how atoms and dimensions are sampled at training/inference time is provided in Appendix 8.1.

3.1.2 Synthesis Decoder↩︎

We define an auxiliary decoder \(g_{\textrm{syn}}(\boldsymbol{Z})\), which, given a latent \(\boldsymbol{Z}\), generates a synthesis pathway for the molecule represented in 3D by \((\boldsymbol{A}, \boldsymbol{X}) = g_{\textrm{3D}}(\boldsymbol{Z})\). A synthesis pathway can be represented using a directed acyclic graph (DAG), where nodes represent molecules and edges represent reactions (see Figure 1). Such a DAG shows how parent building block molecules (source nodes) react to first form more complex intermediate products and eventually a final product (a single sink node). To model such a DAG we can serialize it, in a bottom-up manner similar to [58], into a linear sequence of tokens following the recipe in Appendix 7.2. Once we have such a synthesis sequence \(\mathbf{S}=(s_1, \dots, s_L)\), we train a causal transformer model to generate it by autoregressively modeling the conditional probability \(p(s_i|\boldsymbol{S_{<i}}, \boldsymbol{Z})\). We embed each token in the sequence using a custom embedding scheme (see Appendix 7.3) and train the transformer using teacher forcing conditioned with cross-attention on the latents.

3.1.2.1 Product identity.

In our serialized synthesis plan representation, the identity of each product node is determined by the identities of the corresponding reactant molecules. During training, we use a dataset containing full synthesis plans (see Section 4), allowing the correct reaction products to be inserted directly. At test time, however, product information is not available and must be inferred. We therefore employ a reaction-prediction oracle that can predict the most likely product given a set of reactants. The oracle is a transformer-based encoder-decoder model built off the BART architecture [73], closely related to the Molecular Transformer of [74]. Details of its training are provided in Appendix 7.7.

3.1.3 Training objective↩︎

The two decoders are jointly trained with a weighted loss: the 3D decoder uses cross-entropy and MSE losses to reconstruct atom types and coordinates, while the synthesis decoder is optimized with a cross-entropy loss between its output logits and the ground-truth token sequence:

\[\label{eq:loss} \begin{align} \mathcal{L} &= \lambda \mathcal{L}_{\rm 3D} + (1-\lambda)\mathcal{L}_{\text{synthesis}} \\ &\quad + \beta \, D_{\mathrm{KL}}\!\Bigl( \mathcal{N}( \mu_{\mathbf{Z}}, \sigma_{\mathbf{Z}}) \,\Bigm\|\, \mathcal{N}(0,I_d) \Bigr) \end{align}\tag{1}\]

where:

\[\begin{align} \mathcal{L}_{\rm 3D} &= \frac{1}{N} \sum_{i=1}^{N} H(a_i, \hat{a}_i) + \frac{1}{3N} \sum_{i=1}^{N} \frac{\|\tilde{x}_i - \hat{\tilde{x}}_i\|^2}{\sigma^2}, \\ \mathcal{L}_{\text{synthesis}} &= \frac{1}{L} \sum_{i=1}^{L} H(s_i, \hat{s}_i). \end{align}\]

Here, \(\tilde{x}\) represents zero-centered coordinates, \(\sigma^2\) is \(1.0\), \(0\leq \lambda \leq 1\) is the linear interpolation parameter between \(\mathcal{L}_{\text{synthesis}}\) and \(\mathcal{L}_{\rm 3D}\), \(L\) is the total sequence length, and \(s_i\) is the ground-truth synthesis token at step \(i\). Following [75], [29], and [30] we introduce a per-channel KL-penalty weighted by \(\beta\) that regularizes the learned latent toward a standard normal distribution.

3.2 Conditional generative modeling of latent representations↩︎

We train a flow matching model [30], [68], [76], [77] on the latent space learned by our autoencoder, by denoising latent samples from a base Gaussian distribution towards a target distribution represented by latents of ground-truth molecules.

3.2.0.1 Pharmacophore embeddings.

A pharmacophore abstracts components of a molecule into pre-defined key interaction features, providing a compact representation of interaction-relevant moieties and their spatial coordinates. We represent pharmacophores using six discrete pharmacophore types \(\boldsymbol{P}=\{p_i\}_{i=1}^{N_p}\in \mathbb{Z}^{1\times N_p}\) (hydrogen bond donor, hydrogen bond acceptor, cation, anion, aromatic ring and hydrophobe) and coordinates for each pharmacophore feature present in the molecule: \(\boldsymbol{X^{\text{ph}}}=\{x^{\text{ph}}_{i}\}_{i=1}^{N_p}\in \mathbb{R}^{3\times N_p}\). These are extracted from the molecule conformation through a deterministic mapping of the pharmacophore’s atom positions. The features are embedded into a unified learned representation using a standard transformer (see Appendix 7.6), and this embedding serves as conditioning for the diffusion model via cross-attention. The parameters of the pharmacophore transformer encoder and those of the DiT are jointly optimized via the conditional flow matching loss:

\[\begin{gather} \mathcal{L}_{\text{CFM}} = \mathbb{E}_{\substack{ (\boldsymbol{P}, \boldsymbol{X}^{\text{ph}}, \boldsymbol{Z}_1)\sim p_{\text{data}},\\ t \sim \mathcal{U}(0,1),\; \boldsymbol{Z}_0 \sim \mathcal{N}(0,I) }} \Big[ \big\| u_\theta\!\big( \boldsymbol{Z}_t,\; t,\; \tau_\theta(\boldsymbol{P}, \boldsymbol{X}^{\text{ph}}) \big) \\ - v_t\!\big(\boldsymbol{Z}_t \mid \boldsymbol{Z}_0,\boldsymbol{Z}_1\big) \big\|_2^2 \Big], \end{gather}\]

where \(u_\theta\) is the denoiser backbone, which predicts the vector field based on the current state \(z_t\), the time step \(t\), and the pharmacophore conditional embedding \(\tau_\theta(\boldsymbol{P}, \boldsymbol{X^{\text{ph}}})\). \(v_t\) is the ground-truth vector field and the state \(\boldsymbol{Z_t} = (1-t)\boldsymbol{Z_0} + t\boldsymbol{Z_1}\) is a linear interpolation between a clean latent sample \(\boldsymbol{Z_1}\) and a noise sample \(\boldsymbol{Z_0}\) drawn from a standard normal distribution \(\mathcal{N}(0,I)\). The time step \(t\) is sampled uniformly from \(\mathcal{U}(0,1)\). Leveraging the equivalence between the velocity and endpoint flow matching formulations [76], the model receives \(\boldsymbol{Z_t}\) as input and predicts the terminal points of the trajectory during training.

Following [30], we use a Diffusion Transformer (DiT) [69] as our denoiser architecture with self-conditioning [78] and adaptive layer norm for time-step \(t\) modulation. We additionally introduce a cross-attention conditioning mechanism (trained with a conditioning dropout probability of 20%) to sample using classifier-free guidance during inference [79]:

\[\begin{gather} \label{eq:cfg} \hat{u}_\theta(\boldsymbol{Z_t}, t, \tau_\theta(\boldsymbol{P}, \boldsymbol{X^{\text{ph}}})) = (1+w)u_\theta(\boldsymbol{Z_t}, t, \tau_\theta(\boldsymbol{P}, \boldsymbol{X^{\text{ph}}})) \\ - w u_\theta(\boldsymbol{Z_t}, t, \emptyset) \end{gather}\tag{2}\]

where \(\hat{u}_\theta\) is the final vector field prediction. This is a linear combination of the conditional and unconditional predictions, where \(\emptyset\) represents the null token used for the dropped-out condition. The guidance scale \(w\) controls the strength of the conditioning signal and we experiment with both \(w=0\) and \(w>0\).

4 Experiments↩︎

4.0.0.1 Dataset.

We train SynLaD on a set of reaction pathways extracted from the USPTO dataset [80]. We detail the construction of the dataset in Appendix 7.1. Our dataset contains 67,512 synthesis pathways, each with 1 to 6 intermediary reactions. For each final product of the pathways, we compute 10 low-energy conformers using OpenEye’s conformer generation software Omega [81], [82]. Our dataset contains synthesis products with 10 to 40 heavy atoms.

4.0.0.2 Metrics.

We evaluate SynLaD on its ability to generate valid, novel, and realistic molecules. Sampled molecules generated by SynLaD are assessed using validity and internal diversity [83] (definitions provided in Appendix 7.8). Diversity is further measured by the number of unique Bemis-Murcko scaffolds [84], defined as the molecule’s ring systems and connecting linkers with all side chains removed. Synthesizability is additionally evaluated using AiZynthFinder [85], and pose quality is assessed with the PoseBusters benchmark [86] for molecules generated by our 3D decoder. For pharmacophore-conditioned experiments, we use OpenEye’s ROCS scoring functions [5] to evaluate shape and pharmacophore overlap via Tanimoto Shape and Tanimoto Color scores, each ranging from 0 to 1 (with 0 representing no overlap and 1 representing perfect overlap), whose sum defines the Tanimoto Combo score. Following [2], we define a hit as a molecule with a Tanimoto Combo score of at least 1.2 relative to a query.

As described in Section 3.1.2, we use a reaction-prediction oracle to infer product information from sets of reactant tokens sampled by our decoder. The oracle is trained separately; implementation details and accuracy results are provided in Appendix 7.7, where it achieves a top-1 product prediction accuracy of 84.6%.

4.1 Joint conformer and synthesis latent space↩︎

We first investigate reconstruction performance of the joint autoencoder. We randomly select 1000 molecules from the test set, encode them into our latent space, and decode each latent using both decoders. For the 3D decoder, we report match rate and RMSD between ground-truth and predicted coordinates (see Appendix 7.8); for the synthesis decoder, we report match rate and 2D fingerprint Tanimoto similarity between the predicted end product and the ground-truth synthesis outcome. We ablate the latent dimension during training and compare inference strategies for the synthesis decoder (see Appendix 7.4). Results are summarized in Table 1 for both sampling and beam search decoding. Sampling generates \(N\) independent candidate sequences by temperature-scaled stochastic sampling at each step, optionally using top-k truncation by setting logits of tokens outside of the top \(k\) to \(-\infty\). In contrast, beam search deterministically expands the \(b\) (number of beams) most probable partial sequences. In our experiments, we set \(N=1\) and \(k=10\) for sampling and compare against beam search with \(b=5\). Beam search consistently outperforms top-\(k\) sampling and is therefore used for latent-space sampling in subsequent experiments.

2.5pt

4.2 Unconditional generation↩︎

We next train an unconditional latent DiT denoiser using latents produced by our joint VAE (variational autoencoder) described in Section 4.1; model parameters and training details are provided in Appendix 7.5. We sample using \(T=100\) ODE integration steps and decode the resulting latent with both decoders. Results are reported in Table ¿tbl:tab:uncond95results?, where we compare against (1) a latent diffusion model trained without a synthesis decoder head (equivalent to [30]) and (2) an unconditional autoregressive synthesis-decoder-only head. The results show that jointly training the two decoders does not degrade synthesizability or PoseBusters performance, and yields good validity, diversity, and novelty overall. Interestingly, joint training improves the synthesizability of molecules decoded by the 3D decoder, suggesting that the synthesis decoder steers generation toward latent regions more likely to decode to synthesizable molecules. Both SynLaD decoders significantly outperform the synthesis-unconstrained baseline in synthesizability although trained on the same data.

4.2.0.1 Consistency analysis.

To enable pharmacophore-conditioned sampling of synthesizable molecules that also match a desired 3D profile, we evaluate how well the two decoders remain aligned on generated samples. Starting from a latent \(\boldsymbol{Z}\) sampled from the trained diffusion model, we decode both a 3D molecule and a synthesis pathway and quantify how often these outputs correspond to the same, or closely matching, chemical features. We report cross-decoder agreement via exact identity and chemical/shape similarity in Figure 2. The synthesis-decoded products preserve the 3D shape and pharmacophore characteristics implied by the 3D-decoded structures. We show examples of output pairs in Appendix 7.

4.3 Pharmacophore-conditioned molecule generation↩︎

4.3.0.1 In-distribution pharmacophore conditioning.

We evaluate the conditional generation performance of SynLaD by selecting 50 random query molecules from the test set and extracting their pharmacophores. For each query, we generate 100 candidate molecules with SynLaD conditioned on the extracted pharmacophore, and, as a baseline, compare against 5,000 randomly selected molecules from the training set. All candidates—whether sampled from the training set or generated by SynLaD—are evaluated using conformer enumeration with Omega, followed by shape and pharmacophore overlay scoring against the query molecule (see Appendix 7.9). For consistency, we apply the same conformer enumeration to 3D-decoded outputs, even though these samples already exhibit good conformer alignment to the query (see Appendix Figure 10). For each molecule, we retain the highest scoring conformer and assess whether it is a hit or not (i.e., with a Tanimoto combo score \(\ge1.2\) to the query). We report the median number of hits and unique scaffold hits in Figure 3 and distributions of Tanimoto shape and color scores, as well as synthesizability scores assessed using AiZynthFinder (further metrics in Appendix Table 7). Results show that SynLaD generates samples with significantly higher 3D pharmacophore and shape similarities to the query than the baseline, while the synthesis decoder yields molecules with high synthesizability.

4.3.0.2 Screening case study.

To assess SynLaD in a practical drug-discovery setting, we compare it against brute-force ROCS screening over our dataset, which serves as a proxy for a molecular library and reflects standard but time-consuming ligand-based screening. We randomly select 10 query molecules with more than 15 heavy atoms from our test set and calculate ROCS scores against all other molecules in the dataset (\(\sim\)​67k molecules), excluding the query molecules themselves. In parallel, we perform pharmacophore-conditional generation with SynLaD, generating 1000 samples per query. For this experiment, we consider only outputs of the synthesis decoder of SynLaD, ensuring that outputs are synthesizable and directly applicable for experimental testing. We report hit counts for both methods in Table 2, and per-query results in Appendix Figure 8. SynLaD achieves higher hit counts than the typical brute-force baseline in ligand-based screening, while using nearly two orders of magnitude fewer samples (1k molecules) than the size of the screened library (\(\sim\)​67k).

3pt

4.3.0.3 Bioactive hit diversification.

We next investigate the out-of-distribution performance of SynLaD on a bioactive hit diversification task, which aims to find structural analogues of known bioactive compounds. This mirrors real-world drug discovery practices, where diverse candidates are needed to iteratively optimize activity while satisfying synthesizability, developability, and other constraints. We select ten targets from the Lit-PCBA benchmark [35] (see Appendix 8.2) and compare against ShEPhERD [52], SynFormer [25], and REINVENT [37], [87]. ShEPhERD is an SE(3)-equivariant diffusion model that jointly denoises 3D molecular graphs and representations of their shapes and interaction profiles; we compare against their inpainting setting which conditions on electrostatic potential surfaces and pharmacophores. SynFormer is a transformer-based encoder-decoder framework that generates synthesizable 2D molecules in reaction-template and building-block space; we specifically compare against their SynFormerED variants for amortized analogue generation. As an orthogonal baseline, following [88], we fine-tune REINVENT [37] using a ROCS Tanimoto color score as reward function. Unlike SynLaD and the other baselines, REINVENT does not perform amortized sampling and instead requires reinforcement learning-based optimization per query, making the comparison methodologically orthogonal (similar techniques can also be applied to autoencoders [89]). REINVENT is also substantially more computationally expensive, requiring approximately 6 hours per query to generate the optimized samples here, compared to \(\sim\)​1 min for SynLaD, \(\leq\)​1 min for SynFormer, and \(\sim\)​28 min for ShEPhERD (Appendix 8.4) per 100 samples, but we include it as it is a strong baseline in molecular generation.

For each query, we generate 500 molecules by conditioning on the pharmacophore profile of the native ligand in its bound conformation. We apply the same evaluation pipeline as in the previous section and report distributions of aggregated Tanimoto shape and color scores across all ten queries, along with averaged molecule validity, synthesizability, and hit counts, in Figure 4 and Table ¿tbl:tab:pharm95cond95results?. Note that AiZynthFinder is reported for subsets (100 samples due to compute cost) of all generated molecules, not just hits. SynLaD achieves a strong overall balance, producing a high number of synthesizable hits, succeeding on a large fraction of queries, and maintaining high validity. Compared to the other synthesis-constrained baseline, SynFormer, SynLaD’s synthesis decoder produces candidates with consistently higher ROCS shape and color overlap, suggesting that the explicit 3D representation our model learns improves conditional design; it also yields more hits and higher retrosynthesis success. The 3D outputs of SynLaD outperform ShEPhERD on most queries in terms of pharmacophore overlap, hit count, and synthesizability, while achieving comparable shape overlap. The synthesis-decoder output also outperforms ShEPhERD in average number of synthesizable hits.

Although REINVENT achieves the highest average number of synthesizable hits, this result is dominated by a single query with an unusually large hit count; across the benchmark, it generates synthesizable hits for only three queries and shows limited scaffold diversity. We show examples of query-generated molecule pairs with highlighted pharmacophore and shape overlap in Figure 5 and further examples of synthesis and 3D decoder outputs in Appendix Figure 9.

While we report results for both decoders for completeness, the key strength of our approach lies in the synthesis decoder, which—despite operating on linear synthesis trajectories—learns to satisfy shape and pharmacophore constraints more effectively than baselines that do not learn 3D information. This highlights the benefit of our strategy: learning latent representations that encode rich 3D conditioning signals and using them to condition synthesis generation. We further support this claim in Appendix 8.3, where we perform two ablations: 1) we remove the structure-aware latent component (3D reconstruction heads and diffusion module) from SynLaD and keep the synthesis architecture and pharmacophore-conditioning fixed and 2) we switch the 3D reconstruction task for SMILES reconstruction, and keep the rest of the framework identical. We see that both ablations lead to a substantial performance drop, demonstrating that the structure-aware latent component is critical.

5 Conclusion↩︎

In this work, we introduce a conditional latent diffusion framework that decodes latent representations into both 3D molecular structures and reaction-based synthesis pathways. We show that jointly training these two decoders is mutually beneficial at sampling time: it increases the synthesizability of molecules produced by the 3D decoder and improves the synthesis decoder’s ability to satisfy shape/pharmacophore constraints under conditioning, while also enhancing sample diversity compared to a synthesis decoder conditioned directly on pharmacophores. Overall, our results demonstrate that conditioning an autoregressive synthesis decoder through a 3D-informed latent representation yields high numbers of diverse, synthesizable, pharmacophore-aligned hits. Looking ahead, an exciting direction is to scale our method to larger datasets (e.g., Pistachio, [90]), which would enable access to more novel transformations and more combinatorially complex regions of chemical space.

Impact Statement↩︎

This paper presents a generative modeling approach aimed at advancing small-molecule drug discovery by jointly addressing 3D design objectives and synthetic feasibility. By enabling the generation of molecules that are both pharmacophore-consistent and synthesizable, our work has the potential to accelerate early-stage drug discovery and reduce reliance on costly and time-consuming experimental screening pipelines.

At the same time, we acknowledge that generative models for molecular design could, in principle, be misused to propose harmful or toxic compounds. While our work does not explicitly target such applications and is focused on drug discovery settings, the underlying methodology could be adapted in unintended ways.

We emphasize that our approach operates within standard drug discovery pipelines, where generated candidates are subject to extensive downstream validation, including safety, synthesizability, and biological evaluation. As such, we believe the primary impact of this work is to support responsible and beneficial applications in medicinal chemistry.

6 Limitations and Future Work↩︎

The focus of this work is to introduce and validate SynLaD’s ability to generate small molecules that preserve key phramacophore interactions of ligands and have high synthesizability. To this end, we restrict our methodology and experiments to ligand-based pharmacophore conditioning, reaction-constrained synthesis-pathway generation, and evaluation on small-molecule analogue generation tasks, as well as evaluation of our method’s core modules. A current limitation of our framework is that it relies on external reaction predicitons, and its performance is therefore affected by predictor bias or failure models. Future work could evaluate alternative predictors or ensembles of predictors. Second, our experiments here are limited to the USPTO dataset, and scaling to larger and more diverse reaction datasets may further test the generality of the approach. Additionally, further work could also focus on improving cross-decoder agreement (e.g., perhaps via considering techniques from contrastive learning), as stronger consistency across decoding heads may lead to better molecule feature alignment.

7 Methods↩︎

7.1 Dataset↩︎

For our dataset, we make use of the provided synthesis DAGs from [58]. Briefly, these DAGs were obtained from a cleaned version of the USPTO reaction dataset after stripping out reagents [40], [80]; first, by building up a reaction network from the USPTO reactions and an initial set of “building block” nodes (picked by frequency), before extracting a possible synthesis plan for each non-building block node by tracing reactions backward to a loop-free subgraph that terminates in building blocks (selecting a single route when multiple alternatives exist). Using this data, we first filter the reactions such that the final products contain between 10 and 40 heavy atoms, and obtain 67512 synthesis DAGs, 90% of which we keep as our training set and use the remaining data for our validation and test sets. Note that we randomly split by final molecule (i.e., the different train, validation, and test sets lead to different final molecules but may share reactions/sub-networks); analyzing SynLaD’s ability to extrapolate to different data splits (for instance, a “Reachable” vs “Hard” split akin to [91, p. 7]) is an interesting future direction to explore. For each final product of the synthesis DAGs we compute a maximum of 10 low-energy conformers using Omega [81], which results in a total of 584896 3D molecules.

7.2 Reaction network serialization↩︎

A synthesis pathway here is defined as a directed acyclic graph (DAG), where nodes represent molecules (each unique molecule maps to a unique node) and the edges represent reactions. In order to describe these pathways using an autoregressive language model we need to serialize them into a linear sequence of tokens. Our proposed approach follows [58] (alternative serialization schemes have also been proposed, [28], [33], [34]) and is shown in Figure 6. The scheme relies on two kinds of tokens: action tokens and molecular tokens. Action tokens define a new operation to take with respect to the graph and are generally followed by a series of molecular nodes that specify the arguments required to define the action’s effect. There are three types of action tokens:

• ‘B’ adds a building block node to the network (these are easily purchasable compounds and the token is followed by a single molecule token defining the identity of the chosen building block).

• ‘F’ adds a reaction edge (and a product node if the molecule does not yet exist in the network) to our network through a single forward reaction. This action is followed by selecting one or more molecule nodes (representing the existing molecules in the network) to use as reactants, a sub-action ‘F_r’ (which indicates the reaction should be run) and a final product node (followed by ‘P’) that represents the product of the reaction.

• ‘S’ is the stop action, which halts the prediction.

7.3 Synthesis action embeddings↩︎

The input tokens to the transformer are obtained via summing a series of embeddings as shown in Figure 1C: a learned embedding which specifies the token type, a molecule embedding (if the action involves a molecule) which is a projection of the molecule’s 2048-bit Morgan fingerprint, and a learned positional embedding.

7.4 VAE training↩︎

7.4.0.1 Autoencoder ablations

We report ablations for autoencoder hyperparameters in Tables 3 and ¿tbl:tab:app95autoencoder95reconstruct?. Increasing the weight on the synthesis decoder loss (see Eq. 1 ) harms 3D coordinate reconstruction and adds little benefit for the synthesis sequence reconstruction. Note that we report the accuracy of builiding block prediction in its order-invariant form (i.e., building blocks might be added earlier or later on in the sequence, with the end product being the same); this differs during training in which we randomly pick an order (for building blocks at the same level) and train on this fixed sequence using teacher forcing. We use \(\lambda=0.92\) in our experiments and \(\beta=\num{e-5}\), as larger KL weight harms reconstruction accuracy for both decoders. We use a latent dimension \(d=16\) which shows the best molecule and synthesis reconstruction.

2.8pt

1.5pt

7.4.1 Synthesis decoder training↩︎

The synthesis decoder is implemented as a decoder-only transformer. Each input token is embedded by summing (i) a learned embedding for special tokens or retrieved molecular embeddings, (ii) a learned token-type embedding, and (iii) a learned positional embedding. Graph-level information is provided via a fingerprint encoder. The embedded sequence is processed by transformer blocks with RMSNorm pre-normalization, multi-head self-attention, and a SwiGLU feed-forward network (further details can be found in our attached code). When conditioning is enabled, each block additionally applies a cross-attention sublayer over an external conditioning vector. During training, we use teacher forcing [92], i.e., for each prediction step the decoder is conditioned on the ground-truth previous tokens. We compute the cross-entropy loss only on positions where the model is expected to generate tokens, and mask out product tokens (and any padded positions) so they do not contribute to the loss. We report ablations for synthesis decoder inference in Tables ¿tbl:tab:app95beam95search95ablations? and [tbl:tab:app95sampling95ablations].

7.5 DiT training↩︎

7.5.0.1 Sampling ablations.

We investigate the effect of the number of sampling steps used during inference on molecule quality and synthesis/3D decoder output match rate in Table 4. Unless otherwise stated, we use 100 sampling steps in all experiments. We additionally ablate the classifier-free guidance weight \(w\) in Eq. 2 (Table 5). Increasing \(w\) improves control over latent quality and yields more valid decoded molecules, but reduces scaffold diversity. Since our primary objective is to maximize the number of generated hits, we therefore use \(w=0\) (i.e., no classifier-free guidance) throughout.

6pt

6pt

7.6 Training and hyperparameters↩︎

We jointly train the VAE and synthesis decoder, then train the DiT in a second stage, using AdamW with a constant learning rate of \(\num{e-4}\), zero weight decay, and a batch size of 256. We maintain an exponential moving average of the DiT parameters with decay \(0.9999\). Each stage is run to convergence, for up to 3000 epochs or 3 days on a single B200 GPU.

For the VAE, we use standard transformers for both the encoder and 3D decoder with hidden dimension \(d_{\text{model}}=256\), 4 attention heads and 6 layers. For the synthesis decoder we used \(d_{\text{model}}=384\), 8 attention heads, and 8 layers with causal self-attention, cross-attention conditioning to the latents, and SwigLU FFNs.

For the second stage training, the pharmacophore encoder is a standard transformer with \(d_{\text{model}}=256\), 4 heads, and 4 layers, and the DiT denoiser has \(d_{\text{model}}=768\), 12 attention heads, and 12 layers.

7.6.0.1 Molecule and pharmacophore embeddings.

Molecules and pharmacophores are represented by their respective atom or pharmacophore types together with 3D coordinates (see Sections 3.1 and 3.2). We use separate embedding layers for atom / pharmacophore types and atomic positions, following [30]. Atom types are embedded via a learned embedding lookup, while positions are projected using a linear layer. The resulting representations are combined by summation.

7.7 Reaction predictor oracle training↩︎

As explained in the main text, when running our synthesis decoder at inference time we require access to a reaction predictor oracle to predict the molecular identity of the product molecule given the set of reactants (we ignore conditions and reagents in this current work). For this, we use a model based around the BART architecture [73], using an implementation adapted from the transformers library [93] for use in reaction prediction [94]. This is similar in spirit to other encoder-decoder models used for reaction prediction/one-step retrosynthesis [74], [95]–[97]. We train our reaction prediction oracle on reactant-product pairs from the same USPTO dataset we extract the synthesis plans from for training SynLaD; teacher forcing is used at training time and beam search is used during inference for making predictions. We leave to future work the exploration of using alternative reaction predictors [40], [98]–[109], or even multiple at a time [110] for more robust predictions. We report accuracy results from training on a train/validation/test split of approximately 409k/30k/40k data points extracted from the USPTO dataset in Table 6.

4pt

7.8 Evaluation metrics↩︎

We define the metrics used thoroughout the paper as follows:

• Validity (Valid.): % of molecules that can be processed by RDKit.

• Internal diversity (IntDiv.): One minus the value of the averaged pairwise Tanimoto similarities of the Morgan fingerprints (size 2048, radius 2) of the molecules in the generated set.

• Murcko Scaffold count: number of unique Murcko scaffolds in a generated set.

• Novelty (Nov.): % of valid molecules not present in the training set.

• PoseBusters (PB): % of molecules passing all PoseBusters filters [86]. This was only evaluated for models that generated 3D poses.

• SA: Averaged synthetic accessibility score from [111].

• AiZynthFinder (AiZynth.): % of molecules (out of 100) for which AiZynthFinder finds valid synthetic pathways given its default set of reactions and ZINC building blocks.

• RMSD and match rate: For 3D molecule reconstruction we report RMSD and match rate computed using MoleculeMatcher from PyMatGen [112].

• Hit rate: We compute hit rate by counting the number of unique molecules with ROCS Tanimoto Combo score \(\geq 1.2\).

7.9 ROCS evaluation pipeline↩︎

We scored 3D similarity with OpenEye ROCS [5] using a standard OMEGA→ROCS workflow. We kept the query conformation intact (from our test set or Lit-PCBA benchmark pose) and for each candidate molecule, we first generated an ensemble of 100 low-energy 3D conformations with OMEGA [81], which constructs initial geometries from fragment conformers and then samples torsions using a knowledge-based torsion library, retaining a pruned set of energetically and geometrically reasonable conformers. We then used ROCS to perform fast rigid overlays between all query-candidate conformer pairs, representing molecular shape as a sum of atom-centered Gaussian functions and maximizing shared volume and pharmacophoric (“color”) feature overlap. For each candidate, we reported the best overlay across conformers, and computed ShapeTanimoto, ColorTanimoto, and their sum (TanimotoCombo) as the final shape, color, and combined similarity scores. The same pipeline was applied for SynLaD and all baselines, regardless of the methods’ output types.

8 Experiments↩︎

8.1 Pharmacophore conditioned generation↩︎

During training, we load the pharmacophore features for a conformer and randomly drop \(k\) features, where \(k\) is sampled uniformly from \([0, \max(0,N_{\text{ph}}-3)]\), and \(N_{\text{ph}}\) denotes the number of pharmacophore features in the sample. During sampling, we set the latent to contain \(M\) atoms, with \(M\) sampled uniformly from \([\max(1, N-3), N+3]\), where \(N\) is the number of heavy atoms in the query (i.e., excluding hydrogens).

8.2 Bioactive hit diversification experiment↩︎

We used ten PDB targets from [35], with PDB IDs 4LDO, 5L2M, 2P15, 2IOK, 5FV7, 2V3D, 4ZZN, 6B73, 3ZME, 3A2I. We extracted the corresponding ligand for each of these as the ‘query’. We show examples of decoded outputs, together with the conditioning query in Figure 9. We show additional results (including the per-query number of hits for the 3D decoder) in Figure 11.

For the hit diversification experiment, we quantify the extent to which molecules are out-of-distribution from the training set by plotting Tanimoto distances of the Lit-PCBA ligands to all training molecules in Figure 12 (computed using 2048-bit Morgan fingerprints).

In addition to the results in Table ¿tbl:tab:pharm95cond95results?, we provide a ranking table that summarizes the relative performance of each method across all evaluated metrics. For each metric, methods are ranked per task and the average rank is reported, using a joint ranking scheme. This view highlights the consistency of each method across objectives, rather than emphasizing performance on any single metric. As shown in Table 8, SynLaD achieves the best or tied-best average rank for validity, unique scaffold hits, and synthesizable hits.

8.3 Ablations↩︎

8.3.0.1 Removing the 3D VAE and diffusion module.

In principle one could condition a synthesis decoder like the one in SynLaD, or equivalently others from the literature [33], [34] on pharmacophores, to have the model learn to generate synthesizable molecules that preserve those pharmacophores. To see whether there is benefit in modeling the conditioning input to the synthesis decoder using a latent diffusion framework as we do in SynLaD, we perform an ablation where we remove the 3D latent component from SynLaD and condition the synthesis decoder directly on the pharmacophore. We preserve the pharmacophore embedder described in Section 3.2 and use the resulting embedding to condition on with an otherwise identical synthesis decoder via cross-attention. This yields a simplified pipeline, pharmacophore \(\rightarrow\) pharmacophore embedding \(\rightarrow\) synthesis pathway (see Figure 13), which isolates the contribution of the 3D branch while keeping the synthesis decoder architecture and training procedure unchanged. We run the ablation on the in-distribution (ID) and out-of-distribution (OOD) conditioning tasks reported in Section 4 and show in Table 9 that although the Tanimoto Combo scores do not suffer a drastic drop, the ablated version of the model appears to have collapsed to generating few samples per input (see Figure 14), showing that the entire framework is crucial to generating diverse hits.

8.3.0.2 Using a SMILES VAE.

We further investigate the effects of learning a conditional diffusion model for phramacophore-conditioned generation in a less expressive latent space. We replace the 3D autoencoder with a transformer VAE that operates on SMILES strings, to remove the 3D-reasoning component of SynLaD, while keeping the latent-conditioning component. SMILES strings are tokenized at the atom level using a regular-expression tokenizer adapted from [96]. The SMILES VAE preserves the encoder-decoder architecture of the 3D autoencoder, and the same per-token diagonal-Gaussian latent. \(\mathcal{L}_{3D}\) in Eq. 1 is replaced by a cross-entropy loss over tokens, and \(\lambda\) is set to 0.5. We keep an identical synthesis decoder, and train a second-stage latent diffusion model conditioned on pharmacophore embeddings. We run the ablation on both the ID and OOD conditioning tasks in Section 4 and report results in Table 9. We note that the diversity and uniqueness advantages of SynLaD are preserved, supporting our hypothesis that learning a pharmacophore-conditional distribution over latents via the diffusion module helps, promoting variability in the synthesis decoder’s conditioning. In contrast, the number of hits and Tanimoto scores decrease, suggesting that 3D information helps better capture the latent space for the pharmacophore-conditioning task.

3.2pt

8.4 Comparison to baselines↩︎

8.4.0.1 ShEPhERD

We used the official implementation at https://github.com/coleygroup/shepherd [52] and use the \(p(x_1|x_3, x_4)\) conditional setting provided in the repository where \(x_1\) denotes molecular structure, \(x_3\) is the query electrostatic potential surface, and \(x_4\) is the query’s pharmacophore. We generate 500 analogues with number of atoms uniformly sampled from the interval \([\max(3, N-6), N+6]\), where \(N\) is the number of atoms in the query molecule (including hydrogens), to best match SynLaD’s settings.

8.4.0.2 SynFormer

We used the official implementation at https://github.com/wenhao-gao/synformer and changed the following inference settings to allow for higher quality designs compared to the default: search_width=32, exhaustiveness=128, time_limit=300. We generate and evaluate 500 molecules per query.

8.4.0.3 REINVENT

We used the official implementation of REINVENT 4 available at https://github.com/MolecularAI/REINVENT4 [37], [87], [113]. We start from the default REINVENT prior (from https://zenodo.org/records/15641297), which was trained on ChEMBL [114] to cover a broad region of chemical space. We followed the repository’s example reinforcement learning configuration (including hyperparameters), using a Tanimoto color scorer (with up to 5 conformers for computational speed) and a maximum of 300 steps (we do not include any other scoring components). We use the final 500 sampled molecules for evaluation. Although additional steps and hyperparameter tuning may further improve results, each target run already takes approximately 6 hours and can involve scoring up to 30000 proposed molecules; therefore, we defer a comprehensive hyperparameter study (for both REINVENT and SynLaD) to future work.

References↩︎