Generative Flows on Synthetic Pathway for Drug Design

## Abstract

Abstract Generative models in drug discovery have recently gained attention as efficient alternatives to brute-force virtual screening. However, most existing models do not account for synthesizability, limiting their practical use in real-world scenarios. In this paper, we propose RxnFlow , which sequentially assembles molecules using predefined molecular building blocks and chemical reaction templates to constrain the synthetic chemical pathway. We then train on this sequential generating process with the objective of generative flow networks (GFlowNets) to generate both highly rewarded and diverse molecules. To mitigate the large action space of synthetic pathways in GFlowNets, we implement a novel action space subsampling method. This enables RxnFlow to learn generative flows over extensive action spaces comprising combinations of 1.2 million building blocks and 71 reaction templates without significant computational overhead. Additionally, RxnFlow can employ modified or expanded action spaces for generation without retraining, allowing for the introduction of additional objectives or the incorporation of newly discovered building blocks. We experimentally demonstrate that RxnFlow outperforms existing reaction-based and fragment-based models in pocket-specific optimization across various target pockets. Furthermore, RxnFlow achieves state-of-the-art performance on CrossDocked2020 for pocket-conditional generation, with an average Vina score of –8.85 kcal/mol and 34.8% synthesizability.

## 1 Introduction

Structure-based drug discovery (SBDD) has emerged as a pivotal paradigm for early drug discovery, facilitated by the increasing accessibility of protein structure prediction tools (Jumper et al., ) and high-resolution crystallography (Liu et al., ).
However, traditional brute-force virtual screening is computationally expensive (Graff et al., ), prompting the development of deep generative models that can bypass this inefficiency.
In this context, various approaches such as deep reinforcement learning (Zhavoronkov et al., ), variational autoencoders (Zhung et al., ) generative adversarial network (Ragoza et al., ), and diffusion models (Guan et al., ; ) have been proposed to directly sample candidate molecules against a given protein structure.

While generative models have shown success in molecular discovery with desirable biological properties, most overlook synthesizability which is a crucial factor for wet-lab validation (Gao & Coley, ).
One line to improve synthesizability is multi-objective optimization using cheap functions to estimate the synthesizability (Ertl & Schuffenhauer, ), but this is too simplified to reflect complex synthetic principles (Cretu et al., ).
Other efforts aim to project molecules from generative models into a synthesizable space (Gao et al., ; Luo et al., ), but chemical modifications in this process can degrade the optimized properties.

To address this issue, recent works formulate the generation of synthetic pathways as a Markov decision process (MDP) for molecular design (Gottipati et al., ).
These approaches return a synthesizable molecule by assembling purchasable building blocks and reaction templates according to the generated synthetic pathway.
Notably, the emergence of virtual libraries—created by combinatorially enumerating building blocks and reaction templates, such as Enamine REAL (Grygorenko et al., )—allows the generated molecules to be readily synthesizable on demand with their synthetic pathways.
Recent studies (Cretu et al., ; Koziarski et al., ) advanced this approach by training the decision-making policy using the objective of generative flow networks (GFlowNets; Bengio et al., ).
This objective encourages the policy to sample in proportion to the reward function, enabling the retrieval of samples from a diverse range of modes.

Figure: Figure 1: Overview of RxnFlow. (a) Synthetic action space which is represented in a continuous action space. Each colored box corresponds to a reaction template and the molecules in the box are reactant blocks. (b) Policy estimation using the action space subsampling in a manner of importance sampling. (c) Molecular generation process and model training.
Refer to caption: extracted/5905292/Figures/action_space_subsampling2.png

Unlike atom-based or fragment-based models, the synthetic action spaces are massive and composed of millions of building blocks and tens of reaction templates.
While the large action spaces offer opportunities to discover novel hit candidates by expanding an explorable chemical space (Sadybekov et al., ), it incurs significant computational overhead.
Thus, prior works have restricted the action spaces to trade off the size of the search space for efficiency.
However, reducing the search space leads to a decrease in diversity and synthetic complexity.
While one could compensate for the reduced number of building block candidates by adding more reaction steps, this leads to the increase in synthetic complexity, negatively influencing synthesizability, yield, and cost (Coley et al., ; Kim et al., ).

In response to this challenge, we propose RxnFlow, a synthesis-oriented generative framework that allows training generative flows over a large action space to generate synthetic pathways for drug design. The distinctive features of this method are as follows.
First, we introduce an action space subsampling (Figure 1) to handle massive action spaces without significant memory overhead, enabling us to explore a broader chemical space with fewer reaction steps than existing models.
Then, we train the generative policy with a GFlowNet objective to sample both diverse and potent molecules from the expanded search space.
We demonstrate that RxnFlow effectively generates drug candidates, outperforming existing reaction-based, atom-based, and fragment-based baselines across various SBDD tasks, while ensuring the synthesizability of generated drug candidates.
We also achieve a new state-of-the-art Vina score, drug-likeness, and synthesizability on the CrossDocked2020 pocket-conditional generation benchmark (Luo et al., ).

Furthermore, we formulate an adaptable MDP (Figure 2) for consistent flow estimation on modified building block libraries, which can be highly practical in real-world applications.
By combining the proposed MDP with action embedding (Dulac-Arnold et al., ), which represents actions in a continuous space instead of a discrete space, RxnFlow can achieve further objectives or incorporate newly discovered building blocks without retraining.
We experimentally show that RxnFlow can achieve an additional solubility objective and behave appropriately for unseen building blocks.
This capability makes RxnFlow highly adaptable to real-world drug discovery pipeline, where new objectives frequently arise (Fink et al., ) and building block libraries are continuously expanding (Grygorenko et al., ).

## 2 Related Works

Structure-based drug discovery.
The first type of SBDD involves pocket-specific optimization methods to enhance docking scores against a single pocket, including evolutionary algorithms (Reidenbach, ), reinforcement learning (RL) (Zhavoronkov et al., ), and GFlowNets (Bengio et al., ).
However, these require individual optimizations for each pocket, limiting scalability.
The second type is based on pocket-conditional generation, which generates molecules against arbitrary given pockets without additional training.
This can be achieved by distribution-based generative models (Ragoza et al., ; Peng et al., ; Guan et al., ; Schneuing et al., ) trained on protein-ligand complex datasets to model ligand distributions for given pockets.
On the other hand, Shen et al. () formulated a pocket-conditioned policy for GFlowNets to generate samples from reward-biased distributions in a zero-shot manner.

Syntheis-oriented generative models.
To ensure the synthesizability of generated molecules, synthesis-oriented de novo design approaches incorporate combinatorial chemistry principles into generative models.
Bradshaw et al. () represented the synthetic pathways as directed acyclic graphs (DAG) for generative modeling.
Horwood & Noutahi () formulated the synthetic pathway generation as an MDP and optimize the molecules with RL.
Similarly, Gao et al. () employed a genetic algorithm to optimize synthesis trees to generate molecules with the desired properties.
Seo et al. () proposed
a conditional generative model to directly sample molecules with desired properties without optimization.
Recently, Cretu et al. (); Koziarski et al. () trained the models with the GFlowNets objective to generate diverse and potent molecules.

Action embedding for large action spaces.
To handle large action spaces in the synthesis-oriented generation, Gottipati et al. () employed action embedding (Dulac-Arnold et al., ) that represents building blocks in a continuous action space with their chemical information.
Later, Seo et al. (); Koziarski et al. () experimentally demonstrated that it can enhance the model training and generative performance.
The continuous action space provides the benefit of reducing the computational complexity for sampling from large space of actions and the memory requirement for parameterizing the categorical distribution over the large action space.

Generative Flow Networks.
GFlowNets are a learning framework for a stochastic generative policy that constructs an object through a series of decisions, where the probability of generating each object is proportional to a given reward associated with that object (Bengio et al., ).
Unlike other optimization methods that maximize rewards and often converge to a single solution, GFlowNets aim to sample a diverse set of high-rewarded modes, which is vital for novel drug design (Shen et al., ; Jain et al., ).
To this end, the generative policy is trained using objectives such as flow matching (Bengio et al., ), detailed balance (Bengio et al., ), and trajectory balance (Malkin et al., ).
Extending the GFlowNets to various applications is an active area of research,
e.g., GFlowNets have been applied to designing crystal structures (Nguyen et al., ), phylogenetic inference (Zhou et al., ), finetuning diffusion models (Venkatraman et al., ), and causal inference (Nguyen et al., ).

## 3 Method

### 3.1 GFlowNet Preliminaries

GFlowNets (Bengio et al., ) are the class of generative models that learn to sample objects $x\in\mathcal{X}$ proportional to a given reward function, i.e., $p(x)\propto R(x)$.
This is achieved by sequentially constructing a compositional object $x$ through a series of state transitions $s\rightarrow s^{\prime}$, forming a trajectory $\tau=(s_{0}\rightarrow\ldots\rightarrow s_{n}=x)\in\mathcal{T}$.
The set of all complete trajectories from the initial state $s_{0}$ can be represented as a directed acyclic graph $\mathcal{G}=(\mathcal{S},\mathcal{A})$ with a reachable state space $\mathcal{S}$ and an action space $\mathcal{A}$.
Each action $a$ induces a transition from the state $s$ to the state $s^{\prime}$, expressed as $s^{\prime}=T(s,a)$ and represented as $s\rightarrow s^{\prime}$.
Then, we define the trajectory flow $F(\tau)$, which flows along the trajectory $\tau=(s_{0}\rightarrow\ldots\rightarrow s_{n}=x)$, as the reward of the terminal state, $R(x)$.
The edge flow $F(s\rightarrow s^{\prime})$, or equivalently $F(s,a)$, is defined as the total flow along the edge $a:s\rightarrow s^{\prime}$:

$$ $F(s\rightarrow s^{\prime})=F(s,a)=\sum_{\tau\in\mathcal{T}\text{ s.t. }(s \rightarrow s^{\prime})\in\tau}F(\tau).$ (1) $$

The state flow $F(s)$ for the intermediate state is defined as the total flow through the state $s$:

$$ $F(s)=\sum_{\tau\in\mathcal{T}\text{ s.t. }s\in\tau}F(\tau)=\underbrace{\sum_{( s^{\prime\prime}\rightarrow s)\in\mathcal{A}}F(s^{\prime\prime}\rightarrow s)= \sum_{(s\rightarrow s^{\prime})\in\mathcal{A}}F(s\rightarrow s^{\prime}).}_{ \text{Intermediate flow matching condition}}$ (2) $$

In addition to the flow matching condition for intermediate states, there are two boundary conditions for the states.
First, the flow of a terminal state $x$ must equal the reward of the objective: $F(x)=R(x)$.
Second, the partition function, $Z$, is equivalent to the sum of all trajectory flows and the sum of all rewards: $Z=F(s_{0})=\sum_{\tau\in\mathcal{T}}F(\tau)=\sum_{x\in\mathcal{X}}R(x)$.
These three conditions—one for intermediate states and two for boundary states—are known as the flow matching conditions and ensure that GFlowNets generate objectives proportional to their rewards.

To convert the flow network into a usable policy, we define the forward policy as the forward transition probability $P_{F}(s^{\prime}|s)$ and the backward policy as the backward transition probability $P_{B}(s|s^{\prime})$:

$$ $P_{F}(s^{\prime}|s):=P(s\rightarrow s^{\prime}|s)=\frac{F(s\rightarrow s^{ \prime})}{F(s)},\quad P_{B}(s|s^{\prime}):=P(s\rightarrow s^{\prime}|s^{\prime })=\frac{F(s\rightarrow s^{\prime})}{F(s^{\prime})}.$ (3) $$

### 3.2 Action Space for Synthetic Pathway Generation

Following Cretu et al. (), we treated a chemical reaction as a forward transition and a synthetic pathway as a trajectory for molecular generation.
For the initial state $s_{0}$, the model always chooses AddFirstReactant to sample a building block $b$ from the entire building block set $\mathcal{B}$ as a starting molecule.
For the later states $s$, the model samples actions among ReactUni, ReactBi, or Stop.
When the action type is ReactUni, the model performs in silico uni-molecular reactions with an assigned reaction template $r\in\mathcal{R}_{1}$.
When the action type is ReactBi, the model performs bi-molecular reactions with a reaction template $r\in\mathcal{R}_{2}$ and a reactant block $b$ in the possible reactant set for the reaction template $r$: $\mathcal{B}_{r}\subseteq\mathcal{B}$.
If Stop is sampled, the trajectory is terminated.
To sum up, the allowable action space $\mathcal{A}(s)$ for the state $s$ is:

$$ $\mathcal{A}(s)=\displaystyle\begin{cases}\mathcal{B}&\text{if~{}}s=s_{0}\\ \{\texttt{Stop}\}\cup\mathcal{R}_{1}\cup\{(r,b)|r\in\mathcal{R}_{2},b\in \mathcal{B}_{r}\}&\text{otherwise}\end{cases}$ (4) $$

where unavailable reaction templates to the molecule of the state $s$ are masked.

### 3.3 Flow Network on Action Space Subsampling

We propose a novel memory-efficient technique called the action space subsampling, that estimates the state flow $F_{\theta}(s)$ from a subset of the outgoing edge flows $F_{\theta}(s\rightarrow s^{\prime})$ for forward policy estimation.
First, we implement an auxiliary policy, termed subsampling policy $\mathcal{P}(\mathcal{A})$, which samples a subset of the action space $\mathcal{A}^{\ast}\subseteq\mathcal{A}$.
This reduces both the memory footprint and the computational complexity from $\mathcal{O}(|\mathcal{B}||\mathcal{R}_{2}|)$ to $\mathcal{O}(|\mathcal{B}^{\ast}||\mathcal{R}_{2}|)$ with the controllable size $|\mathcal{B}^{\ast}|$.
We then estimate the forward policy by importance sampling.
In contrast to the parameterized forward policy, we formulate a fixed backward policy since it is hard to force invariance to molecule isomorphism (Malkin et al., ).
Theoretical backgrounds are provided in Sec. A.

Subsampling policy.
Subsampling policy $\mathcal{P}(\mathcal{A})$ performs uniform sampling for the initial state and importance sampling for the later states.
For the initial state, the allowable action space $\mathcal{A}(s_{0})=\mathcal{B}$ is homogeneous since all of them are AddFirstReactant actions.
For the later state, the action space is comprised of one Stop, tens of ReactUni actions, and millions of ReactBi actions.
To capture rare-type actions in the inhomogeneous space, we use all Stop and ReactUni actions.
The partial action space $\mathcal{A}^{\ast}(s)\sim\mathcal{P}(\mathcal{A}(s))$ comprises the uniform subset $\mathcal{B}^{\ast}\subseteq\mathcal{B}$ or $\mathcal{B}_{r}^{\ast}\subseteq\mathcal{B}_{r}$:

$$ $\mathcal{A}^{\ast}(s)=\displaystyle\begin{cases}\mathcal{B}^{\ast}&\text{if~{} }s=s_{0}\\ \{\texttt{Stop}\}\cup\mathcal{R}_{1}\cup\{(r,b)|r\in\mathcal{R}_{2},b\in \mathcal{B}^{\ast}_{r}\}&\text{otherwise~{}}\end{cases}$ (5) $$

Forward policy.
To estimate the forward policy $P_{F}(s^{\prime}|s)=F(s\rightarrow s^{\prime})/F(s)$ from the partial action space $\mathcal{A}^{\ast}$, we estimate the state flow $F(s)$ with a subset of outgoing edge flows $F(s\rightarrow s^{\prime})$.
Since we introduce importance sampling for action types, we weight the edge flow $F(s\rightarrow s^{\prime})$ of edge $a:s\rightarrow s^{\prime}$ according to the subsampling ratio:

$$ $w_{a}=w_{(s\rightarrow s^{\prime})}=\begin{cases}|\mathcal{B}|/|\mathcal{B}^{ \ast}|&\text{if}~{}a\in\mathcal{B}\\ |\mathcal{B}_{r}|/|\mathcal{B}^{\ast}_{r}|&\text{if}~{}a~{}\text{is}~{}(r,b)~{ }\text{where}~{}r\in\mathcal{R}_{2}\\ 1&\text{if}~{}a~{}\text{is}~{}\texttt{Stop}~{}\text{or}~{}a\in\mathcal{R}_{1} \\ \end{cases}$ (6) $$

By weighting edge flows, we can estimate the state flow $\hat{F}_{\theta}(s;\mathcal{A}^{\ast})$ as:

$$ $\hat{F}_{\theta}(s;\mathcal{A}^{\ast})=\sum_{(s\rightarrow s^{\prime})\in \mathcal{A}^{\ast}(s)}w_{(s\rightarrow s^{\prime})}F_{\theta}(s\rightarrow s^{ \prime}),$ (7) $$

which the estimated forward policy is $\hat{P}_{F}(s^{\prime}|s;\mathcal{A}^{\ast};\theta)=F_{\theta}(s\rightarrow s^
{\prime})/\hat{F}_{\theta}(s;\mathcal{A}^{\ast})$.

Action embedding.
In the standard implementation (Bengio et al., ) of the flow function $F_{\theta}$ and its neural network $\phi_{\theta}$, the edge flow of the edge $a:s\rightarrow s^{\prime}$ is computed with the corresponding action-specific parameter $\theta_{a}$: $F_{\theta}(s\rightarrow s^{\prime})=F_{\theta}(s,a)=\phi^{\text{flow}}_{\theta
_{a}}(\phi_{\theta}^{\text{state}}(s))$.

However, large action spaces require numerous parameters which increase model complexity.
To address this, we use an additional network $\phi^{\text{block}}_{\theta}$ for AddFirstReactant and ReactBi, which embeds the building block $b$ into a continuous action space with its structural information:

$$ $\displaystyle F_{\theta}(s_{0},b)=\phi^{\text{flow}}_{\theta}(\phi_{\theta}^{ \text{state}}(s),\phi_{\theta}^{\text{block}}(b)),\quad F_{\theta}(s,(r,b))= \phi^{\text{flow}}_{\theta}(\phi_{\theta}^{\text{state}}(s),\delta(r),\phi_{ \theta}^{\text{block}}(b))$ (8) $$

where $\delta(r)$ is the one-hot encoding for a bi-molecular reaction template $r$.

GFlowNet training.
In this work, we use the trajectory balance (TB; Malkin et al., ) as the training objective of GFlowNets from Eq. 9 and train models following Algorithm 1.
The action space subsampling is performed for each transition $s_{t}\rightarrow s_{t+1}$: $\mathcal{A}^{\ast}_{t}\sim\mathcal{P}(\mathcal{A})$.

$$ $\hat{\mathcal{L}}_{\text{TB}}(\tau)=\left(\log\frac{Z_{\theta}\prod^{n}_{t=1} \hat{P}_{F}(s_{t}|s_{t-1};\mathcal{A}^{\ast}_{t-1};\theta)}{R(x)\prod^{n}_{t=1 }P_{B}(s_{t-1}|s_{t})}\right)^{2}$ (9) $$

For online training, we use the sampling policy $\pi_{\theta}$ proportional to $\hat{P}_{F}(-|-;\mathcal{A}^{\ast};\theta)$, given by:

$$ $\pi_{\theta}(s^{\prime}|s;\mathcal{A}^{\ast})=\frac{w_{(s\rightarrow s^{\prime })}F_{\theta}(s\rightarrow s^{\prime})}{\sum_{(s\rightarrow s^{\prime\prime}) \in\mathcal{A}^{\ast}(s)}w_{(s\rightarrow s^{\prime\prime})}F_{\theta}(s \rightarrow s^{\prime\prime})}$ (10) $$

Figure: Algorithm 1 Training GFlowNets with action space subsampling

### 3.4 Joint Selection of Templates and Blocks

Figure: Figure 2: Comparison of using modified building block library for the generation: (a) a hierarchical MDP, and (b) a non-hierarchical MDP. More details are in Figure 7.
Refer to caption: extracted/5905292/Figures/mdp.png

For bi-molecular reactions, existing synthesis-oriented methods (Gao et al., ; Cretu et al., ; Koziarski et al., ) formulated a hierarchical MDP which selects a reaction template $r$ first and then the corresponding reactant block $b$ sequentially from Eq. 11.
However, as shown in Figure 2, the probability of selecting each reaction template $r$ is fixed after training in a hierarchical MDP, and this rigidity can lead to incorrect policy estimates in modified block libraries.
Therefore, we formulate a non-hierarchical MDP that jointly selects reaction templates and reactant blocks $(r,b)$ at once, as given by Eq. 12, resulting in more consistent estimates of forward policy $P_{F}$.

$$ $\displaystyle P_{F}(T(s,(r,b))|s;\theta)$ $\displaystyle=\frac{F_{\theta}(s,r)}{\sum_{r^{\prime}\in\mathcal{R}_{1}\cup \mathcal{R}_{2}\cup\{\texttt{Stop}\}}{F_{\theta}(s,r^{\prime})}}\times\frac{F_ {\theta}(s,(r,b))}{\sum_{b^{\prime}\in\mathcal{B}_{r}}F_{\theta}(s,(r,b^{ \prime}))}$ (11) $\displaystyle P_{F}(T(s,(r,b))|s;\theta)$ $\displaystyle=\frac{F_{\theta}(s,(r,b))}{\sum_{r^{\prime}\in\mathcal{R}_{1} \cup\{\texttt{Stop}\}}{F_{\theta}(s,r^{\prime})}+\sum_{r^{\prime}\in\mathcal{R }_{2}}\sum_{b^{\prime}\in\mathcal{B}_{r^{\prime}}}F_{\theta}(s,(r^{\prime},b^{ \prime}))}$ (12) $$

## 4 Experiments

Overview. We validate the effectiveness of RxnFlow in two common SBDD tasks: pocket-specific optimization (Sec. 4.1) and pocket-conditional generation (Sec. 4.2).
To the best of our knowledge, this is the first synthesis-oriented approach for pocket-conditional generation.
We also investigate the applicability of RxnFlow in real-world drug discovery pipelines where new further objectives may be introduced (Sec. 4.3) and the building block libraries are constantly expanded (Sec. 4.4).
Lastly, we conduct an ablation study in Sec. 4.5 and a theoretical analysis in Sec. D.3.

Setup.
We use the reaction template set constructed by Cretu et al. () including 13 uni- and 58 bi-molecular reaction templates.
For the building blocks, we use 1.2M blocks from the Enamine comprehensive catalog.
We use up to 3 reaction steps for generation following Enamine REAL Space (Grygorenko et al., ), while SynFlowNet and RGFN allow 4 steps.
For the subsampling policy, we set a sampling ratio of 1%.
The experimental details are provided in Sec. C.

Synthesizability estimation.
To assess the synthesizability of the generated compounds, we used the computationally intensive retrosynthetic analysis tool AiZynthFinder (Genheden et al., ) with the Enamine building block library.
We note that the molecule is identified as synthesizable only if it can be synthesized using the USPTO reactions (Lowe, ) and given building blocks.

### 4.1 Pocket-Specific Optimization with GPU-accelerated Docking

**Table 1: Hit ratio. Mean and standard deviation over 4 runs. The best results are in bold.**
|  |  | Hit Ratio (%, $\uparrow$) |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Category | Method | ADRB2 | ALDH1 | ESR_ago | ESR_antago | FEN1 |
| Fragment | FragGFN | 04.00 ($\pm$ 3.54) | 03.75 ($\pm$ 01.92) | 00.25 ($\pm$ 00.43) | 00.25 ($\pm$ 0.43) | 00.25 ($\pm$ 0.43) |
| FragGFN+SA | 05.75 ($\pm$ 1.48) | 04.00 ($\pm$ 01.58) | 00.25 ($\pm$ 00.43) | 00.00 ($\pm$ 0.00) | 00.00 ($\pm$ 0.00) |  |
| Reaction | SynNet | 45.83 ($\pm$ 7.22) | 25.00 ($\pm$ 25.00) | 00.00 ($\pm$ 00.00) | 00.00 ($\pm$ 0.00) | 50.00 ($\pm$ 0.00) |
| BBAR | 21.25 ($\pm$ 5.36) | 18.25 ($\pm$ 01.92) | 03.50 ($\pm$ 01.12) | 02.25 ($\pm$ 1.09) | 11.75 ($\pm$ 2.59) |  |
| SynFlowNet | 52.75 ($\pm$ 1.09) | 57.00 ($\pm$ 06.04) | 30.75 ($\pm$ 10.03) | 11.25 ($\pm$ 1.48) | 53.00 ($\pm$ 8.92) |  |
| RGFN | 46.75 ($\pm$ 6.87) | 39.75 ($\pm$ 08.17) | 04.50 ($\pm$ 01.66) | 01.25 ($\pm$ 0.43) | 19.75 ($\pm$ 4.32) |  |
| RxnFlow | 60.25 ($\pm$ 3.77) | 63.25 ($\pm$ 03.11) | 71.25 ($\pm$ 04.15) | 46.00 ($\pm$ 7.00) | 65.50 ($\pm$ 4.09) |  |

Setup.
Since GFlowNets sample a large number of molecules for online training, we employed a GPU-accelerated UniDock (Yu et al., ) with Vina scoring (Trott & Olson, ).
It is well known that docking can be hacked by increasing molecule size (Pan et al., ), so the appropriate constraints are required.
We select QED (Bickerton et al., ) as a comprehensive molecular property constraint, QED$>$0.5, and set the reward function as $R(x)=w_{1}\text{QED}(x)+w_{2}\widehat{\text{Vina}}(x)$
where $w_{1},w_{2}$ are used as the input of multi-objective GFlowNets (Jain et al., ) and are set to 0.5 for non-GFlowNet baselines.
$\widehat{\text{Vina}}$ is a normalized docking score (Eq. 30).

Each method generates up to 64,000 molecules for each of the 15 proteins in the LIT-PCBA dataset (Tran-Nguyen et al., ).
We then filter the molecules with the property constraint and select the top 100 diverse candidates based on the docking score, using a Tanimoto distance threshold of 0.5 to ensure structural diversity.
The selected molecules are evaluated with the following metrics:
Hit ratio (%) measures the fraction of hits, defined as the molecules that are identified as synthesizable by AiZynthFinder and having better docking scores than known active ligands (Lee et al., ).
Vina (kcal/mol) measures the average docking score.
Synthesizability (%) is the fraction of synthesizable molecules.
Synthetic complexity, which is highly correlated to yield and cost, is evaluated as the average number of synthesis steps (Coley et al., ).

Baselines.
We perform comparisons to various synthetic-oriented approaches: genetic algorithm (SynNet) (Gao et al., ), conditional generative model (BBAR(^1^11Since BBAR requires labeled training data with QED and docking score, we perform docking with random 62,720 ZINC molecules for training and evaluate 1,280 samples according to the reported splitting ratio.)) (Seo et al., ), and GFlowNets (SynFlowNet, RGFN) (Cretu et al., ; Koziarski et al., ).
For SynFlowNet and RGFN, we used 6,000 and 350 blocks, respectively, and set the maximum reaction step to 4 following the original papers.
Moreover, we consider fragment-based GFlowNets (FragGFN) to analyze the effects of synthetic constraints on the performance.
For FragGFN, we also consider additional synthesizability objectives with commonly-used synthetic accessibility score (SA; Ertl & Schuffenhauer, ) (FragGFN+SA).

Results.
The results for the first five targets are shown in Tables 4 and 4, and additional results for the 10 remaining targets are reported in Sec. D.1.
RxnFlow outperforms the baselines across all test proteins, demonstrating that the expanded sample space with the large action space enabled the model to generate more potent and diverse molecules.
Additionally, as shown in Tables 4 and 4, RxnFlow ensures the synthesizability of the generated molecules more effectively than the other synthesis-oriented methods and GFlowNets which employ the same reactions as ours.
These results support our primary assertion that existing synthesis-oriented approaches using more synthetic steps with smaller building block libraries can increase overall synthesis complexity and reduce synthesizability.
Furthermore, FragGFN+SA does not show meaningful improvement in synthesizability, implying that optimization of a cheap synthesizability estimation function is suboptimal.

Furthermore, it is noteworthy that RxnFlow outperformed FragGFN which does not consider synthesizability.
This improvement can be attributed to two key factors.
First, the Enamine building block library is specifically curated for drug discovery, limiting the search space $\mathcal{X}$ to a drug-like chemical space and simplifying the optimization complexity.
Second, RxnFlow needs shorter trajectories compared to FragGFNs since it assembles molecules with large building blocks instead of atoms or small fragments.
This is beneficial for trajectory balance objectives where stochastic gradient variance tends to increase over longer trajectories (Malkin et al., ).

### 4.2 Zero-Shot Sampling via Pocket Conditioning

**Table 5: CrossDocked2020 benchmark. We report the average and median values over the average properties for each test pocket. The best results are in bold and the second ones are in underlined. We denote the Validity as Valid, Synthesizability as Synthesiz, and Diversity as Div. Reference means the known binding active molecules in the test set.**
|  |  | Valid ($\uparrow$) | Vina ($\downarrow$) | QED ($\uparrow$) | Synthesiz ($\uparrow$) | Div ($\uparrow$) | Time ($\downarrow$) |  |  |  |
| --- | --- | --- | --- | --- | --- | --- | --- | --- | --- | --- |
| Category | Model | Avg. | Avg. | Med. | Avg. | Med. | Avg. | Med. | Avg. | Avg. |
| Reference | - | - | -7.71 | -7.80 | 0.48 | 0.47 | 36.1% | - | - | - |
| Atom | Pocket2Mol | 098.3% | -7.60 | -7.16 | 0.57 | 0.58 | 29.1% | 22.0% | 0.83 | 2504^a |
| TargetDiff | 091.5% | -7.37 | -7.56 | 0.49 | 0.49 | 09.9% | 03.2% | 0.87 | 3428^a |  |
| DecompDiff | 066.0% | -8.35 | -8.25 | 0.37 | 0.35 | 00.9% | 00.0% | 0.84 | 6189^a |  |
| Fragment | TacoGFN | 100.0% | -8.24 | -8.44 | 0.67 | 0.67 | 01.3% | 01.0% | 0.67 | 0004^b |
| Reaction | RxnFlow | 100.0% | -8.85 | -9.03 | 0.67 | 0.67 | 34.8% | 34.5% | 0.81 | 0004^b |

Figure: Figure 3: Visualization of generated molecules in a zero-shot manner. (a-b) Docking results of generated molecules and known reference ligands of TBK1 (PDB Id: 1FV, SU6). (c) Generative trajectory, which is the generated synthetic pathway of the left molecule in (a).
Refer to caption: extracted/5905292/Figures/result/generated_molecules_tiny.png

Setup.
We extend our works to a pocket-conditional generation problem to design binders for arbitrary pockets without additional training oracles (Ragoza et al., ; Liu et al., ; Peng et al., ; Guan et al., ; ; Schneuing et al., ).
To address this challenge, we follow TacoGFN (Shen et al., ), which is a fragment-based GFlowNet for pocket-conditional generation.
Since it requires more training oracles to learn pocket-conditional policies than target-specific generation, TacoGFN used pre-trained proxies that predict docking scores against arbitrary pockets using pharmacophore representation (Seo & Kim, ).
Since RxnFlow explicitly considers synthesizability, we exclude the SA score from the TacoGFN’s reward function as described in Sec. C.2.

We generate 100 molecules for each of the 100 test pockets in the CrossDocked2020 benchmark (Francoeur et al., ) and evaluate them with the following metrics:
Vina (kcal/mol) measures the average docking score from QuickVina (Alhossary et al., ).
QED measures the average drug-likeness of molecules.
Synthesizability (%) is the fraction of synthesizable molecules.
Diversity measures an average pairwise Tanimoto distance of ECFP4 fingerprints (Morgan, ).
Moreover, we report the Validity (%) which is the percentage of unique RDKit-readable molecules, and Time (sec.) which is the average sampling time to generate 100 molecules.

Baselines.
We compare RxnFlow with state-of-the-art distribution learning-based models trained on a synthesizable drug set: an autoregressive model (Pocket2Mol (Peng et al., )) and diffusion models (TargetDiff (Guan et al., ), DecompDiff (Guan et al., )).
We also perform a comparison with an optimization-based TacoGFN (Shen et al., ).
For a fair comparison with the distribution learning-based approaches, we used the docking proxy trained on CrossDocked2020.

Result.
As shown in Table 5, RxnFlow achieves significant improvements in drug-related properties.
In particular, RxnFlow attains the best average docking score of -8.85 kcal/mol, significantly outperforming TacoGFN at -8.24 kcal/mol and the state-of-the-art diffusion model (DecompDiff) at -8.35 kcal/mol while achieving high drug-like properties.
Moreover, RxnFlow ensures the synthesizability comparable to known active ligands, outperforming the fragment-based TacoGFN trained on the SA score objective and the distributional learning-based models trained on synthesizable drug molecules.
Figure 3 illustrates generated molecules against TANK-binding kinase 1 (TBK1) which is not included in the training set.

An important finding is that RxnFlow maintains high structural diversity (0.81) despite the typical trade-off between optimization power and diversity (Gao et al., ).
This is a significant improvement over the fragment-based TacoGFN, which scored 0.67 and is comparable to the distributional learning-based models that range from 0.83 to 0.87.
We attribute this enhancement to our action space, which contains chemically diverse building blocks, in contrast to the small and limited fragment sets used in fragment-based GFlowNets.
This suggests that our model can effectively balance the potency and diversity of generated molecules.

Figure: Figure 4: Property distribution of sampled molecules with “all” building blocks and “low”-TPSA building blocks. Vina score was calculated against the KRAS-G12C target.
Refer to caption: extracted/5905292/Figures/result/ana2_further_objective.png

### 4.3 Introducing Additional Objective without Retraining

In drug discovery, new objectives often arise during the research process, such as enhancing solubility, reducing toxicity, or improving selectivity (Fink et al., ; Joshi et al., ).
These additional objectives typically not only require retraining models but also increase the optimization complexity.
In this context, RxnFlow can achieve some additional objectives by simply introducing constraints to MDP without retraining thanks to the non-hierarchical action space (Sec. 3.4).

As shown in Figure 5, we explore the scenario of adding a solubility objective to a pre-trained GFlowNet in Sec. 4.2.
Specifically, we target the generation of hydrophobic molecules by restricting the building blocks with topological polar surface area (TPSA) in the bottom 15% (“low”) and sampled 500 molecules for the KRAS-G12C mutant (PDB Id: 6oim).
While there are slight differences in the QED distributions due to the correlation between TPSA and QED, the generated molecules are more hydrophobic and retain similar overall reward distributions.

### 4.4 Scaling Action Space without Retraining

The building block libraries for drug discovery continue to grow, from 60,000 in 2020 to over 1.2 million blocks today, to enhance chemical diversity and novelty (Grygorenko et al., ).
However, existing generative models require retraining to accommodate newly discovered building blocks, limiting their scalability and adaptability.
On the other hand, RxnFlow can integrate new building blocks without retraining by understanding the chemical context of actions through action embedding.
We first divide the 1M-sized block library (“all”) into two sets: 500,000 blocks for training (“seen”) and the remaining blocks (“unseen”).
After training with the QED objective on various reward exponent settings ($R^{\beta}$), we generate 5,000 molecules from each set (“seen”, “unseen”, and “all”).
Figure 5 shows that the reward distributions of samples are nearly identical, demonstrating that RxnFlow performs robustly with unseen building blocks.
This result highlights the generalization ability and scalability of RxnFlow, a significant advantage for real-world applications.

Figure: Figure 6: Optimization power, diversity, and generation time according to building block library size. (a-c) Average and standard deviation of properties of the top-1000 high-affinity molecules over 4 runs on pocket-specific generation. (a) Average docking score. (b) The uniqueness of Bemis-Murcko scaffolds. (c) Average Tanimoto distance. (d) Average runtime to generate 100 molecules over the CrossDocked2020 test pockets in a zero-shot manner.
Refer to caption: extracted/5905292/Figures/result/ana3_ablation.png

### 4.5 Ablation Study: the effect of building block library size

Expanding the size of building block libraries provides an opportunity to discover more diverse and potent drug candidates (Grygorenko et al., ).
In Sec. 4.1, RxnFlow outperforms SynFlowNet and RGFN which use smaller block libraries, but differences in model architectures may have contributed to these results.
To isolate the effect of the building block library size, we conduct an ablation study using partial libraries with a pocket-specific optimization task against the kappa-opioid receptor (PDB Id: 6b73), as illustrated in Figure 6(a-c).
The results indicate that increasing the library size enhances both optimization power, in terms of docking scores, and diversity, in terms of a higher number of unique Bemis-Murcko scaffolds (Bemis & Murcko, ) and an increased Tanimoto diversity of the generated molecules.
Additionally, as shown in Figure 6(d), the generation time only doubles on the 100-fold larger action space, highlighting the efficiency of RxnFlow.
These results demonstrate the forte of RxnFlow in navigating a broader chemical space to discover novel drug candidates by overcoming the computational limitations for expanding the action space.

## 5 Conclusion

In this work, we introduce RxnFlow, a synthesis-oriented generative framework designed to explore broader chemical spaces, thereby enhancing both diversity and potency for drug discovery.
Our framework efficiently handles massive action spaces to expand the search space without significant computational or memory overhead by employing a novel action space subsampling technique.
RxnFlow can effectively identify diverse drug candidates with desired properties and synthetic feasibility by learning the objective of generative flow networks on synthetic pathways.
Additionally, by formulating a non-hierarchical MDP, RxnFlow can model generative flows on modified action spaces, allowing it to achieve additional objectives and incorporate newly discovered building blocks without retraining.
These results highlight the potential of RxnFlow as a practical and versatile solution for real-world drug discovery.

## Appendix A Theoretical Analysis

In this section, we provide the theoretical background of action subsampling.
We define $\mathcal{U}$ as the uniform subsampling policy, i.e., $\mathcal{B}^{\ast}\sim\text{Uniform}(\{\mathcal{B}^{\ast}\subseteq\mathcal{B}~
{}|~{}|\mathcal{B}^{\ast}|=k\})$.

#### Bias of log forward policy.

In this section, we prove that the log forward policy estimation is unbiased (Eq. 7) with the following weights:

$$ $w_{a}=\begin{cases}|\mathcal{B}|/|\mathcal{B}^{\ast}|&\text{if}~{}a\in\mathcal {B}\\ |\mathcal{B}_{r}|/|\mathcal{B}^{\ast}_{r}|&\text{if}~{}a~{}\text{is}~{}(r,b)~{ }\text{and}~{}r\in\mathcal{R}_{2}\\ 1&\text{if}~{}a~{}\text{is}~{}\texttt{Stop}~{}\text{or}~{}a\in\mathcal{R}_{1} \\ \end{cases}$ $$

For readability, we express the edge flow $F$ and forward policy $P_{F}$ where the action $a$ is $s\rightarrow s^{\prime}$ as follows:

$$ $F_{\theta}(s\rightarrow s^{\prime};\theta)=F_{\theta}(s,a;\theta),\quad P_{F}( s^{\prime}|s;\theta)=P_{F}(a|s;\theta)$ $$

Then, the forward policy $P_{F}$ and the estimated policy $\hat{P}_{F}$ (Eq. 7) can be rewritten as follows:

$$ $\displaystyle P_{F}(s^{\prime}|s;\theta)$ $\displaystyle=P_{F}(a|s;\theta)=\frac{F_{\theta}(s,a)}{F_{\theta}(s)}=\frac{F_ {\theta}(s,a)}{\sum_{a^{\prime}\in\mathcal{A}(s)}F_{\theta}(s,a^{\prime})}$ $\displaystyle\hat{P}_{F}(s^{\prime}|s;\mathcal{A}^{\ast};\theta)$ $\displaystyle=\hat{P}_{F}(a|s;\mathcal{A}^{\ast};\theta)=\frac{F_{\theta}(s,a) }{\hat{F}_{\theta}(s;\mathcal{A}^{\ast})}=\frac{F_{\theta}(s,a)}{\sum_{a^{ \prime}\in\mathcal{A}^{\ast}(s)}w_{a^{\prime}}F_{\theta}(s,a^{\prime})}$ $$

The expectation of the estimated initial state flow $\hat{F}_{\theta}(s_{0};\mathcal{A}^{\ast})=\sum_{a\in\mathcal{A}^{\ast}}w_{a}F
_{\theta}(s_{0},a)$ is given by

$$ $\displaystyle\mathbb{E}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}[\hat{ F}_{\theta}(s_{0};\mathcal{A}^{\ast})]=\mathbb{E}_{\mathcal{B}^{\ast}\sim \mathcal{U}(\mathcal{B})}\left[\frac{|\mathcal{B}|}{|\mathcal{B}^{\ast}|}\sum_ {b\in\mathcal{B}^{\ast}}F_{\theta}(s_{0},b)\right]=\sum_{b\in\mathcal{B}}F_{ \theta}(s_{0},b)=F_{\theta}(s_{0}),$ (13) $$

For a later state $s\neq s_{0}$, the expectation of the state flow $\hat{F}_{\theta}(s;\mathcal{A}^{\ast})$ is given by

$$ $\displaystyle\mathbb{E}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}[\hat{ F}_{\theta}(s;\mathcal{A}^{\ast})]$ $\displaystyle=\sum_{a\in\mathcal{R}_{1}\cup\{\texttt{Stop}\}}F_{\theta}(s,a)+ \sum_{r\in\mathcal{R}_{2}}\mathbb{E}_{\mathcal{B}^{\ast}_{r}\sim\mathcal{U}( \mathcal{B}_{r})}\left[\frac{|\mathcal{B}_{r}|}{|\mathcal{B}^{\ast}_{r}|}\sum_ {b\in\mathcal{B}^{\ast}_{r}}F_{\theta}(s,(r,b))\right]$ $\displaystyle=\sum_{a\in\mathcal{R}_{1}\cup\{\texttt{Stop}\}}F_{\theta}(s,a)+ \sum_{r\in\mathcal{R}_{2}}\sum_{b\in\mathcal{B}_{r}}F_{\theta}(s,(r,b))$ $\displaystyle=P_{F}(s^{\prime}|s;\theta)$ (14) $$

#### Variance of log forward policy.

We define the standard deviation of the probability $P_{F}(-|s;\theta)$ as $\sigma_{\theta}(-|s)$.
For the initial state $s_{0}$, the variance of $\log\hat{P}_{F}(s^{\prime}|s_{0};\mathcal{A}^{\ast};\theta)$ is given by

$$ Var ${}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}\left[\log\hat{P}_{F}(s^{ \prime}|s_{0};\mathcal{A}^{\ast};\theta)\right]$ $\displaystyle=\text{Var}_{\mathcal{B}^{\ast}\sim\mathcal{U}(\mathcal{B})}\left [\cancel{\log F_{\theta}(s_{0}\rightarrow s^{\prime})}-\log\left(\cancel{\frac {|\mathcal{B}|}{|\mathcal{B}^{\ast}|}}\sum_{b\in\mathcal{B}^{\ast}}F_{\theta}( s_{0},b)\right)\right]$ $\displaystyle=\text{Var}_{\mathcal{B}^{\ast}\sim\mathcal{U}(\mathcal{B})}\left [\log\sum_{b\in\mathcal{B}^{\ast}}P_{F}(b|s_{0};\theta)\right]$ $\displaystyle\approx\frac{|\mathcal{B}|^{2}(|\mathcal{B}|-|\mathcal{B}^{\ast}| )}{|\mathcal{B}^{\ast}|(|\mathcal{B}|-1)}\sigma_{\theta}(-|s_{0})^{2}\quad \text{(by \lx@cref{creftypecap~refnum}{eq-appendix: variance-log})}$ (15) $$

For the later state $s\neq s_{0}$, the variance of $\log\hat{P}_{F}(s^{\prime}|s;\mathcal{A}^{\ast};\theta)$ is:

$$ $\displaystyle\text{Var}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}\left[ \log\hat{P}_{F}(s^{\prime}|s;\mathcal{A}^{\ast};\theta)\right]$ $\displaystyle=\text{Var}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}\left [\log\left(\sum_{a\in\mathcal{R}_{1}\cup\{\texttt{Stop}\}}P_{F}(a|s;\theta)+ \sum_{r\in\mathcal{R}_{2}}\frac{|\mathcal{B}_{r}|}{|\mathcal{B}^{\ast}_{r}|} \sum_{b\in\mathcal{B}^{\ast}_{r}}P_{F}((r,b)|s;\theta)\right)\right]$ $\displaystyle\approx\frac{\text{Var}_{\mathcal{A}^{\ast}\sim\mathcal{P}( \mathcal{A})}\left[\cancel{\sum_{a\in\mathcal{R}_{1}\cup\{\texttt{Stop}\}}P_{F }(a|s;\theta)}+\sum_{r\in\mathcal{R}_{2}}\frac{|\mathcal{B}_{r}|}{|\mathcal{B} ^{\ast}_{r}|}\sum_{b\in\mathcal{B}^{\ast}_{r}}P_{F}((r,b)|s;\theta)\right]}{ \mathbb{E}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}\left[\sum_{a\in \mathcal{R}_{1}\cup\{\texttt{Stop}\}}P_{F}(a|s;\theta)+\sum_{r\in\mathcal{R}_{ 2}}\frac{|\mathcal{B}_{r}|}{|\mathcal{B}^{\ast}_{r}|}\sum_{b\in\mathcal{B}^{ \ast}_{r}}P_{F}((r,b)|s;\theta)\right]^{2}}$ $\displaystyle=\frac{1}{\cancel{\left(\sum_{a\in\mathcal{A}(s)}P_{F}(a|s;\theta )\right)^{2}}}\text{Var}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A}(s))} \left[\sum_{r\in\mathcal{R}_{2}}\frac{|\mathcal{B}_{r}|}{|\mathcal{B}^{\ast}_{ r}|}\sum_{b\in\mathcal{B}^{\ast}_{r}}P_{F}((r,b)|s;\theta)\right]$ $\displaystyle=\sum_{r\in\mathcal{R}_{2}}\frac{|\mathcal{B}_{r}|^{2}(|\mathcal{ B}_{r}|-|\mathcal{B}^{\ast}_{r}|)}{|\mathcal{B}^{\ast}_{r}|(|\mathcal{B}_{r}|- 1)}\sigma_{\theta}((r,-)|s)^{2}\quad\text{(by \lx@cref{creftypecap~refnum}{eq- appendix: variance-mc})}$ (16) $$

Finally, we get the variance of $\log\hat{P}_{F}(\tau;\theta)$ where $\tau=(s_{0}\rightarrow...\rightarrow s_{n}=x)$:

$$ $\displaystyle\text{Var}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}$ $\displaystyle\left[\log\hat{P}_{F}(\tau;\mathcal{A}^{\ast};\theta)\right]$ $\displaystyle\approx$ $\displaystyle\frac{|\mathcal{B}|^{2}(|\mathcal{B}|-|\mathcal{B}^{\ast}|)}{k| \mathcal{B}^{\ast}|(|\mathcal{B}|-1)}\sigma_{\theta}(-|s_{0})^{2}+\sum_{t=1}^{ n-1}\sum_{r\in\mathcal{R}_{2}}\frac{|\mathcal{B}_{r}|^{2}(|\mathcal{B}_{r}|-| \mathcal{B}^{\ast}_{r}|)}{|\mathcal{B}^{\ast}_{r}|(|\mathcal{B}_{r}|-1)}\sigma _{\theta}((r,-)|s_{t})^{2}$ (17) $$

#### Expectation of trajectory balance loss.

Due to the variance of forward policy, there is a bias of the trajectory balance loss equal to the variance of $\log\hat{P}_{F}(\tau)$:

$$ $\displaystyle\mathbb{E}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}[ \mathcal{\hat{L}}_{\text{TB}}(\tau)]$ $\displaystyle=\mathbb{E}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}\left [\log\frac{Z_{\theta}\hat{P}_{F}(\tau;\theta)}{R(x)P_{B}(\tau|x)}\right]^{2}+ \text{Var}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}\left[\log\frac{Z_{ \theta}\hat{P}_{F}(\tau;\theta)}{R(x)P_{B}(\tau|x)}\right]$ $\displaystyle=\mathcal{L}_{\text{TB}}(\tau)+\text{Var}_{\mathcal{A}^{\ast}\sim \mathcal{P}(\mathcal{A})}\left[\log\hat{P}_{F}(\tau;\theta)\right]$ (18) $\displaystyle\approx\mathcal{L}_{\text{TB}}(\tau)+\frac{|\mathcal{B}|^{2}(| \mathcal{B}|-|\mathcal{B}^{\ast}|)}{|\mathcal{B}^{\ast}|(|\mathcal{B}|-1)} \sigma_{\theta}(-|s_{0})^{2}+\sum_{t=1}^{n-1}\sum_{r\in\mathcal{R}_{2}}\frac{| \mathcal{B}_{r}|^{2}(|\mathcal{B}_{r}|-|\mathcal{B}^{\ast}_{r}|)}{|\mathcal{B} ^{\ast}_{r}|(|\mathcal{B}_{r}|-1)}\sigma_{\theta}((r,-)|s_{t})^{2}$ (19) $$

Therefore, the gradient of trajectory balance loss is:

$$ $\displaystyle\mathbb{E}_{\mathcal{A}^{\ast}\sim\mathcal{P}(\mathcal{A})}\left[ \nabla_{\theta}\hat{\mathcal{L}}_{\text{TB}}(\tau)\right]$ $\displaystyle\approx\nabla_{\theta}\mathcal{L}_{\text{TB}}(\tau)+\frac{| \mathcal{B}|^{2}(|\mathcal{B}|-|\mathcal{B}^{\ast}|)}{|\mathcal{B}^{\ast}|(| \mathcal{B}|-1)}\nabla_{\theta}\sigma_{\theta}(-|s_{0})^{2}+\sum_{t=1}^{n-1} \sum_{r\in\mathcal{R}_{2}}\frac{|\mathcal{B}_{r}|^{2}(|\mathcal{B}_{r}|-| \mathcal{B}^{\ast}_{r}|)}{|\mathcal{B}^{\ast}_{r}|(|\mathcal{B}_{r}|-1)}\nabla _{\theta}\sigma_{\theta}((r,-)|s_{t})^{2}$ (20) $$

Given that $\sigma_{\theta}(\cdot)$, which is the standard deviation of probability to select actions, is inversely proportional to $|\mathcal{B}|$ or $|\mathcal{B}_{r}|$, the bias is highly dependent on the subsampling size (e.g. $|\mathcal{B}^{\ast}|$) rather than the subsampling ratio (e.g. $|\mathcal{B}^{\ast}|/|\mathcal{B}|$).
Moreover, we can decrease the bias by MC sampling for state flow estimation, and the bias is reversely proportional to the number of samples $k$.
However, in the same computational cost (proportional to $\mathcal{O}(k\times\sum_{r}|\mathcal{B}_{r}^{\ast}|)$), using the larger partial action spaces without MC sampling is more precise than using smaller partial action spaces with multiple MC samples.
In Sec. D.3, we expermentally show that the bias is relatively small to trajectory balance loss with the toy experiment.

### A.1 Theoretical backgrounds

#### Variance of uniformly sampled subset

For the set $\mathcal{X}$ with a size of $n$, we define its uniformly sampled subset with size $m$ as $\mathcal{X}^{\prime}\sim\mathcal{P}(\mathcal{X})$.
For the function $F(\mathcal{X})=\sum_{x\in\mathcal{X}}f(x)$, we define the unbiased estimation:

$$ $\hat{F}(\mathcal{X})=\frac{n}{m}F(\mathcal{X}^{\prime})=\frac{n}{m}\sum_{x\in \mathcal{X}^{\prime}}f(x)$ $$

When the mean and standard deviation of $f(x)$ is $\mu_{f}$ and $\sigma_{f}$, variance of $\hat{F}(\mathcal{X})$ is:

$$ $\displaystyle\text{Var}_{\mathcal{X}^{\prime}\sim\mathcal{P}(\mathcal{X})} \left[\hat{F}(\mathcal{X})\right]$ $\displaystyle=\mathbb{E}_{\mathcal{X}^{\prime}\sim\mathcal{P}(\mathcal{X})} \left[\hat{F}(\mathcal{X})^{2}\right]-\mathbb{E}_{\mathcal{X}^{\prime}\sim \mathcal{P}(\mathcal{X})}\left[\hat{F}(\mathcal{X})\right]^{2}$ $\displaystyle=\frac{1}{{}^{n}C_{m}}\sum_{\mathcal{X}^{\prime}\in\mathcal{P}( \mathcal{X})}\left(\frac{n}{m}\sum_{x\in\mathcal{X}^{\prime}}f(x)\right)^{2}-F (\mathcal{X})^{2}$ $\displaystyle=\frac{1}{{}^{n}C_{m}}\frac{n^{2}}{m^{2}}\left({}^{n-1}C_{m-1} \sum_{i=1}^{n}f(b_{i})^{2}+^{n-2}C_{m-2}\sum_{i=1}^{n}\sum_{j>i}^{n}2f(b_{i})f (b_{j})\right)-\left(\sum_{x\in\mathcal{X}}f(x)\right)^{2}$ $\displaystyle=\frac{n}{m}\sum_{i=1}^{n}f(x_{i})^{2}+\frac{2n(m-1)}{m(n-1)}\sum _{i=1}^{n}\sum_{j>i}^{n}f(x_{i})f(x_{j})-\left(\sum_{i=1}^{n}f(x_{i})^{2}+2 \sum_{i=1}^{n}\sum_{j>i}^{n}f(x_{i})f(x_{j})\right)$ $\displaystyle=\frac{n-m}{m(n-1)}\left((n-1)\sum_{i=1}^{n}f(x_{i})^{2}-2\sum_{i =1}^{n}\sum_{j>i}^{n}f(x_{i})f(x_{j})\right)$ $\displaystyle=\frac{n-m}{m(n-1)}\sum_{i=1}^{n}\sum_{j>i}^{n}(f(x_{i})-f(x_{j}) )^{2}$ $\displaystyle=\frac{n^{2}(n-m)}{m(n-1)}\sigma_{f}^{2}$ (21) $$

For MC sampling with $k$ samples, the variance is:

$$ $\displaystyle\text{Var}_{\mathcal{X}^{\prime}\sim\mathcal{P}(\mathcal{X})} \left[\hat{F}(\mathcal{X})\right]=\frac{n^{2}(n-m)}{km(n-1)}\sigma_{f}^{2}$ (22) $$

The variance of $\log\hat{F}(\mathcal{X})$ is:

$$ $\text{Var}_{\mathcal{X}^{\prime}\sim\mathcal{P}(\mathcal{X})}\left[\log\hat{F} (\mathcal{X})\right]\approx\frac{\text{Var}_{\mathcal{X}^{\prime}\sim\mathcal{ P}(\mathcal{X})}\left[\hat{F}(\mathcal{X})\right]}{\mathbb{E}_{\mathcal{X}^{ \prime}\sim\mathcal{P}(\mathcal{X})}\left[\hat{F}(\mathcal{X})\right]^{2}}= \frac{(n-m)}{km(n-1)}(\sigma_{f}/\mu_{f})^{2}$ (23) $$

When the $F(\mathcal{X})=1$, the variance is:

$$ $\text{Var}_{\mathcal{X}^{\prime}\sim\mathcal{P}(\mathcal{X})}\left[\log\hat{F} (\mathcal{X})\right]\approx\frac{m^{2}(n-m)}{km(n-1)}\sigma_{f}^{2}$ (24) $$

## Appendix B RxnFlow Architecture

### B.1 Forward policy

For the model architecture, we used a graph transformer (Yun et al., ) as the backbone $f_{\theta}$ following Bengio et al. () and a multi-layer perceptron (MLP) for action embedding $g_{\theta}$.
The graph embedding dimension is $d_{1}$, and the building block embedding dimension is $d_{2}$.
The molecular graph for a state is $s$, and the GFlowNet condition vector is $c$ which includes a reward exponent and multi-objective optimization weights (Jain et al., ).
$a\|b$ means the concatenation of two feature vectors $a$ and $b$.

#### Initial block selection.

For the first action, the model always selects AddFirstReactant action which selects $b\in\mathcal{B}$ for the starting molecule with $\text{MLP}_{\texttt{AddFirstReactant}}:\mathbb{R}^{d_{1}+d_{2}}\rightarrow
\mathbb{R}$.

$$ $F_{\theta}(s_{0},b,c)=\text{MLP}_{\texttt{AddFirstReactant}}(f_{\theta}(s_{0}, c)\|g_{\theta}(b))$ (25) $$

#### Reaction selection.

For the later states $s\neq s_{0}$, the model calculates the logits for the Stop action, ReactUni actions $r_{1}\in\mathcal{R}_{1}$, and ReactBi actions ($r_{2},b)\in\mathcal{R}_{2}\times\mathcal{B}$.

The logit for Stop is calculated by $\text{MLP}_{\texttt{Stop}}:\mathbb{R}^{d_{1}}\rightarrow\mathbb{R}$:

$$ $F_{\theta}(s,\texttt{Stop},c)=\text{MLP}_{\texttt{Stop}}(f_{\theta}(s,c)).$ (26) $$

The logit for a ReactUni action $r_{1}\in\mathcal{R}_{1}$ is calculated by $\text{MLP}_{\texttt{ReactUni}}^{r_{1}}:\mathbb{R}^{d_{1}}\rightarrow\mathbb{R}$:

$$ $F_{\theta}(s,r_{1},c)=\text{MLP}_{\texttt{ReactUni}}^{r_{1}}(f_{\theta}(s,c)).$ (27) $$

Finally, the logit for a ReactBi action $(r_{2},b)$ is calculated by the one-hot embedding of reaction template $\delta(r_{2}):\{0,1\}^{|\mathcal{R}_{2}|}$ and $\text{MLP}_{\texttt{ReactBi}}:\mathbb{R}^{d_{1}+|\mathcal{R}_{2}|+d_{2}}
\rightarrow\mathbb{R}$:

$$ $F_{\theta}(s,(r_{2},b),c)=\text{MLP}_{\texttt{ReactBi}}(f_{\theta}(s,c)\| \delta(r_{2})\|g_{\theta}(b)).$ (28) $$

#### Pocket conditioning.

For a pocket-conditional generation, the model uses a K-NN pocket residual graph $\mathcal{G}^{P}$ and encodes GVP-GNN (Jing et al., ) according to Shen et al. ().
The pocket conditions are included in the GFlowNet condition vector $c$.

### B.2 Backward policy.

Malkin et al. () introduced a uniform backward policy for small-scale drug discovery.
However, in a directed acyclic graph (DAG) on synthetic pathways, where each state has numerous outgoing edges but few incoming ones, there is a significant imbalance in the number of trajectories along each incoming edge, depending on the distance from the initial state.
In a uniform backward policy, where the flow of all incoming edges is equal, this imbalance diminishes the flow of trajectories that reach a state via shorter routes, i.e., those with fewer reaction steps.
To facilitate shorter synthetic pathways as trajectories, we set the backward transition probability proportional to the expected number of trajectories along each incoming edge of the state, $a:s^{\prime\prime}\rightarrow s$:

$$ $P_{B}(s^{\prime\prime}|s)=P((s^{\prime\prime}\rightarrow s)\in\tau|s\in\tau)= \frac{\sum_{\tau\in\mathcal{T}:\tau=(s_{0}\rightarrow...\rightarrow s^{\prime \prime}\rightarrow s)}|\mathcal{A}(s)|^{N-|\tau|}}{\sum_{\tau\in\mathcal{T}: \tau=(s_{0}\rightarrow...\rightarrow s)}|\mathcal{A}(s)|^{N-|\tau|}}$ (29) $$

where $N$ is the maximum length of trajectories.

### B.3 Building block representation for action embedding.

To represent building blocks, we used both physicochemical properties, chemical structural properties, and topological properties.
For physicochemical properties, we used 8 molecular features: molecular weight, the number of atoms, the number of H-bond acceptors/donors (HBA, HBD), the number of aromatic/non-aromatic rings, LogP, and TPSA.
For chemical properties, we used the MACCS fingerprint (Durant et al., ), which represents the composition of the chemical functional groups in the molecule.
For topological information, we used the Morgan ECFP4 fingerprint (Morgan, ) with a dimension of 1024, which is widely used in fingerprint-based deep learning researches.
All properties are calculated with RDKit (Landrum et al., ).

### B.4 Non-hierarchical Markov decision process

**Table 6: Comparison with SynFlowNet and RGFN. RGFN investigates scaling up to 64,000 building blocks, but their experimental validation and proof-of-principle implementations use only 350.**
| Methods | Hierarchical | Action Embedding | ReactUni | Num Blocks | Max Reaction Steps |
| --- | --- | --- | --- | --- | --- |
| Enamine REAL | - | - | $\bigcirc$ | $>$1M | 3 |
| SynFlowNet | $\bigcirc$ | $\times$ | $\bigcirc$ | 6,000 | 4 |
| RGFN | $\bigcirc$ | $\bigcirc$ | $\times$ | 350$\sim$64,000 | 4 |
| RxnFlow | $\times$ | $\bigcirc$ | $\bigcirc$ | $>$1M | 3 |

While RxnFlow shares similar rules with SynFlowNet (Cretu et al., ) and RGFN (Koziarski et al., ), there are several major differences as described in Table 6.
In particular, we formulate the non-hierarchical MDP instead of the hierarchical MDP for bi-molecular reactions.

To describe the difference, we supposed the situation that the toxicity is observed in molecules containing a functional group $-\text{NR}_{3}$, thereby excluding blocks with the function group from the building block library.
In trained hierarchical MDP, the modification of the block library cannot change the probability of reaction templates, leading to the overestimation or underestimation of edge flows (the red color of Figure 7(a)).
However, in the non-hierarchical MDP, the exclusion of some reactants does not affect actions that share the same reaction template (Figure 7(b)).

Figure: Figure 7: Illustration of a situation where an additional objective is introduced: excluding building blocks (reactants) containing the functional group $-\text{NR}_{3}$. The Gray dashed lines mean masked actions. (a) GFlowNet on hierarchical MDP. (b) GFlowNet on non-hierarchical MDP.
Refer to caption: extracted/5905292/Figures/appendix/non-hierarchical.png

## Appendix C Experimental details

### C.1 Action space details

#### Reaction templates.

In this work, we used the reaction template set constructed by Cretu et al. () from two public collections Hartenfeller et al. (); Button et al. ().
The entire reaction template set includes 71 reaction templates, 13 for uni-molecular reactions and 58 for bi-molecular reactions.
In in silico reactions with bi-molecular reaction templates, the products depend on the order of the input reactants.
To ensure the consistency of the action, we consider two templates for each bi-molecular template according to the order of reactants, i.e. $|\mathcal{R}_{1}|=13,|\mathcal{R}_{2}|=116$.
We note that our template set does not contain templates that have the same first and second reactant patterns.

#### Building blocks.

We used the Enamine comprehensive catalog with 1,309,385 building blocks released on 2024.06.10 (Grygorenko et al., ).
We filtered out building blocks that are not RDKit-readable, have no possible reactions, or contain unallowable atoms(^2^22The allowable atom types are B, C, N, O, F, P, S, Cl, Br, I.), resulting in 1,193,871 remaining blocks.

### C.2 GFlowNet Training details

To minimize optimization performance influencing factors, we mostly followed the standard GFlowNet’s model architecture and hyperparameters(^3^33Default hyperparameters in [https://github.com/recursionpharma/gflownet/blob/trunk/src/gflownet/tasks/seh_frag.py](https://github.com/recursionpharma/gflownet/blob/trunk/src/gflownet/tasks/seh_frag.py)) except for some parameters in Table 7.
All experiments were performed on a single NVIDIA RTX A4000 GPU.

**Table 7: Default hyperparameters used in RxnFlow training.**
| Hyperparameters | Values |
| --- | --- |
| Minimum trajectory length | 2 (minimum reaction steps: 1) |
| Maximum trajectory length | 4 (maximum reaction steps: 3) |
| GFN temperature $\beta$ | Uniform(0, 64) |
| Train random action probability | 0.05 (5%) |
| Action space subsampling ratio | 1% |
| Building block embedding size | 64 |

For action space subsampling, we randomly subsample 1% actions for AddFirstReactant and each bi-molecular reaction template $r\in\mathcal{R}_{2}$.
However, for bi-molecular reactions with small possible reactant block sets $\mathcal{B}_{r}\in\mathcal{B}$, the memory benefit from the action space subsampling is small while a variance penalty is large.
Therefore, we set the minimum subsampling size to 100 for each bi-molecular reaction, and the action space subsampling is not performed when the number of actions is smaller than 100.

The number of actions for each action type is imbalanced, and the number of reactant blocks ($\mathcal{B}_{r}$) for each bi-molecular reaction template $r$ is also imbalanced.
This can make some rare action categories not being sampled during training.
We empirically found that ReactBi action were only sampled during 20,000 iterations (1.28M samples) in a toy experiment that uses one bi-molecular reaction template and 10,000 building blocks in some random seeds.
Therefore, we set the random action probability as the default of 5%, and the model uniformly samples each action category in the random action sampling.
This prevents incorrect predictions by ensuring that the model experiences trajectories including rare actions.
We note that this random selection is only performed during model training.

### C.3 Experimental details

#### Pocket-specific optimization.

For pocket-specific optimization with GPU-accelerated UniDock, we normalize the Vina and SA scores for multi-objective optimization:

$$ $\widehat{\text{Vina}}(x)=-0.1\max(\text{Vina}(x),0),\quad\widehat{\text{SA}}(x )=(10-\text{SA}(x))/9.$ (30) $$

For FragGFN, we set the maximum trajectory length to 9, and for SynFlowNet and RGFN, we used the same hyperparameters as our framework except for the maximum trajectory length and building block library size.
According to their default setting, we set a the maximum trajectory length to 5 rather than 4, and we randomly sampled 350 building blocks for RGFN and 6000 building blocks for SynFlowNet.
Moreover, we do not perform action subsampling for SynFlowNet and RGFN.

#### Pocket-conditional generation in a zero-shot manner.

We used the modified version of the TacoGFN’s reward function and training set.
Since we don’t need to optimize SA (Ertl & Schuffenhauer, ), we excluded the SA term from the reward functions:

$$ $r_{\text{affinity}}(x)=\begin{cases}0&\text{if}~{}0\leq\text{Proxy}(x)\\ -0.04\times\text{Proxy}(x)&\text{if}~{}-8\leq\text{Proxy}(x)\leq 0\\ -0.2\times\text{Proxy}(x)-1.28&\text{if}~{}-13\leq\text{Proxy}(x)\leq-8\\ 1.32&\text{if}~{}\text{Proxy}(x)\leq-13\\ \end{cases}$ $$

$$ $r_{\text{QED}}(x)=\begin{cases}\text{QED}(x)/0.7&\text{if}~{}\text{QED}(x)\leq 0 .7\\ 1&\text{otherwise}\\ \end{cases}$ $$

$$ $\text{Reward}(x)=\frac{r_{\text{affinity}}(x)\times r_{\text{QED}}(x)}{\sqrt[3 ]{\text{HeavyAtomCounts}(x)}}$ (31) $$

For hyperparameters, we set the pocket embedding dimension to 128 and the training GFN temperature to Uniform(0, 64) which are used in TacoGFN.
We trained the model with 40,000 oracles whereas TacoGFN is trained for 50,000 oracles.

#### Introducing further objectives without retraining.

For the restricted block library (TPSA$<$30), we set the action space subsampling ratio as 10% for both AddFirstReactant and ReactBi, and we set that as 1% for an entire library.

#### Scaling action space without retraining.

For action space subsampling, we set the subsampling ratio as 2% for both AddFirstReactant and ReactBi for the “seen” and “unseen” libraries. For the “all” library (“seen” + “unseen”), we set the subsampling ratio as 1%.

#### Ablation study.

We set different subsampling ratios according to the building block library size.
For 100-sized, 1k-sized, and 10k-sized libraries, we do not perform the action space subsampling.
We set an action space subsampling ratio of 10% for a 100k-sized one and 1% for a 1M-sized one.

### C.4 Softwares

#### Molecular docking software.

For a fair comparison with the baseline model, we used UniDock (Yu et al., ) for target-specific generation and QuickVina 2.1 (Alhossary et al., ) for SBDD.
The initial ligand conformer is generated with srETKDG3 (Wang et al., ) in RDKit (Landrum et al., ).
For QuickVina, we converted the molecule format to pdbqt with OpenBabel (O’Boyle et al., ) and AutoDock Tools (Huey et al., ).
To set up an exhaustive search, we set the search mode to balance for UniDock and the exhaustiveness to 8 for QuickVina.
We kept the seed fixed at 1 throughout the ETKDG and the whole docking process.

#### Docking proxy.

We used the QuickVina proxy proposed by Shen et al. () which is implemented in PharmacoNet (Seo & Kim, ).
We used a proxy model trained on the CrossDocked2020 training set rather than the model trained on the ZINCDock15M training set.

#### Synthetic accessibility estimation.

To evaluate the synthetic accessibility of molecules, we used the retrosynthesis planning tool AiZynthFinder (Genheden et al., ).
AiZynthFinder uses MCTS to find synthesis paths and estimate the number of steps, search time, success rate, and synthetically accessible score as metrics to indicate synthesis complexity or synthesizability.

### C.5 Baselines

#### SynNet.

For SynNet (Gao et al., ), we perform multi-objective optimization with the following reward function:

$$ $R(x)=0.5\text{QED}(x)+0.5\widehat{\text{Vina}}(x).$ (32) $$

According to the standard setting for optimization, we set the number of offspring to 512 and the number of oracles to 125.
To use pre-trained models, we used SynNet’s template set (91 templates) instead of a template set (71 templates) used in our work.

#### RGFN.

Since the code of RGFN is not released, we reimplement the RGFN.

#### BBAR.

Since BBAR (Seo et al., ) allows multi-conditional generation, we directly used QED and docking scores without any processing.
We split 64,000 ZINC20 molecules according to the reported splitting of BBAR: 90% for the training set, 8% for the validation set, and 2% for the test set (in our case, the number of sampling molecules).
We performed UniDock for training and validation set to prepare the label of the molecules.
Since BBAR requires the desired property value, we used the average docking score of the top 100 diverse modes from our model.

#### Pocket2Mol, TargetDiff, DecompDiff, TacoGFN.

We followed the reported generative setting to generate 100 molecules for each CrossDocked test pocket.
We set the center of the pockets with the reference ligands in the CrossDocked2020 database.
We reuse reported runtime in Shen et al. (), which is measured on NVIDIA A100 for Pocket2Mol, TargetDiff, and DecompDiff, and NVIDIA RTX3090 for TacoGFN.

### C.6 LIT-PCBA Pockets

Table 8 describes the protein information used in pocket-specific optimization with UniDock, which is performed on Sec. 4.1.

**Table 8: The basic target information of the LIT-PCBA dataset and PDB entry used in this work.**
| Target | PDB Id | Target name |
| --- | --- | --- |
| ADRB2 | 4ldo | Beta2 adrenoceptor |
| ALDH1 | 5l2m | Aldehyde dehydrogenase 1 |
| ESR_ago | 2p15 | Estrogen receptor $\alpha$ with agonist |
| ESR_antago | 2iok | Estrogen receptor $\alpha$ with antagonist |
| FEN1 | 5fv7 | FLAP Endonuclease 1 |
| GBA | 2v3d | Acid Beta-Glucocerebrosidase |
| IDH1 | 4umx | Isocitrate dehydrogenase 1 |
| KAT2A | 5h86 | Histone acetyltransferase KAT2A |
| MAPK1 | 4zzn | Mitogen-activated protein kinase 1 |
| MTORC1 | 4dri | PPIase domain of FKBP51, Rapamycin |
| OPRK1 | 6b73 | Kappa opioid receptor |
| PKM2 | 4jpg | Pyruvate kinase muscle isoform M1/M2 |
| PPARG | 5y2t | Peroxisome proliferator-activated receptor $\gamma$ |
| TP53 | 3zme | Cellular tumor antigen p53 |
| VDR | 3a2i | Vitamin D receptor |

## Appendix D Additional results

### D.1 Additional results for Pocket-specific generation task

We reported the additional results of Sec. 4.1 for the remaining 10 pockets on the LIT-PCBA benchmark.

**Table 9: Hit ratio (%). Average and standard deviation for 4 runs. The best results are in bold.**
|  |  | Hit ratio (%, $\uparrow$) |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Category | Method | GBA | IDH1 | KAT2A | MAPK1 | MTORC1 |
| Fragment | FragGFN | 5.00 ($\pm$ 4.24) | 4.50 ($\pm$ 1.66) | 1.25 ($\pm$ 0.83) | 0.75 ($\pm$ 0.83) | 0.00 ($\pm$ 0.00) |
| FragGFN+SA | 3.00 ($\pm$ 1.00) | 4.50 ($\pm$ 4.97) | 1.50 ($\pm$ 0.50) | 2.00 ($\pm$ 1.73) | 0.00 ($\pm$ 0.00) |  |
| Reaction | SynNet | 50.00 ($\pm$ 0.00) | 50.00 ($\pm$ 0.00) | 29.17 ($\pm$ 18.16) | 37.50 ($\pm$ 21.65) | 0.00 ($\pm$ 0.00) |
| BBAR | 17.75 ($\pm$ 2.28) | 19.50 ($\pm$ 1.50) | 18.75 ($\pm$ 1.92) | 16.25 ($\pm$ 3.49) | 0.00 ($\pm$ 0.00) |  |
| SynFlowNet | 58.00 ($\pm$ 4.64) | 59.00 ($\pm$ 4.06) | 55.50 ($\pm$ 10.23) | 47.25 ($\pm$ 6.61) | 0.00 ($\pm$ 0.00) |  |
| RGFN | 48.00 ($\pm$ 1.22) | 43.00 ($\pm$ 2.74) | 49.00 ($\pm$ 1.22) | 38.00 ($\pm$ 4.12) | 0.00 ($\pm$ 0.00) |  |
| RxnFlow | 66.00 ($\pm$ 1.58) | 64.00 ($\pm$ 5.05) | 66.50 ($\pm$ 2.06) | 63.00 ($\pm$ 4.64) | 0.00 ($\pm$ 0.00) |  |
|  |  | OPRK1 | PKM2 | PPARG | TP53 | VDR |
| Fragment | FragGFN | 0.50 ($\pm$ 0.50) | 7.25 ($\pm$ 1.92) | 0.75 ($\pm$ 0.43) | 4.25 ($\pm$ 1.64) | 0.00 ($\pm$ 0.00) |
| FragGFN+SA | 0.50 ($\pm$ 0.87) | 4.50 ($\pm$ 1.50) | 1.00 ($\pm$ 0.71) | 2.25 ($\pm$ 1.92) | 0.00 ($\pm$ 0.00) |  |
| Reaction | SynNet | 0.00 ($\pm$ 0.00) | 0.00 ($\pm$ 0.00) | 33.33 ($\pm$ 20.41) | 8.33 ($\pm$ 14.43) | 0.00 ($\pm$ 0.00) |
| BBAR | 2.50 ($\pm$ 1.12) | 20.00 ($\pm$ 0.71) | 10.50 ($\pm$ 2.69) | 14.00 ($\pm$ 3.94) | 0.00 ($\pm$ 0.00) |  |
| SynFlowNet | 23.50 ($\pm$ 5.94) | 50.75 ($\pm$ 1.09) | 53.50 ($\pm$ 5.68) | 55.50 ($\pm$ 9.94) | 0.00 ($\pm$ 0.00) |  |
| RGFN | 2.50 ($\pm$ 2.06) | 34.75 ($\pm$ 6.57) | 29.00 ($\pm$ 6.52) | 37.00 ($\pm$ 6.60) | 0.00 ($\pm$ 0.00) |  |
| RxnFlow | 72.25 ($\pm$ 2.05) | 62.00 ($\pm$ 3.24) | 65.50 ($\pm$ 4.03) | 67.50 ($\pm$ 2.96) | 1.75 ($\pm$ 0.83) |  |

**Table 11: Synthesizability. Average and standard deviation for 4 runs. The best results are in bold.**
|  |  | Percentage of Synthesizable Molecules (%, $\uparrow$) |  |  |  |  |
| --- | --- | --- | --- | --- | --- | --- |
| Category | Method | GBA | IDH1 | KAT2A | MAPK1 | MTORC1 |
| Fragment | FragGFN | 5.00 ($\pm$ 4.24) | 4.50 ($\pm$ 1.66) | 1.25 ($\pm$ 0.83) | 0.75 ($\pm$ 0.83) | 2.75 ($\pm$ 1.30) |
| FragGFN+SA | 3.00 ($\pm$ 1.00) | 4.50 ($\pm$ 4.97) | 1.50 ($\pm$ 0.50) | 3.25 ($\pm$ 1.48) | 3.50 ($\pm$ 2.50) |  |
| Reaction | SynNet | 50.00 ($\pm$ 0.00) | 50.00 ($\pm$ 0.00) | 45.83 ($\pm$ 27.32) | 50.00 ($\pm$ 0.00) | 54.17 ($\pm$ 7.22) |
| BBAR | 17.75 ($\pm$ 2.28) | 19.50 ($\pm$ 1.50) | 18.75 ($\pm$ 1.92) | 16.25 ($\pm$ 3.49) | 18.75 ($\pm$ 3.90) |  |
| SynFlowNet | 58.00 ($\pm$ 4.64) | 59.00 ($\pm$ 4.06) | 55.50 ($\pm$ 10.23) | 47.25 ($\pm$ 6.61) | 57.00 ($\pm$ 7.58) |  |
| RGFN | 48.00 ($\pm$ 1.22) | 43.00 ($\pm$ 2.74) | 49.00 ($\pm$ 1.22) | 42.00 ($\pm$ 3.00) | 44.50 ($\pm$ 4.03) |  |
| RxnFlow | 66.00 ($\pm$ 1.58) | 64.00 ($\pm$ 5.05) | 66.50 ($\pm$ 2.06) | 63.00 ($\pm$ 4.64) | 70.50 ($\pm$ 2.87) |  |
|  |  | OPRK1 | PKM2 | PPARG | TP53 | VDR |
| Fragment | FragGFN | 2.50 ($\pm$ 2.29) | 8.75 ($\pm$ 3.11) | 0.75 ($\pm$ 0.43) | 4.25 ($\pm$ 1.64) | 3.50 ($\pm$ 2.18) |
| FragGFN+SA | 3.25 ($\pm$ 1.79) | 9.75 ($\pm$ 2.28) | 1.25 ($\pm$ 1.09) | 2.25 ($\pm$ 1.92) | 3.75 ($\pm$ 2.77) |  |
| Reaction | SynNet | 54.17 ($\pm$ 7.22) | 50.00 ($\pm$ 0.00) | 54.17 ($\pm$ 7.22) | 29.17 ($\pm$ 18.16) | 45.83 ($\pm$ 7.22) |
| BBAR | 13.75 ($\pm$ 3.11) | 20.00 ($\pm$ 0.71) | 15.50 ($\pm$ 2.29) | 18.50 ($\pm$ 3.28) | 12.25 ($\pm$ 3.34) |  |
| SynFlowNet | 56.50 ($\pm$ 7.63) | 50.75 ($\pm$ 1.09) | 53.50 ($\pm$ 5.68) | 55.50 ($\pm$ 9.94) | 53.50 ($\pm$ 1.80) |  |
| RGFN | 48.00 ($\pm$ 2.55) | 48.50 ($\pm$ 3.20) | 47.00 ($\pm$ 5.83) | 53.25 ($\pm$ 3.63) | 46.50 ($\pm$ 2.69) |  |
| RxnFlow | 72.25 ($\pm$ 2.05) | 62.00 ($\pm$ 3.24) | 65.50 ($\pm$ 4.03) | 67.50 ($\pm$ 2.96) | 66.75 ($\pm$ 2.28) |  |

### D.2 Additional results for scaling action space without retraining

We reported the additional results for additional reward exponent settings ($R^{\beta}$).

Figure: Figure 8: QED reward distribution of generated molecules.
Refer to caption: extracted/5905292/Figures/result/ana1-ext.png

### D.3 Theoretical analysis

To assess the impact of action space subsampling in GFlowNet training, we conduct a toy experiment using a simplified setup with 10,000 blocks, one uni-molecular reaction template, one bi-molecular reaction template, and the QED objective.
We used the Hell-Volhard-Zelinsky reaction as a uni-molecular reaction and the Amide reaction as a bi-molecular reaction, which are illustrated in Fig. 9.
We used a minimum trajectory length of 1, max trajectory length of 2, constant GFN temperature of 1.0, and learning rate decay of 3,000 for $P_{F}$ and $\log Z$.
For the GFlowNet sampler, we used the same weights of the proxy model, i.e. EMA factor of 0.
We performed optimization for 30,000 oracles with a batch size of 64.

As shown in Figures 10(a) and 10(b), we compare a baseline GFlowNet trained without subsampling (“base”) to models using various subsampling ratios and Monte Carlo (MC) sampling.
The differences in $\log Z_{\theta}$ are relatively small ($<$0.005) across all settings, and increasing MC samples for state flow estimation $F_{\theta}$ further reduced the bias.
In Figures 10(c) and 10(d), we also evaluate the bias in the trajectory balance loss ($\hat{\mathcal{L}}_{\text{TB}}$) and its gradient norm ($\|\nabla_{\theta}\hat{\mathcal{L}}_{\text{TB}}\|$) during training, finding negligible differences compared to the true values.
These results indicate that our importance sampling reweighting approach effectively mitigates bias from action space subsampling, enabling efficient and accurate policy estimation.

Figure: Figure 9: Reaction templates employed in toy experiments. (a) Hell-Volhard-Zelinsky reaction. (b) Amide reaction.
Refer to caption: extracted/5905292/Figures/appendix/toy-templates.png

Figure: Figure 10: Bias estimation. (a) $\log Z_{\theta}$ according to the action space subsampling ratio (left) and the number of MC samples where the subsampling ratio is 1/9 (right). (b) The trajectory balance loss ($\mathcal{L}_{\text{TB}},\hat{\mathcal{L}}_{\text{TB}}$) where the subsampling ratio is 1/9 under 4 MC samples. (c) The loss gradient norms ($\|\nabla_{\theta}\mathcal{L}_{\text{TB}}\|,\|\nabla_{\theta}\hat{\mathcal{L}}_ {\text{TB}}\|$) where the subsampling ratio is 1/9 under 4 MC samples.
Refer to caption: extracted/5905292/Figures/result/ana4_toy.png

## Peer Reviews

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: - Table 5: Values in time column have ‘a’ and ‘b’ superscripts. What do these mean?
- Section 4.2: Can the authors clarify what is ‘zero-shot’ about this method?
- Do you use the same reward functions for methods when comparing to previous GFlowNets?
soundness: 2
strengths: - RxnFlow performs well across all tasks/targets and the comparison of a reaction-based GFlowNet to SBDD models is welcome.
- The subspace sampling technique substantially reduces memory and computational complexity.
- The use of non-hierarchical action base combining reactions and build blocks is novel and well motivated.
- Article is generally very well written and the quality of figures is high.
summary: This paper introduces RXNFLOW, a generative framework that integrates synthesizability considerations into molecular generation for drug design. Following previous work, RXNFLOW addresses synthesizability by generating molecules using molecular building blocks, chemical reaction templates and a GFlowNet model. RXNFLOW aims to efficiently learn to sample from a large and complex action space, employing an action space subsampling technique to handle the combinatorial explosion of possible synthetic pathways. The authors claim that RXNFLOW achieves strong performance in pocket-specific and pocket-conditional generation tasks, with high synthesizability and diversity scores compared to previous models.
weaknesses: **Weakness/points to be worked on:**

- The main contribution of the works is the non-hierarchical and continuous action space. There are many theoretical benefits to this, but the benefit of their method is not concretely controlled with respect to the building block data used and compute budgets.
    - I would recommend Figure 6 be amended with scaling laws for the SynFlowNet [1], RGFN [2] and SyntheMol [3] methods.
- Another fundamental claim is that the model can generalise to ‘’unseen” building blocks. While in theory their method should allow this, it has not been justified with experimental evidence.
    - As far as I can tell, this is only benchmarked in Figure 5 which looks at QED. This is a highly uninformative metric given how well generative models do on QED. Benchmarking the seen and unseen molecules on one of the targets for instance would be a better example.
    - Furthermore, a random split of the building blocks are used, given the high degree of redundancy in the Enamine Building Block set it is not surprising that performance does not change as these sets are highly similar. How does the model do on highly novel building blocks?
- While technically impressive, I see no practical utility gained in the case on not retraining on new building block sets as (i) these libraries at not updated regularly and (ii) the costs of training these models are so small one would prefer to retrain a model anyways.
    - The y-axis scales in Figure 6 are extremely narrow for many metrics and does not represent a meaningful change in values for Vina and Structural Diversity. For Figure 6b there is only a meaningful reduction in scaffold diversity when number of available building blocks is 10^2, an extremely smaller number and other works do not use.
    - Furthermore, RGFN [2] introduces a similar method to scale to large build block spaces and a proper comparison is not conducted.
- Some aspects of the evaluation could be improved:
    - What is the with per run variance for the docking scores for molecules in Tables 1-4? Is there really a significant difference between the methods? I would recommend including per run box plots of Vina scores in the appendix.
    - There is next to no discussion on the limitations of their model and possible future work
    - L458: The authors do not include statistical tests or error bars/standard deviation values in Table 5 yet claim ”RXNFLow achieves significant improvements in drug-related properties” in the CrossDocked experiments.
    - No effort is made to disentangle the technical changes of the authors model v.s. building block library size when comparing to baselines.
- The article is generally very well written but parts are not clear, for instance:
    - The fact that the continuous embedding space is based on chemical fingerprints is not mentioned in the main text.
    - Table 6:  No legend/keys are provided to indicate what the symbols mean. 
- L429: “Since RXNFLOW explicitly considers synthesiz- ability, we exclude the SA score from the TacoGFN’s reward function as described in Sec. C.2.”.  I Could not find reference to SA score in Section C.2. Am I right in thinking the authors removed the SA reward from TacoGFN? In which case they ablated another model they are comparing to. Or have they removed it when training their model? Please clarify.
- No code is provided.

**Minor points:**
- L122: “...trained models with the GFlowNet objective”, is it not more accurate to say these works propose reaction-based GFlowNet models?
- I struggle to see the added value in Algorithm 1 in the main text. It does not aid in understanding the method. I recommend removing.
- Table 6 seems wrong? 
- In Table 3, calling this metric “Percentage of Synthesizable Molecules” is misleading, as these are predicted not ground truth. I would prefer AIZynFinder success as done in the previous literature [1].
- I appreciate the use of the LIT-PCBA targets, but is there a reason the SEH Proxy used in the rest of the GFlowNet literature [1,2] was not used?
- Section B.1: The formulation sees to be highly inspired by SynFlowNet, if so please cite accordingly.




1. Cretu, Miruna, et al. "Synflownet: Towards molecule design with guaranteed synthesis pathways." arXiv preprint arXiv:2405.01155 (2024).
2. Koziarski, Michał, et al. "RGFN: Synthesizable Molecular Generation Using GFlowNets." arXiv preprint arXiv:2406.08506 (2024).
3. Swanson, Kyle, et al. "Generative AI for designing and validating easily synthesizable and structurally novel antibiotics." Nature Machine Intelligence 6.3 (2024): 338-353.

---

confidence: 4
contribution: 2
initial_rating: 5
presentation: 3
questions: Refer to above
soundness: 3
strengths: 1. The motivation of the paper is clear. I believe that apart from the generative approaches that formulate the problem as a constraint generation/projection, the proposed methods focus on an alternative perspective, i,e,  explicitly limiting the action space of gflownets, which should also be explored.

2. The paper introduces a simple yet effective approach, which is referred to as subspace sampling. The method takes a very simple formulation while it enables the reduction of complexity and enables feasible training of the algorithm. 

3. The empirical study of the proposed framework is extensive, with both pocket conditioned/ specific ligand generation being considered.
summary: This paper introduces a new gflownets-based method RxnFlow for conducting the synthesizability-aware generation. To this end, the authors proposed to redefine the action space with Reaction templates and building blocks. To make the training procedure flexible, the authors involve a subsampling method for training RxnFlow in huge spaces.  The experiments over the SBDD benchmark have demonstrated the effectiveness of the RxnFlow for generating molecules with both good binding affinity and synthesizability.
weaknesses: 1. Though I appreciate the methods with simplicity and effectiveness, I believe that a more systematic investigation and overview of the proposed methods is needed.  Based on the bias/variance discussion in the appendix, does the proposed approach conduct a variance/efficiency tradeoff, i.e. large variance for high efficiency? This is a little counterintuitive for me, could the authors discuss this further. From the ablation in Fig. 6, with sufficient steps, a larger subspace shows better performance. Does the key benefit of the proposed method is accelerating the optimization?

2. Some important baselines are suggested to be included in Table 5, for example [1] for atom-level SBDD. And the paper mentioned the recent relevant work [2], so why not include it? 

3. I would like to suggest the author include the SA score as another metric for evaluating the synthesizability. I am curious to see whether SA is aligned with the used metrics in the paper. 

I would like to adjust my score based on the response.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. What limitations do you perceive in the current method, and how might it be improved in the future?
2. Regarding the pocket-conditional generation task, I am curious about the method's ability to synthesize molecules with specific binding characteristics, such as calculating the 'delta score' metric propsed in Paper [1].


[1]. Gao, Bowen, et al. "Rethinking Specificity in SBDD: Leveraging Delta Score and Energy-Guided Diffusion." ICML 2024.
soundness: 3
strengths: 1. The method enables the generation of synthetic pathways for molecules, allowing for the sampling of highly synthesizable compounds while maintaining a significant level of diversity, which is meaningful for drug discovery.
2. With the enhancement of building blocks, the method demonstrates good scalability.
3. The experiments were conducted thoroughly, and the presentation is relatively clear.
summary: This paper proposes  a GFlowNets based method 'RXNFLOW' for drug discovery. By employing space subsampling technique, RXNFLOW can expand the search space and handle massive action spaces. This method can be adopted for pocket-specific optimization task as well as pocket-conditional generation task.
weaknesses: 1. Regarding the pocket-conditional generation task, to my knowledge, the more advanced methods [1] have not been compared.  This has somewhat affected the persuasive power of the experiments.
2. Compared to other SBDD methods, it seems that direct generation of conformations combined with the pocket is not achievable.
3. Regarding the pocket-specific optimization task, I notice that the reward function consits of Vina Score which is also used for evaluation. I have concerns about whether the method may be overfitting to Vina, which could result in inflated evaluation metrics.


[1]. Qu, Yanru, et al. "MolCRAFT: Structure-Based Drug Design in Continuous Parameter Space." ICML 2024.

---

confidence: 4
contribution: 3
initial_rating: 5
presentation: 3
questions: 1. For conventional SBDD methods, they are generally trained on molecules that can bind to protein pockets, resulting in generated small molecules with a specific distribution. I am curious about how the molecules generated in this space of synthesizable reaction-based small molecules differ from those generated by SBDD methods, such as pocket2mol. Because drug-likeness is not solely determined by synthesizability and vina score alone. For instance, the bond angle, number of rings. These facts are hard to encode into the reward function.
2. Please compare the speed of different methods. 					
3. Why use 3 reaction steps for your method and 4 steps for SynFlowNet and RGFN?
soundness: 3
strengths: Overall, this paper is interesting. The GFlowNet model is suitable to be adopted to solve the synthesizability problem of small molecules. The problem with GFlowNet is that it lacks a method to construct a graph that covers the entire space of small molecules. This paper cleverly utilizes chemical synthesis reactions to construct this graph, while also adopting a down-sampling approach to enhance the effectiveness of the method.
summary: In this manuscript, the authors adopted the GFlowNet to generate synthesizable molecules based on the protein pocket. They were able to demonstrate that RXNFLOW outperforms existing reaction-based and fragment-based models in pocket-specific synthetic pathways, and their method achieves state-of-the-art performance on CrossDocked2020 for pocket-conditional generation.
weaknesses: The main issue is the insufficient characterization of 3D interactions; the model does not explicitly model the interactions between atoms in 3D space and only uses the generated molecule's vina score with the pocket as the reward. In addition, a careful discussion of the differences from SynFlowNet is needed. Currently, the authors have avoided mentioning SynFlowNet in the introduction. According to my understanding, the method proposed in this paper has no structural difference from SynFlowNet in terms of the model framework, with some improvements in sampling techniques. I hope the authors can discuss the differences from SynFlowNet in detail from the perspectives of graph construction for the flow, reward settings, and the sampling algorithm.

---

confidence: 3
contribution: 2
initial_rating: 5
presentation: 2
questions: 1. Table 5, what is a and b?
2. line 45, what does it mean by 'chemical modification can degrade in the optimized propeties', can you give some examples?
3. It's unclear to me how you introduce the additional objective without retraining, can you clarity that?
soundness: 2
strengths: 1. The proposed methods outperforms baseline in standard benchmark.
2. The methods exhibit nice scaling behavior on unseen blocks.
summary: The author proposed a GFlowNet for on synthetic pathway, aiming to leverage the large action space for drug synthesis while retaining computational efficiency. The author proposed an action embedding which embeds the template and building blocks to parameterize the policy network. They further reduces the computational complexity by introducing a subsampling technique on the action space via importance sampling. The author validate their approach on LIT-PCBA, CrossDocked2020 dataset and outperformed baselines.
weaknesses: 1. **Conceptual Novelty:** Utilizing reinforcement learning (RL) to determine reaction pathways is not a novel approach. The authors have simply applied this existing methodology to GFlowNets, which does not represent a significant conceptual advancement.

2. **Technical Novelty:** The proposed action embedding and importance sampling techniques do not appear to introduce new concepts. Essentially, the authors are reparameterizing the policy network with template and building block embeddings— a standard practice in RL when addressing large action spaces.

3. **Figure Clarity:** Figure 2 is currently confusing and requires further refinement. The use of diamonds to represent templates and squares for reactions may initially cause some misunderstanding. In contrast, Figure 7 is clearer and more effectively conveys the intended information.

4. **Efficiency Discussion:** The paper's main claim centers on reducing the action space to enhance computational tractability. However, there is a need for a more in-depth discussion and comparison of the training and inference times between the proposed network and existing baselines to substantiate this claim.

5. **Ablation Study:** The manuscript lacks an ablation study comparing Hierarchical Markov Decision Processes (MDPs) with non-hierarchical MDPs. Such a study is essential to evaluate the effectiveness and advantages of the non-hierarchical approach presented.

Meta-review: The paper proposes RxnFlow, a model based on generative flow networks that is shown to scale to large compound spaces. 

The main contribution of the work is in scaling up GFlowNets to handle very realistic and comprehensive molecular spaces defined by 71 reaction templates and 1.2 million building blocks. Beyond its class of models, I understand that relatively few generative methods were shown to effectively scale to such space. RxnFlow is able to find (in its space) highly active molecules (e.g. comparable to MolCRAFT) that are significantly more synthetically accessible (as judged by open-source retrosynthesis software).

This feat is achieved using very well executes techniques such as adding action embeddings. Concurrently, SynFlowNet paper has been updated using related techniques.

Three reviewers voted for marginally accepting the paper, while two maintained their marginally negative opinion of the paper. Reviewers have generally engaged in the rebuttal, in certain cases increasing their score. 

The strong empirical results and the overall soundness of the approach warrant acceptance of the work. The paper opens a practical avenue for building de novo generative models on large building block spaces, enriching the practical toolbox in the early stage drug discovery. Thank you for your submission and it is my pleasure to recommend accepting the paper.

Decision: Accept (Poster)