RGFN: Synthesizable Molecular Generation Using GFlowNets

Authors: michał  koziarski
andrei  rekesh
dmytro  shevchuk
almer  van der sloot
piotr  gai ński
yoshua  bengio
cheng-hao  liu
mike  tyers
robert a batey
Venue: NeurIPS_cc_2024_Conference
Published: 2024-05-15 19:08:36

Abstract: Generative models hold great promise for small molecule discovery, significantly increasing the size of search space compared to traditional in silico screening libraries. However, most existing machine learning methods for small molecule generation suffer from poor synthesizability of candidate compounds, making experimental validation difficult. In this paper we propose Reaction-GFlowNet (RGFN), an extension of the GFlowNet framework that operates directly in the space of chemical reactions, thereby allowing out-of-the-box synthesizability while maintaining comparable quality of generated candidates. We demonstrate that with the proposed set of reactions and building blocks, it is possible to obtain a search space of molecules orders of magnitude larger than existing screening libraries coupled with low cost of synthesis. We also show that the approach scales to very large fragment libraries, further increasing the number of potential molecules. We demonstrate the effectiveness of the proposed approach across a range of oracle models, including pretrained proxy models and GPU-accelerated docking.

Main paper: 
 introduction: 
 In this paper, we propose Reaction-GFlowNet (RGFN), an extension of the GFlowNet framework (Bengio et al., 2023) that generates molecules by combining basic chemical fragments using a chain of reactions. We propose a relatively small collection of cheap and accessible chemical fragments using a chain of reactions, based upon established high-yield chemical transformations, that together still produce a search space orders of magnitude larger than existing chemical libraries. We additionally propose several domain-specific extensions of the GFlowNet framework for state representation and scaling to a larger space of possible actions. We experimentally evaluate RGFN on a set of diverse screening tasks, including docking score approximation with a trained proxy model for soluble epoxide hydrolase (sEH), GPU-accelerated direct docking score calculations for multiple protein targets (Mpro, ClpP, TBLR1 and sEH), and biological activity estimation with a trained proxy model for senolytic classification (Wong et al., 2023) . We demonstrate that RGFN produces similar optimization quality and diversity to existing fragment-based approaches while ensuring synthesizability out-of-the-box. Our analysis further indicates that the generated molecules exhibit a diverse range of chemical properties, as well as target specificity, demonstrating sufficient expressivity of the proposed fragments and synthesis pathways. 
related work: 
 Generative models for molecular discovery. There exists a plethora of methods for molecular generation (Meyers et al., 2021; Bilodeau et al., 2022) using machine learning. They can be categorized depending on the molecular representation used: textual representation such as SMILES (Kang & Cho, 2018; Arús-Pous et al., 2020; Kotsias et al., 2020) , molecular graphs (Jin et al., 2018; Maziarka et al., 2020; Pedawi et al., 2022) or 3D atom coordinate representations (O Pinheiro et al., 2024) ; as well as the underlying methodology, e.g., variational autoencoders (Jin et al., 2018; Maziarka et al., 2020) , reinforcement learning (Pedawi et al., 2022; Korablyov et al., 2024) or diffusion models (Runcie & Mey, 2023) . Recently, Generative Flow Networks (GFlowNets) (Bengio et al., 2021; Nica et al., 2022; Roy et al., 2023; Shen et al., 2023; Volokhova et al., 2024) have emerged as a promising paradigm for molecular generation due to their ability to sample diverse candidate molecules, which is crucial in the drug discovery process. Traditionally, GFlowNets operated on the graph representation level, and candidate molecules were generated as a sequence of actions in which either individual atoms or small molecular fragments were combined to form a final molecule. While using graph representations, as opposed to textual or 3D representations, allows the enforcement of the validity of the generated molecules, it doesn't guarantee a valid route by which to synthesize them. This work expands on the GFlowNet framework by modifying the space of actions to consist of choosing molecular fragments and executing compatible chemical reactions/transformations, in turn guaranteeing both validity and synthesizability.Synthesizability in generative models. One approach to ensuring the synthesizability of generated molecules is by using a scoring function, either utilizing it as one of the optimization criteria (Korablyov et al., 2024) , or as a postprocessing step by which to filter generated molecules. Multiple scoring approaches, both heuristic (Ertl & Schuffenhauer, 2009; Genheden et al., 2020) and ML-based (Liu et al., 2022) , exist in the literature. However, synthesizability estimation is difficult in practice. It can fail to generalize out of distribution in the case of ML models, may significantly reduce the number of high-scoring candidates, and does not necessarily account for the cost of synthesis. Because of this, a preferable approach might be to constrain the space of possible molecules to those easily synthesized by operating in a predefined space of chemical reactions and fragments. Several existing works employ this approach (Bradshaw et al., 2019; Gao et al., 2021; Swanson et al., 2024) , including reinforcement learning-based methods (Gottipati et al., 2020; Horwood & Noutahi, 2020) , which are conceptually closest to this paper. We extend this line of work not only by translating it to the GFlowNet framework but also by proposing a curated set of robust chemical reactions and fragments that ensure efficient synthesis at lower total costs. 
generative flow networks: 
 GFlowNets are amortized variational inference algorithms that are trained to sample from an unnormalized target distribution over compositional objects. GFlowNets aim to sample objects from a set of terminal states X proportionally to a reward function R : X → R + . GFlowNets are defined on a pointed directed acyclic graph (DAG), G = (S, A), where:• s ∈ S are the nodes, referred to as states in our setting, with the special starting state s 0 being the only state with no incoming edges, and the terminal states X have no outgoing edges,• a = s → s ′ ∈ A are the edges, referred to as actions in our setting, and correspond to applying an action while in a state s and landing in state s ′ .A state sequence τ = (s 0 → s 1 → . . . → s n = x), with s n = x ∈ X and a i = (s i → s i+1 ) ∈ A for all i, is called a complete trajectory. We denote the set of trajectories as T .Reaction-GFlowNet generates molecules by combining basic chemical fragments using a chain of reactions. The generation process is illustrated in Figure 1 :In the rest of this section, we describe the design of the Reaction-GFlowNet in detail.Preliminaries. Reaction-GFlowNet uses a predefined set of reaction patterns and molecules introduced in Section 3.3. We denote them as R and M respectively. As a backbone for our forward policy P F , we use a graph transformer model f from (Yun et al., 2019) . The graph transformer takes as an input a molecular graph m and outputs the embedding f (m) ∈ R D , where D is the embedding dimension. In particular, f can embed an empty graph . It can additionally be conditioned on the reaction r ∈ R which we denote as f (m, r). The reaction in this context is represented as its index in the reaction set R.Select an initial fragment. At the beginning of each trajectory, Reaction-GFlowNet selects an initial fragment from M . The probability of choosing i-th fragment m i is equal to:EQUATIONwhere MLP M : R D → R |M | is a multi-layer perceptron (MLP). The σ k is a standard softmax over the logits vector s ∈ R k of the length k:σ k (s) i = exp(s i ) k j=1 exp(s j ).Select the reaction template. The next step is to select a reaction that can be applied to the molecule m. The probability of choosing i-th reaction from R is described as:EQUATIONwhere MLP R : R D → R |R|+1 is an MLP that outputs logits for reactions from R and an additional stop action with index |R| + 1. Choosing the stop action in this phase ends the generation process. Note that not all the reactions may be applied to the molecule m. We appropriately filter those reactions and assume that the score s i for them is equal to -∞.Select another reactant. We want to find a molecule m i ∈ M that will react with m in the reaction r. The probability for selecting m i is defined as:p(m i |m, r) = σ |M | (s) i , s i = MLP M (f (m, r)) i (3)where MLP M is shared with the initial fragment selection phase. As in the previous phase, not all the fragments can be used with the reaction r, so we filter them out.Perform the reaction and select one of the resulting molecules. In this step, we apply the reaction r to the two fragment molecules chosen in previous steps. As the reaction pattern can be matched to multiple parts of the molecules, the result of this operation is a set of possible outcomes M ′ . We choose the molecule m ′ i ∈ M ′ by sampling from the following distribution:EQUATIONwhereMLP M ′ : R D → R scores the embedded m ′ i molecule.Backward Policy. A backward policy in RGFN is only nondeterministic in states corresponding to a molecule m which is a result of performing some reaction r ∈ R on molecule m ′ and reagent m ′′ ∈ R. We denote the set of such tuples (r, m ′ , m ′′ ) that may result in m as T . We override the indexing and let (r i , m ′ i , m ′′ i ) be the i-th tuple from T . The probability of choosing the i-th tuple is:EQUATIONwhere MLP B : R D → R and f is a backbone transformer model similar to the one used in the forward policy. To properly define T , we need to implicitly keep track of the number of reactions performed to obtain m (denoted as k).Only those tuples (r, m ′ , m ′′ ) are contained in the T for which we can recursively obtain m ′ in k -1 reactions.Action Embedding. While the MLP M used to predict the probabilities of selecting a molecule m i ∈ M works well for our predefined M , it underperforms when the size of possible chemical fragments is increased. Our intuition is that such an MLP M struggles to reconstruct the relationship between the molecules. Intuitively, when a molecule m i is chosen in some trajectory, the training signal from the loss function should also influence the probability of choosing a structurally similar m j . However, the MLP M disregards the structural similarity by construction and it intertwines the probabilities of choosing m i and m j only with the softmax function. To incorporate the relationship between molecules into the model, we embed the molecular fragments with a simple machine learning model g and reformulate the probability of choosing a particular fragment m i :EQUATIONwhere ϕ is some activation function (we use GELU) and W ∈ R D×D is a learnable linear layer. Note that if we define g(m i ) as an index embedding function that simply returns a distinct embedding for every m i , we will obtain a formulation equivalent to Equation (3). To leverage the structure of molecules during the training, we use g that linearly embeds a (MACSS) fingerprint of an input molecule m i along with the index i. Note that this approach does not add any additional computational costs during the inference as the embeddings g(m i ) can be cached. In Appendix C.2, we show that this method greatly improves the performance when scaling to larger sets of fragments.Seventeen reactions and 350 building blocks were selected for our first-generation model. The reactions used include amide bond formation, nucleophilic aromatic substitution, Michael addition, isocyanate-based urea synthesis, sulfur fluoride exchange (SuFEx), sulfonyl chloride substitution, alkyne-azide and nitrile-azide cycloadditions, esterification reactions, urea synthesis using carbonyl surrogates, Suzuki-Miyaura, Buchwald-Hartwig, and Sonogashira cross-couplings, amide reduction, and peptide terminal thiourea cyclization reactions to produce iminohydantoins and tetrazoles. The chosen reactions are known to be typically quite robust and generally high-yielding (75-100%), thus enforcing reliable synthesis pathways when sampling molecules from our model. During the construction of the curated building block database, only affordable reagents (building blocks) were considered. For the purposes of this study we define affordable reagents to be those priced at less than or equal to $200 per gram. The mean cost per gram of reagents is $22.52, the lowest cost $0.023 per gram, and the highest cost $190 per gram.A crucial consideration when choosing the set of reactions and fragments used is the state space size. This is difficult to compute precisely since a different set of reactions or fragments is valid for every state in trajectory. We estimate this based on 1,000 random trajectories instead (details can be found in Appendix A). In addition to our 350 low-cost fragments, we also perform this analysis with 8,000 additional random Enamine fragments. Comparison for different numbers of maximum reactions is presented in Figure 2 . As can be seen, even with curated low-cost fragments and limiting the number to a maximum of four reactions, state space size is an order of magnitude greater than Enamine REAL (Enamine, 2024) . This size can increase significantly with the addition of more fragments and/or an increase in the maximum number of reactions. Additional discussion of scaling can be found in Appendix C.2. 350) indicates a variant using 350 hand-picked inexpensive building blocks, while RGFN (8350) also uses 8,000 randomly selected Enamine building blocks. Enamine REAL (6.5B compounds) is shown as a reference. 
experimental study: 
 In the conducted experiments we compare oracle scores and synthesizability scores of RGFN with several state-ofthe-art reference methods. Then, we perform an in-depth examination of produced ligands across several biologically relevant targets. Experimental details can be found in Appendix B. Additional experiments examining the number of modes discovered by different methods and the capabilities of RGFN to scale to larger fragment libraries can be found in Appendix C.We begin experimental evaluation with a comparison to several state-of-the-art methods for molecular discovery. Specifically, we consider a genetic algorithm operating on molecular graphs (GraphGA) (Jensen, 2019) as implemented in (Brown et al., 2019) , which has been demonstrated to be a very strong baseline for molecular discovery (Gao & Coley, 2020) , Monte Carlo tree search-based SyntheMol (Swanson et al., 2024) , and a fragment-based GFlowNet (FGFN) (Bengio et al., 2021) as implemented in (Recursion, 2024) . Training details can be found in Appendix B.6. It is worth noting that besides SyntheMol, which also operates in the space of chemical reactions and building blocks derived from the Enamine database, our remaining benchmarks do not explicitly enforce synthesizability when generating molecules. Because of this, in this section, we will examine not only the quality of generated molecules in terms of optimized properties but also their synthesizability. We consider only a single reactionbased approach, as other existing methods employing this paradigm (Horwood & Noutahi, 2020; Gottipati et al., 2020) do not share code or curated reactions and fragments, making reproduction difficult.We first examine the distributions of rewards found by each method across three different oracles used for training: sEH proxy, senolytic proxy, and GPU-accelerated docking for ClpP. The results are presented in Figure 3 . As can be seen, while RGFN underperforms in terms of average reward when compared to the method not enforcing synthesizability (GraphGA), it outperforms SyntheMol's reaction-based sampling. Interestingly, when compared to standard FGFN, RGFN either performs similarly (ClpP docking) or achieves higher average rewards. This is most striking in the case of the challenging senolytic discovery task, in which a proxy is trained on a severely imbalanced dataset with less than 100 actives, resulting in a sparse reward function. We suspect that this, possibly combined with a lack of compatibility between the FGFN fragments and known senolytics, led to the failure to discover any high-reward molecules. However, RGFN succeeds in the task and finds a wide range of senolytic candidates.Finally, we focus on the synthesizability of generated compounds.We present average values of several synthesizability-related metrics, computed over top-k modes generated for each method, in Table 1 . For completeness, we also include SAScores (Ertl & Schuffenhauer, 2009) , but note that they are only a rough approximation of ease of synthesis. For a better estimate of synthesizability we perform retrosynthesis using AiZynthFinder (Genheden et al., 2020) and count the average number of molecules for which a valid retrosynthesis pathway was found. Note that to reduce variance, we compute SAScores over top-500 modes, but due to high computational cost, AiZynth scores are computed only over top-100 modes. As can be seen, while there is some variance across tasks, RGFN performs similarly to SyntheMol in terms of both synthesizability scores, and significantly outperforms GraphGA and FGFN. All RGFN modes were additionally inspected manually by a chemist and confirmed as synthesizable, which indicates that AiZynth scores are likely underestimated. In the final stage of experiments we examine the capabilities of RGFN to produce high quality binders to a set of diverse protein targets. The aim is to evaluate whether 1) the chemical language used is expressive enough to produce structurally diverse molecules for different targets, and 2) whether generated ligands form realistic poses in the binding pockets. We first demonstrate the diversity of ligands across targets on a UMAP plot of extended-connectivity fingerprints Figure 4 . Ligands assigned to specific targets form very distinct clusters, showcasing their diversity. Interestingly, we observe structural differences between sEH proxy and sEH docking, possibly indicating poor approximation of docking scores by the proxy model. Secondly, we examine the docking poses of the highest scoring generated ligands.As can be seen, the generated molecules produce realistic docking poses, closely resembling those of known ligands (Appendix D). Overall, this demonstrates the usefulness of the proposed approach in the docking-based screens.  
conclusions: 
 In this paper, we present RGFN, an extension of the GFlowNet framework that operates in the action space of chemical reactions. We propose a curated set of high-yield chemical reactions and low-cost molecular fragments that can be used with the method. We demonstrate that even with this small set of reactions and fragments, the proposed approach produces a state space with a size orders of magnitude larger than typical screening libraries while providing high synthesizability of generated compounds. We also show that the size of the search space can be further increased by including additional fragments and that the proposed action embedding mechanism improves scalability to very large fragment spaces.We show that RGFN achieves roughly comparable average rewards to state-of-the-art methods, and it outperforms another approach operating directly in the space of chemical re-actions and, crucially, standard fragment-based GFlowNets. At the same time, it significantly improves the synthesizability of generated compounds when compared to a fragmentbased GFlowNet. Conducted analysis of ligands produced across the set of diverse tasks demonstrates sufficient diversity of proposed chemical space to generalize to various targets. While difficult to demonstrate experimentally, ease of synthesis (due to the small stock of cheap fragments and high-yield chemical reactions used) combined with reasonably high optimization quality offer a promising direction for high-throughput screening applications. 
 
 Appendix: We estimate the state space size by first sampling 1,000 random trajectories, masking out the end-of-sequence action unless the maximum trajectory length max is reached. Then, for every i-th reaction or fragment in the trajectory, we count the average number of valid fragments f rag i and reactions react i from a given state in the trajectory, as well as the average number of unique trajectories traj i into which a state can be decomposed using backward policy. We estimate the state space size asEQUATION)Experimentally derived average values of these parameters can be found in Table 2 . Note that in the second setting we randomly picked 8,000 fragments from the Enamine stock (with the same balancing procedure as in Appendix C.2), which after merging with our own fragments, canonization and duplicate removal yielded a total of 8,317 fragments.Throughout the course of the conducted experimental study, we aim to evaluate the performance of the proposed approach across several diverse biological oracles of interest. This includes proxy models (machine learning oracles, pretrained on the existing data and used for higher computational efficiency): first, the commonly used sEH proxy as described in (Bengio et al., 2021) . Second, a graph neural network trained on the biological activity classification task of senolytic (Wong et al., 2023) recognition. Details of proxy models are provided in Appendix B.2.Per the GFlowNet training algorithm, the reward is calculated for a batch of dozens to hundreds of molecules at each training step, rendering traditional computational docking score algorithms like AutoDock Vina (Trott & Olson, 2010) infeasible for very large training runs. As a result, previous applications of GFlowNets to biological design (Bengio et al., 2021; Shen et al., 2023) employed a fast pre-trained proxy model trained on docking scores instead. These proxies, while lightweight, present potential issues should the GFlowNet generate molecules outside their training data distributions and require receptor-specific datasets. To circumvent this, we use the GPU-accelerated Vina-GPU 2.1 (Tang et al., 2023) implementation of the QuickVina 2 (Alhossary et al., 2015) docking algorithm to calculate docking scores directly in the training loop of RGFN. This approach allows for drastically increased flexibility in protein target selection while eliminating proxy generalization failure. We select human soluble epoxy hydrolase (sEH), ATP-dependent Clp protease proteolytic subunit (ClpP), SARS-CoV-2 main protease (Mpro), and transducin β-like-related protein 1 (5NAF) as targets for evaluating RGFN using a docking reward.Training details for all of the generative methods can be found in Appendix B.6.The sEH proxy is described in (Bengio et al., 2021) . It is an MPNN trained on a normalized docking score data.We utilize the exact same model checkpoint as provided in (Recursion, 2024) .Senolytic classification model is a graph neural network trained on the biological activity classification task of senolytic (Wong et al., 2023) recognition. Specifically, it was trained on two combined, publicly available senolytic datasets (Wong et al., 2023; Smer-Barreto et al., 2023) . Reward is given by the predicted probability of a compound being a senolytic. It is worth noting that due to the low amount of data and high imbalance (< 100 active compounds, a high proportion of which contained macrocycles and were infeasible to construct with fragment-based generative models), this is expected to be a difficult task with sparse reward.The senolytic proxy model consisted of 5 GIN layers (Xu et al., 2018a) with hidden dimensionality of 500, utilized Jumping Knowledge shortcuts (Xu et al., 2018b) , and had a single output MLP layer. Pretraining was done in an unsupervised fashion on the ZINC15 dataset (Sterling & Irwin, 2015) . The training was done for 30 epochs using Adam optimizer with a learning rate of 5 × 10 -5 and batch size of 50.Our docking oracle first accepts canonized SMILES strings as input. These are then converted to RDKit Molecules, protonated, and a low-energy conformer is generated and minimized with the ETKDG (Riniker & Landrum, 2015) conformer generation method and UFF (Rappe et al., 1992) force field, respectively. For computational efficiency, we generate one initial conformer per ligand. Each conformer is converted to a pdbqt file and docked against a target with Vina-GPU 2.1 using model defaults: exhaustiveness (denoted by "thread" in the implementation) of 8000 and a heuristically determined search depth d given byEQUATIONwhere N atom and N rot are the number of atoms and the number of rotatable bonds, respectively, in the generated molecule. Box sizes were determined individually to encompass each target binding site and centroids were calculated to be the average position of ligand atoms in the receptor PDB. A negative score is calculated and returned as a reward.Each target was prepared by removing its complexed inhibitor and atoms of other solvent or solute molecules. We selectively prepared the ClpP 7UVU protein structure by retaining only two monomeric units to ensure the presence of a single active site available for ligand binding and similarly prepared the Mpro 6W63 protein structure by retaining only one monomeric unit.To ensure the diversity, specificity, and conformer validity of top-generated molecules for each target, we initially categorized our molecules into distinct modes, each representing any SMILES string with a Tanimoto similarity of 0.5 or lower with all other modes. Subsequently, we selected the top 100 modes based on their Vina-GPU 2.1 scores and filtered their docked poses using PoseBusters (Buttenschoen et al., 2023) in "mol" mode, where any pose failing any PoseBusters check was excluded from consideration. As a final precaution, we selected only modes with Tanimoto coefficients to known aggregators of 0.4 or lower using UCSF's Aggregation Advisor (Irwin et al., 2015) dataset. This process resulted in 35, 68, 31, and 15 top modes for sEH, ClpP, Mpro, and TBLR1 binders, respectively (Appendix E). Comparative analyses of docked top RGFN modes and confirmed sEH, ClpP, and Mpro ligand poses can be found in Appendix D. TBLR1 was omitted from the analysis due to a lack of known small-molecule ligands.Both RGFN and FGFN were trained with trajectory balance loss (Malkin et al., 2022) using Adam optimizer with a learning rate of 1×10 -3 , logZ learning rate of 1×10 -1 , and batch size of 100. The training lasted 4,000 steps. A random action probability of 0.05 and a replay buffer of 20 samples per batch were used. Both methods use graph transformer policy with 5 layers, 4 heads, and 64 hidden dimensions. Exponentiated reward R(x) = exp(β * score(x)) was used, with β dependent on the task: 8 for sEH proxy, 0.5 for senolytic proxy, and 4 for all docking runs. Note that due to different ranges of score values, this resulted in a roughly comparable range of reward values.All sampling algorithms were outfitted with the Vina GPU-2.1 docking, senolytic proxy, and sEH proxy scoring func-tions. While model architecture hyperparameters and batch sizes were kept consistent between FGFN and RGFN, we allowed FGFN a maximum fragment count of 6 as opposed to RGFN's 5 due to RGFN's larger average building block sizes.GuacaMol's Graph GA model was trained with a population size of 100, offspring size of 200, and a mutation rate of 0.01 for 2000 generations for a total of 400,000 visited molecules.For SyntheMol experiments, we used the default building block library of 132,479 compatible molecules and precomputed docking, senolytic, and sEH proxy scores for all prior to executing rollouts to follow the established methodology. Due to CPU constraints, sampling 500,000 molecules with SyntheMol was impractical. Instead, we executed 100,000 rollouts over approximately 72 hours to match the RGFN training time with docking, yielding 111,964 unique molecules. Additionally, we performed 50,000 rollouts each (approximately 30 hours) for sEH and senolytic proxies, resulting in 73,941 and 69,652 unique molecules, respectively.We examine the number of discovered modes for each method, with a mode defined as a molecule with computed reward above a threshold (sEH: 7, senolytics: 50, ClpP docking: 10), and Tanimoto similarity to every other mode < 0.5. The number of discovered modes across tasks as a function of normalized iterations is presented in Appendix C.1. Note that in the case of GraphGA, FGFN, and RGFN this simply translates to the number of oracle calls, but for Synthe-Mol, due to large computational overhead, we impose a maximum number of rollouts such that training time was comparable to RGFN (see Appendix B.6 for details). As can be seen, despite slightly worse average rewards, FGFN still outperforms other methods in terms of the number of discovered modes (with the exception of senolytic discovery task, where it fails to discover any high-reward molecules). This suggests that RGFN samples are less diverse, possibly due to the relatively small number of fragments and reactions used. However, RGFN still outperforms remaining methods across all tasks, suggesting that it preserves some of the benefits of the diversity-focused GFlowNet framework.Next we investigate the influence of a fragment embedding scheme proposed in Section 3.2. In the standard implementation of the GFlowNet policy, actions are represented as independent embeddings in the MLP. These encode actions as indices, effectively disregarding their respective internal structures and all information contained therein. If there is a helpful structure within the actions, the independent embeddings will need to learn it from scratch. While this may be a relatively easy task for small action spaces, it becomes more difficult when the size of the action space increases. To scale RGFN to a larger size of the building block library, we proposed to encode fragment actions using molecular fingerprints, allowing the model to leverage the internal structure of the actions without any additional computational overhead during the inference. In Figure 6 , we observe that our fingerprint embedding scheme allows for drastically faster convergence compared to the standard independent action embedding, especially for large library sizes. The number of discovered scaffolds with sEH proxy value above 7 (a) and 8 (b) as a function of fragment library size. We compare standard independent embeddings of fragment selection actions (blue) with our fingerprint-based embeddings (orange) that account for the fragments' chemical structure. The number of scaffolds is reported after 2k training iterations for 3 random seeds (the solid line is the median, while the shaded area spans from minimum to maximum values). We observe that our approach greatly outperforms independent embedding when scaling to a larger action space. 


Decision: This paper studies generative models for molecular discovery, focusing on incorporating synthesizability by generating molecules within the chemical reaction space. The proposed method, RGFN, performs molecular generation explicitly in this space.

Given that this paper is on the borderline, we had a detailed AC-reviewer discussion to carefully assess pros and cons of the manuscript. The reviewers mainly identified the following concerns:
- **Evaluation on the AiZynth Score:** The evaluation is not entirely convincing. Specifically, it is important to clarify whether the reaction templates used in the GFlowNet framework overlap with those in AiZynthFinder and whether the baseline models are appropriately tuned, as these factors could significantly influence the scores.
- **Novelty:** The novelty of the approach is somewhat limited, as it is a relatively straightforward extension of the existing GFlowNet method.
- **Discussion of Related Work:** The initial submission lacked a sufficient discussion of related work, although this was adequately addressed in the authors' rebuttal.

In my opinion, these issues are not prohibitive for publication. For a new method, especially in empirical research, comprehensive evaluation is crucial, while it is not mandatory for the proposed method to outperform all existing approaches. Rather, it is essential that the method is reliable, which I believe is an important contribution from the ML side, and this paper satisfies it.

Considering these points, I recommend accepting the paper. Please ensure that the final version includes revisions based on the reviewers' suggestions and incorporates the additional empirical results provided in the authors' rebuttal, particularly those in the PDF.

Reviews: review: paper_summary: paper_summary: This work proposes a workflow to synthesize molecules with reaction templates starting from some building blocks. Through this pipeline, the synthesizability of the generated molecules can be improved based on some experimental evidence of molecular drug design tasks.
strength_weakness: strength_weakness: (1) The related work discussion in this work is insufficient. This reviewer just searches a bit on the web and find three highly related articles, which are "Bridging the gap between chemical reaction pretraining and conditional molecule generation with a unified model", "A generative model for molecule generation based on chemical reaction trees" and "Generating molecules via chemical reactions". This reviewer believes other relevant articles are not covered in this article. A more comprehensive and systematic related work section is required;

(2) The proposed method has many limitations, including reliance on reaction templates and powerful pre-trained models. We state these limitations in the following section;

(3) The evaluation of the proposed approach is not comprehensive enough. It does not include the conventional molecular generation benchmarks like novelty, validity, and diversity.
questions: questions: (1) This reviewer is concerned about the pre-trained oracle models. It seems that almost every step of the pipeline relies on a single pre-trained model. So how to derive these pre-trained models independently?

(2) This reviewer is concerned about the ability to generate truly novel molecules. It seems the proposed approach starts with a selected set of fragments and deduces products through reaction templates. However, since reaction templates are extracted from known reactions, it seems nearly impossible to generate purely novel molecules through this pipeline. How do the authors address and evaluate this issue?
limitations: limitations: (1) It seems the whole pipeline of RGFN heavily relies on many different pre-trained models. Each stage can affect the overall process drastically. 

(2) The proposed approach is based on reaction templates. The reaction templates must be updated frequently to cover more types of reactions. However, the template-based method has poor generalization to unseen reactions (new combinations of reactants).
paper_summary: paper_summary: This work proposes a GFlowNet-based framework for synthesizable molecule design, which comprises of building block selection stages and reaction selection stages. The GFlowNet-based model uses chemical oracle functions as reward and the goal is to generate synthesizable molecules with the desired property quantified by the scoring function.
strength_weakness: strength_weakness: - As acknowledged in the paper, the major limitation is scalability to a larger chemical space. It becomes much more difficult to discover hit scaffolds when the chemical space gets larger. In addition, the capability of this framework is limited by oracle functions.
questions: questions: - Molecules generated by building blocks and reactions are guaranteed to be highly likely synthesizable. Is it still necessary to evaluate empirical synthetic accessibility scores of these molecules? as these scores are generally inaccurate and they cannot tell much when a molecule already has synthesizability guarantee.
- Can this framework be applied to bottom-up synthesis planning?
limitations: limitations: See Weaknesses
paper_summary: paper_summary: This paper presents a model for molecule design called Reaction-GFlowNet (RGFN). RGFN is an extension of the GFlowNet framework [4] (which has previously been used to generate molecules by building them out of small fragments) to generate molecules through virtual chemical reactions, ensuring that the molecules RGFN generates are more likely to be synthesizable. The paper evaluates RGFN against GraphGA (a graph genetic algorithm), FGFN (the original fragment-based GFlowNet), and SyntheMol (a recent synthesis-based generative model of molecules proposed in [54] based around Monte Carlo tree search). The paper shows RGFN finds competitively scoring molecules on two proxy* and one docking task (as expected the unconstrained GraphGA does best), while also ensuring the molecules it suggests do better on important synthesizability metrics.

* proxy meaning the task is to optimize against a graph neural network (GNN) oracle trained on a small amount of target data.

(Edited Aug 14 to update the score: see comments below).
strength_weakness: strength_weakness: ## Originality
As mentioned under the originality/significance heading in the section above, RGFN seems a fairly straightforward combination of the ideas of GFlowNets with recent synthesis-based generative models of molecules. It would have been helpful if there had been more focus on why synthesis-based GFlowNets were a better approach to take than the currently proposed alternatives to better promote the paper's significance. (Hence why I have gone with a low contribution score).

## Only one other synthesis-based generative model of molecules is compared against 
The only other synthesis-based generative model for molecules that is compared against in the experiments is SyntheMol, and this baseline has been limited here in the amount of compute it uses. I can understand no comparison is done against [22,19] due to the lack of open-sourced code, but it would have been nice to have had a comparison against some of the other ones cited, e.g., [17] with its open sourced code available on GitHub (https://github.com/wenhao-gao/SynNet). It would also be nice to have a comparison to the other mentioned methods for encouraging synthesizability, e.g. via using scoring methods [34].

(This along with the fact that no code is available -- see below -- is why I have gone with a lower soundness score).

## No code is shared
As far as I'm aware there is no code currently available for RGFN? It does not seem to be provided as part of the submission (despite this box being ticked in the paper checklist)? I would find it impossible to reproduce the results of the paper using the details provided in Appendices B-D alone (e.g., I'm unsure of what building blocks and reaction templates to use or even how the molecules fed into the GNNs were featurized). One of the proxies was initially trained in an unsupervised manner on the ZINC dataset (Appendix B), but there are scant details on how this was done.

On a similar note, the compute resources used to train the models also seems to be missing. This is again marked as included in the paper checklist, so would appreciate being pointed towards this if I have made a mistake and missed this?    


## Method is currently limited to a restricted synthetic space
RGFN uses only 350 building blocks and 17 reactions. This is low compared to the ~150k building blocks and 49 reactions used in [19], or the 5000 building blocks and 90 reactions used in [22]. The experiments done in Section 4.3 suggests the RGFN scales badly to larger building block libraries even with architecture modifications and larger reaction libraries do not seem to have been tried. This, combined with the fact that RFGN can only explore linear synthetic trees (i.e., multiple complex intermediates cannot be created and then combined) and bimolecular reactions restricts the molecules that the approach could explore. This might be an issue in more challenging optimization tasks.
questions: questions: 1. What is the distribution of the number of reaction steps typically used when finding useful molecules?



2. I don't completely follow how the size of the state space was calculated in Figure 2. Perhaps a cartoon would be useful in Appendix A. Is the state space size different to the total number of unique molecules that can be generated?



3. I'm confused as to why RGFN seems to often outperform FGFN in the docking/proxy scores of the molecules it finds. I would have thought the latter model is more flexible and so should be able to find better scoring molecules, even if they did not end up being synthesizable. The paper suggests that the library of fragments available for RGFN might be less suitable than those found in the building blocks picked for RGFN. Have different fragments been tried or other experiments carried out to investigate this hypothesis?



4. RGFN seems to produce molecules with lower QED scores than SyntheMol and FGFN: is there any intuition for why this might be the case?



5. Line 243 says:
> "All RGFN modes were additionally inspected manually by a chemist and confirmed as synthesizable, which indicates that AiZynth scores are likely underestimated."  
   
Does this mean that a chemist took the proposed molecule and synthesized it in the lab? Were the proposed synthetic routes from RGFN used or were alternatives derived?
 


6. Line 282 states:
>  "It is also important to recognize that RGFN does not explicitly generate synthetic routes to the molecules."    

I did not understand what was meant by this. I thought the whole point of RGFN is that it did explicitly generate synthetic routes to the final molecule through the sampled trajectory? 



7. Equation 4 versus 7:
a. My understanding is that the architecture expressed in Equation 4 is used instead of that of Equation 7 (apart from in Section 4.3)?
b. The architecture expressed in Equation 7 seems to do better in Figure 5 so why not always use this?
c. Could also a similar principle be used to scale up the number of reaction templates available?



8. Line 221 discusses a procedure for picking molecular modes from the list of proposed molecules. Given modes have to be a certain Tanimoto distance apart, I assume this procedure is run greedily and is somewhat dependent on how the molecules are ordered?
limitations: limitations: Two limitations are discussed in Section 5: (i) that the model is currently limited to a very small reaction and building block library, and (ii) the issues inherent with docking oracles. Regarding the first of these, the authors discuss and evaluate a modification of the architecture to try to enable the method to scale (see Section 4.3), but performance still degrades as the initial fragment library increases and so this seems to remain an outstanding limitation (see also the weaknesses section above). The second limitation (docking oracles are not perfect) is not limited to RGFN and is a problem among similar models more generally. The authors discuss incorporating additional, more-accurate oracles in a multi-fidelity framework to resolve this which seems reasonable.

Another limitation which is not discussed is the sample efficiency of the proposed approach (which is again also a limitation for the baselines). The number of molecules the non-GFlowNet baselines visit during optimization for the different tasks is described in Section 4. These are high and range from 70k-400k, and I would be interested in how this compares to the two GFlowNet models.
paper_summary: paper_summary: This study modifies prior GFlowNet-based frameworks to produce molecules satisfying synthesizability. While previous GFlowNet-based molecular generation completes molecules by adding fragments or atoms, the proposed method repeats following stems to generate molecules: (1) selects reaction templates, (2) selects reactants, and (3) performs reactions and selects one of the resulting molecules. In the experiments, the proposed method shows competitive performance while preserving synthesizability.
strength_weakness: strength_weakness: - The preliminary is insufficient for understanding the components of the method, particularly mentioning the flow matching condition while the method section actually describes forward and backward policy.
- To improve synthesizability-related metrics, one may consider including them in the rewards. Can authors provide a comparison with FGFN including the SA score in the reward? (e.g., multiobjective GFN)
- The lower number of modes in RGFN (compared to FGFN) is due to constraints in the generation space (ensure high SA scores). Therefore, it would be interesting to compare the number of modes filtered by SA scores.
- Can authors compare the cost of generation time? Additionally, I am curious about time costs for (finding a synthesizability path + generation with FGFN) vs. (generation with RGFN).

**Minor:**
- Error in line 114 $m'i$
questions: questions: See weaknesses.
limitations: limitations: The authors have acknowledged the limitations of their approach.

Comments: title: Thanks for discussion and further results
Requested additional experiments
Thanks for your reply: responding to your explanation and request for me to clarify my point
Reply to Authors
comment: Thank you for following up. Yes, while I agree your curated set of reactions and fragments could be an interesting contribution, currently it is hard for me to judge given that they are not actually included in the submission (related to point W3). I believe (from the initial rebuttal) you have provided this to the AC though (note I cannot see the AC comment) so leave it with them.

Anyway, I have edited my score to reflect many of my original points getting addressed and to reflect the additional experiments carried out (including the one recently posted), showing that the diversity of the produced molecules could be a reason to prefer the approach here over other previously proposed synthesis-based de novo design methods.

Thanks also for the discussion!
We wanted to provide updated results for the experiments requested by Reviewers (added DRD2 task, and two additional baselines: RxnGenerator [1] and FGFN with the SA score included as a reward term). By now, all of the runs except RxnGenerator in the ClpP task were completed. Crucially, similar to SyntheMol, due to very large computational overhead we capped the run time of RxnGenerator at approximately 24 hours per run (which is above the total RGFN run time). For the final paper, we intend to actually further increase the run time of RxnGenerator, but based on the current results (smaller number of molecules visited and significantly lower average rewards), we believe any changes in the conclusions to be very unlikely. Since we are not allowed to provide any figures at this point, we instead show a table with two most relevant metrics (AiZynthFinder scores and the number of discovered modes).

In general, the trends are comparable to the ones observed in the previous experiments: RGFN shows synthesizability scores comparable with other methods providing synthesizability guarantees (SyntheMol and RxnGenerator), illustrating the benefits of amortization of the GFlowNet framework. While in the sEH task AiZynth scores of RGFN are slightly lower, we suspect that it is mostly due to higher average weights of the generated molecules (since larger molecules are more difficult for AiZynthFinder to predict paths for). At the same time RGFN provides a significantly larger number of discovered modes (lower only than FGFN). Crucially, including SA score as a reward term for FGFN only marginally improves synthesizability (due to the unreliability of SA scores), at the same time decreasing the number of discovered modes. **Overall, this illustrates the usefulness of the RGFN framework, which has synthesizability guarantees similar to other synthesizability-focused methods, while producing significantly more diverse candidate molecules** (and does so using cheaper building blocks).

| Task  | Method       | Mol. weight $\downarrow$ | AiZynth $\uparrow$ | Modes $\uparrow$ |
|-------|--------------|--------------------------|--------------------|------------------|
| DRD2  | GraphGA      | 475.4 ± 53.2             | 0.41               | 664              |
|       | SyntheMol    | 365.6 ± 54.3             | 0.66               | 728              |
|       | RxnGenerator | 418.1 ± 53.3             | **0.87**           | 36               |
|       | FGFN         | 386.5 ± 45.0             | 0.76               | 11396            |
|       | FGFN+SA      | 381.1 ± 35.1             | 0.78               | 5903             |
|       | RGFN         | 447.1 ± 45.7             | **0.87**           | 6076             |
| sEH   | GraphGA      | 528.6 ± 42.3             | 0.04               | 1327             |
|       | SyntheMol    | 411.1 ± 66.7             | *0.80*             | 2336             |
|       | RxnGenerator | 421.6 ± 103.4            | **0.82**           | 5                |
|       | FGFN         | 473.4 ± 58.9             | 0.14               | 21548            |
|       | FGFN+SA      | 473.7 ± 62.2             | 0.27               | 5796             |
|       | RGFN         | 495.2 ± 49.6             | 0.56               | 2875             |
| Seno. | GraphGA      | 485.7 ± 75.6             | 0.05               | 347              |
|       | SyntheMol    | 441.4 ± 83.5             | 0.53               | 598              |
|       | RxnGenerator | 431.5 ± 100.9            | **0.65**           | 12               |
|       | FGFN         | 468.9 ± 47.7             | 0.02               | 44               |
|       | FGFN+SA      | 451.8 ± 54.5             | 0.13               | 14               |
|       | RGFN         | 558.7 ± 62.8             | *0.58*             | 2276             |
| ClpP  | GraphGA      | 521.0 ± 31.8             | 0.00               | 1571             |
|       | SyntheMol    | 458.2 ± 60.7             | *0.56*             | 1107             |
|       | RxnGenerator | -                        | -                  | -                |
|       | FGFN         | 548.6 ± 42.9             | 0.25               | 16960            |
|       | FGFN+SA      | 509.2 ± 52.4             | 0.33               | 11959            |
|       | RGFN         | 526.2 ± 37.6             | **0.65**           | 2880             |

[1] Nguyen, Dai Hai, and Koji Tsuda. "Generating reaction trees with cascaded variational autoencoders." The Journal of Chemical Physics 156.4 (2022).
Thank you for your clarification! To address your points: yes, we agree that some of the existing methods (including the ones you mentioned) could potentially be extended to use e.g. our predefined set of high-yield chemical reactions and low-cost fragments. We would argue that in some sense, this makes our contribution of proposing a curated set of reactions and fragments more impactful, as it can benefit development of other methods. Moreover, even if it might seem “straightforward”, it required significant effort and chemical expertise to craft - so in a sense, it’s perhaps a slightly more interdisciplinary paper than a typical NeurIPS submission (but we would like to think that it makes it potentially more useful to propagate within the ML community).

Regarding the dataset argument, we agree to some extent, but would like to point out that even such a large synthetic dataset would not guarantee generalization out of distribution. The ability to enforce hard domain constraints (by masking invalid actions) is one of the strengths of our method that makes it useful in this context. Finally, it's worth mentioning that the authors of ChemBO paper recognize this limitation:

> While outputs of reactions are well known for simple cases, it is impossible to predict outcomes with complex molecules, and in some cases, the outputs may not even be deterministic. Fortunately however, there have been several advances in computational chemistry to predict outcomes of chemical reactions, which can be used in place of the oracle. In our work we use Rexgen [26]. It should be emphasized that since such predictors are not perfect, so in practice, ChemBO could end up recommending unsynthesizable molecules and/or incorrect synthesis recipes

Thank you once again for an insightful discussion!
We thank the Reviewer for answering our rebuttal.

Regarding points 1 and 2, we again apologize for the omission of some mentioned papers and we are keen on improving the writing as per the Reviewer’s valuable suggestions. However, we believe that it does not fundamentally alter the presentation of our main novel contributions, which as discussed in the reply to the Reviewer 5soV are:
- Handling multiple possible products of a given reaction, which results in a de facto non-deterministic environment. This is not an issue in FGFN, as operating on graphs allows us to specify the place of inserting a new fragment. We deal with this by introducing a product selection step (Eq. 5).
- Since GFNs require computation of parent states and masking invalid actions, this required implementing efficient recursive decomposition of molecules into BBs and crafting specific SMARTS templates for reactions that would limit the possible number of parents.
- Improving scalability by introducing fingerprint-based action embeddings for fragment selection.
- Crafting a specific set of high-yield reactions and low-cost fragments that would enable fast and cheap synthesis experimentally, which has been a long-standing challenge in de novo small molecule discovery.

Regarding point 3, we would like to emphasize that based on the (noisy) synthesizability metrics, RGFN does in fact achieve state-of-the-art performance in terms of synthesizability (and cost of synthesis, as illustrated in the Appendix), while preserving high performance in terms of mode discovery (higher than baseline methods other than FGFN). It’s this favorable balance that makes RGFN, in our opinion, an impactful approach, as easy synthesizability is the most important practical factor (unsynthesizable molecules with high property scores are not useful experimentally). As discussed in the paper, QED, molecular weight are included for completeness, but are secondary to AiZynthFinder scores, which in our opinions are the best approximation of synthesizability currently available. We will further clarify this in the updated manuscript with the Reviewer’s comments in mind.

As discussed in the first answer, we would like to clarify that the metrics we utilize here are **chemically more rigorous than the typical VUN metrics**. VUN metrics are more suitable for methods such as 3D diffusion-based generative modeling where the model is trained on a dataset and inference can generate invalid molecules. All molecules generated by RGFN are valid by design (validity), RGFN is not trained on a dataset of molecules but on a reward function (i.e. novelty is not directly applicable, in the sense all of the molecules are “new”), and we use diversity which is a much stricter version of uniqueness. Typical graph-based generation methods, including most GFlowNet papers, use these more strict metrics as adopted here ([1-6]). We would argue that the ones used in our paper are more appropriate for illustrating a good trade-off between synthesizability and diversity of generated molecules.

Regarding the comment that “additional reasonable metric for evaluating synthesizability should be incorporated”, we are eager to improve our manuscript and would kindly ask for the Reviewer to provide some specific suggestions for what other chemically rigorous metrics could be used. We were unable to identify realistic approaches to reliably do this - lack of good synthesizability evaluation metrics is, in fact, the main motivation behind RGFN. As pointed out by Reviewer rpKJ, showing the synthesizability of RGFN might not even be strictly necessary, as it guarantees it out-of-the-box by operating directly in the space of high-yield chemical reactions, and synthesizability metrics are inherently noisy.

We hope you were convinced by our rebuttal to the other points in the first review (e.g. pre-trained models, usage of reaction templates). Please, let us know if you have any questions.

Thank you once again for your detailed review and for considering our revisions. We are active in any further discussions and would really appreciate it if AC and SAC could take our explanations and clarification into consideration!

[1] Zhang, Zaixi, et al. "Molecule generation for target protein binding with structural motifs." ICLR (2023).  
[2] Guo, Jeff, et al. "Link-INVENT: generative linker design with reinforcement learning." Digital Discovery (2023).  
[3] Zhu, Yiheng, et al. "Sample-efficient multi-objective molecular optimization with GFlowNets." NeurIPS (2024).  
[4] Korovina, Ksenia, et al. "ChemBO: Bayesian optimization of small organic molecules with synthesizable recommendations." PMLR (2020).  
[5] Nguyen, Dai Hai, and Koji Tsuda. "Generating reaction trees with cascaded variational autoencoders." The Journal of Chemical Physics (2022).  
[6] Swanson, Kyle, et al. "Generative AI for designing and validating easily synthesizable and structurally novel antibiotics." Nature Machine Intelligence (2024).
Thank you for your response and kind words!

**Re Q6:** This is helpful thanks and answers my question!


**Re my remark:**
Happy to clarify this point. This was in response to the global "Author Rebuttal by Authors" posted above, namely under the section "Originality / Similar Works" and the statements:
> We note that the majority (UniRXN, RXNGenerator, SynNet, ChemBO) rely on pre-trained neural networks to predict synthetic trees and therefore require the curation of large datasets of druglike molecules and generating ground-truth synthesis routes. ... . Moreover, while template-free methods like UniRXN, ChemBO, and MoleculeChef are capable of accurately predicting the products of reactions given reactants, they do not specify exactly the reaction used at each step - valuable information for chemists to prioritize easy-to-synthesize molecules out of a large pool of generated candidates. 

The particular form of the reaction oracle does not seem to be a key component of many current synthesis-based de novo design algorithms. To use methods which initially ran with template-free reaction oracles, such as ChemBO and MoleculeChef, with a template-based reaction oracle (that specified exactly the reaction used at each step) would be a very straightforward switch.

Moreover, I don't believe that the requirement to curate large datasets is a particularly strong argument either. Yes, a model like SynNet uses approximately 200k synthesis plans for pre-training its NNs, but these synthesis plans are generated artificially and so this step seems on par with the 400k synthetic routes generated by GFlowNet in the process of its optimization. Likewise, ​​ChemBO does not actually require pre-training on synthetic routes.

Regarding the significance/originality, I was somewhat more convinced by the arguments you made in the rebuttal to my review above, suggesting work was needed to develop GFlowNets to operate in reaction space that I might have overlooked (bullet points: 1, 2, and 4, starting with "Handling multiple possible products...").
Thank you for providing a detailed rebuttal. Below are our reviewer's responses:

(1) The reviewer expresses concern regarding the omission of related works, which may impact the novelty and significance of the proposed research. Merely including the omitted literature as references is insufficient. Further writing improvements are necessary to elucidate the major contributions of the proposed method compared to previous methods, thereby justifying its significance.

(2) From a machine learning algorithm perspective, the technical contribution of the proposed method appears somewhat limited. The algorithm generates molecules using reaction templates, an approach widely adopted in template-based retrosynthesis analysis and forward reaction prediction. The reviewer does not identify significant algorithmic differences between the proposed method and template-based reaction modeling (both forward and backward).

(3) The evaluations appear somewhat confusing and not comprehensive enough. Table 1 suggests that the proposed RGFN method does not achieve state-of-the-art performance in most tasks, causing the reviewer to question the method's superiority. Moreover, since the title emphasizes "molecular generation," typical empirical comparisons should include molecular generation benchmarks to demonstrate the basic VUN (Validity, Uniqueness, Novelty) of the proposed method. Building upon these fundamentals, an additional reasonable metric for evaluating synthesizability should be incorporated. Although the reviewer genuinely believes that molecular generations based on chemical reactions can be more reliable and synthesizable, the current evaluations are insufficient to demonstrate RetroGFN's superiority.

In conclusion, the reviewer has decided to maintain their rating, and we kindly defer the decision to the AC.
Thank you for your comment. We would also like to take this opportunity to thank you again for a very thorough review, as that didn’t fit in the official answer. It is very appreciated and definitely helped us improve the paper.

Answering the follow up question regarding Q6:

Let us preface by saying that we suspect some of the ambiguity comes from different common understanding of the term “synthesis route” as used in ML and chemistry communities. In the ML community, the synthesis route is a sequence of reactions transforming sets of reactants into products. According to this definition, RGFN framework can output a synthesis route (albeit not necessarily the optimal one). However, in the answers and the paper we refer to the chemistry-based understanding of the term, elaborated on below. We apologize for the lack of clarity, we tried to be precise in chemical terminology, but recognize that given the audience we should clarify this, which we will do in the revised paper.

To elaborate on the additional components needed for actual wet lab synthesis (which we refer to when we discuss a “synthesis route”):

An explicit route is not generated, but rather the likely building blocks necessary to construct a molecule (and the knowledge of the reactions encoded) will be readily identifiable. This does not mean that an explicit synthetic route has been generated. A synthesis is more than just the building blocks and reactions used.

1. Firstly, consider a molecule ABCD constructed from four fragments (building blocks). A synthesis of such a molecule might proceed as A + B -> **AB** + C -> **ABC** + D -> **ABCD**. An alternative linear route might start with a coupling of B + C as the first step. Another alternative might employ a convergent synthesis, i.e., A + B -> **AB**, then C + D -> **CD**, and then **AB** + **CD** -> **ABCD**. There are other possibilities. Even if the same four building blocks are used for each these are not the same syntheses, and in practice they would show different degrees of success. A trained synthetic chemist would be able to suggest the best order of couplings.

2. In some cases “protection group” strategies or the use of surrogate functional groups and reactions may be needed or be more efficient.

3. In some cases alternative strategies to the synthesis might still be preferred. For example, in the convergent 3 step synthesis of **ABCD** shown above. If the final coupling is somehow incompatible (or not optimal) it might be necessary to employ a synthesis such as C + E -> **CE** and then do the coupling **AB** + **CE** -> **ABCE** followed by **ABCE** -> **ABCD**. Although this adds an extra reaction, if it avoids some incompatibility (or inefficiency) between the component **D** and coupling used, then this would be an acceptable solution (indeed issues like this are common – more often than not!). Again a trained synthetic chemist would be able to easily suggest such strategies (or alternatively software packages that explicitly design syntheses).

4. Careful choice of reaction conditions, external reagents, catalysts, etc. are also essential aspects of any synthesis. These are not established with RGFN.

We hope that this answers the question.

Also, regarding another remark:

> Pushing back though slightly on the importance of the template-based approach for reaction prediction used here, many of these other works often seem fairly agnostic to the reaction oracles used.

Could we please ask for a clarification on whether the Reviewer refers to the discussion in the shared rebuttal, Originality / Similar Works paragraph, specifically the reliance on pre-trained neural networks to predict synthetic trees? As we are not entirely sure to which point does that comment refer to.
Thanks to the authors for their detailed rebuttal and responding to each of my points. Particularly appreciate the additional experiments provided to address some of my/other reviewers' concerns.

I have one minor follow up question regarding Q6, the answer to which I do not think I fully understand (partly also due to the response to Reviewer rpKJ's Question 2). I realize that RGFN might not necessarily generate the _best_ synthetic route to a proposed molecule, but I still do not understand why it does not generate _an_ explicit synthetic route. Are these not provided just by tracing back the sampling process for a product back to the initial building blocks?

Some general comments on the other points made:
* I appreciate that restricting to a smaller fragment library increases robustness (global response), although I am glad to see you recognise that the one used here is particularly small for many use cases and are taking steps to address this (and also support more complicated synthesis plans in future work).
* Glad you will mention the related work initially overlooked. Pushing back though slightly on the importance of the template-based approach for reaction prediction used here, many of these other works often seem fairly agnostic to the reaction oracles used.
* Thanks again for the new results (particularly with the quick turnaround), and also for discovering and fixing a mistake with the FGFN baseline. I still think it's a little weird that FGFN does not do particularly well on optimization (Q3), but thanks for the intuition regarding the fragment library/implementation provided in your rebuttal.

Will begin discussing the paper with the other reviewers and (unless anything new comes up from this) increase my score to reflect you addressing several of my points.
Thank you for the detailed and informative clarification. Most of my concerns are resolved. I believe that the paper has a basic value to be accepted. Therefore, I'd like to increase my score (5->6).