MolHF: A Hierarchical Normalizing Flow for Molecular Graph Generation

1 Introduction

Designing novel and valuable molecular structures is a critical yet challenging task with great application potential in drug and material discovery. Owing to the enormously vast chemical space with discrete molecular structures, whose scale was estimated to be on the order of 1033 [Polishchuk et al., 2013], it is time-consuming and costly to search exhaustively in this space. Recently, there has been a surge of interest in developing generative models for graphs to resolve this issue by exploring the chemical space in a data-driven manner and sampling novel molecules [Jin et al., 2018].

∗Corresponding authors.

Much prior work represents molecules as attributed graphs and leverages a variety of deep generative frameworks for molecular graph generation, including variational autoencoders (VAEs) [Simonovsky and Komodakis, 2018], generative adversarial networks (GANs) [De Cao and Kipf, 2018], and normalizing flows [Shi et al., 2020; Zang and Wang, 2020]. However, most of them ignore the inherent hierarchical structure (with rich semantic information) of the molecular graphs. The significance of exploiting hierarchical structure lies in the fact that some particular substructures, such as molecular fragments and functional groups, appear frequently in molecules and are likely to determine molecular properties. Thereby, these models are sub-optimal for molecular generation and property optimization. Additionally, molecular graph generation is prone to errors, e.g., randomly removing or adding a covalent bond probably yields an invalid molecular structure. Thus, whilst existing models have demonstrated promising performance for small molecules, they tend to perform dreadfully for large molecules or even fail to effectively generate more complex molecules of larger size [Jin et al., 2020]. Recently, substructure-based methods [Jin et al., 2020], where the substructures are manually extracted beforehand and leveraged as building blocks, have been proposed to exploit the implicit hierarchical structure. Nevertheless, the substructure extraction strategies are mostly based on heuristics and prior domain knowledge [Jin et al., 2018; Jin et al., 2020], and the extraction process is cumbersome (e.g., counting the frequency of occurrence of all possible substructures). These practical drawbacks limit the flexibility and scalability of such methods.

Advances in hierarchical generative frameworks have made significant progress in high-resolution image generation [Karras et al., 2018; Yu et al., 2020] and superresolution [Liang et al., 2021], with the merits of better quality, scalability, and stability. Their key motivation is to progressively generate higher resolution images based on lower resolution input. However, applying hierarchical generative frameworks to derive molecular graphs is non-trivial. The primary challenge unique to hierarchical graph generation is generating intermediate discrete coarsened graph structures, which leads to the inability to backpropagate through samples. Currently, hierarchical molecular graph generation is still in its infancy, and it remains unclear how exploiting the hierarchical structure of molecular graphs could enhance the performance of generative models.

To overcome the above issues, we cast the tricky problem of hierarchical generation over discrete spaces as the reverse process of hierarchical representation learning resorting to invertible normalizing flows [Dinh et al., 2014; Dinh et al., 2016], and propose a new hierarchical flow-based generative model called MolHF, which generates molecular graphs in a coarse-to-fine manner. By defining an invertible transformation to capture the hierarchical structural information in the inference (encoding) process, MolHF can also exploit it in the reverse generation (decoding) process. Specifically, we propose CC-Glow, a Criss-Cross network [Huang et al., 2019] based multi-scale Glow [Kingma and Dhariwal, 2018] which offers a better relational inductive bias, to generate bonds at first, and a coarsened graph conditional flow to generate atoms by performing graph convolution [Kipf and Welling, 2017; Schlichtkrull et al., 2018] on the pre-generated coarsened graph structure (bonds) of each scale. We evaluate our MolHF in aspects of random generation and property optimization on the standard ZINC250K dataset [Irwin et al., 2012]. Further, we use the Polymer dataset [St. John et al., 2019] to verify that MolHF is capable of modeling larger molecules that are impractical with prior flow-based methods. As demonstrated with extensive experiments, MolHF achieves state-of-the-art performance, indicating the effectiveness of the proposed hierarchical architecture. Finally, we analyze the generative mechanism of MolHF through hierarchical visualization. Our key contributions are summarized as follows.

• We propose MolHF for hierarchical molecular graph generation. MolHF is the first flow-based model that can exploit the hierarchical structural information.
• We cast the non-trivial hierarchical generation problem over discrete spaces as the reverse process of hierarchical representation learning to sidestep the nondifferentiable issue.
• MolHF outperforms state-of-the-art baselines in random generation and property optimization. Meanwhile, experiments on the Polymer dataset show that MolHF is the first flow-based model that can model larger molecules with more than 100 heavy atoms.

2 RelatedWork

2.1 Normalizing Flows

Normalizing flows are a class of probabilistic generative models that aim to learn invertible and differentiable mappings between complex target distributions and simple prior distributions. Compared with the prevailing generative models, such as VAEs and GANs, normalizing flows possess the merits of exact latent variable inference and likelihood estimation, efficient sampling, and stable training [Kingma and Dhariwal, 2018]. Advances in normalizing flows have made impressive progress in various applications, such as density estimation [Dinh et al., 2016], image generation [Kingma and Dhariwal, 2018], and variational inference [Rezende and Mohamed, 2015]. Recently, normalizing flows also found their way into generating more valid and novel molecular graphs [Shi et al., 2020]. Furthermore, much work has been devoted to improving normalizing flows, ranging from expressive power [Ho et al., 2019] to numerical stability [Behrmann et al., 2021]. For a more comprehensive review of normalizing flows, we refer the readers to [Papamakarios et al., 2021].

2.2 Molecular Graph Generation

Over the past few years, there has been a surge in generating molecular graphs resorting to various deep generative frameworks, e.g., VAEs [Kingma and Welling, 2013], GANs [Goodfellow et al., 2014], and autoregressive models. For example, JT-VAE [Jin et al., 2018] first generates a junction tree that represents the scaffold over chemical substructures, then assembles them into a valid molecular graph. MolGAN [De Cao and Kipf, 2018] is the first GAN-based model but shows poor performance in generating unique molecules. GCPN [You et al., 2018] formulates molecular graph generation as a sequential decision process generating the atoms and bonds alternately and leverages reinforcement learning to optimize the generated molecules.

There are also several flow-based models for molecular graph generation, which can be mainly grouped into two categories, i.e., autoregressive methods and one-shot methods. GraphAF [Shi et al., 2020] and GraphDF [Luo et al., 2021] adapt autoregressive flows [Papamakarios et al., 2017] to generate atoms and bonds alternately in a sequential manner. GraphNVP [Madhawa et al., 2019] and MoFlow [Zang and Wang, 2020] are representative one-shot methods, which generate atoms conditioned on the pre-generated bonds in a twostep manner. Our MolHF follows the latter approach to exploit the coarsened graph structure for hierarchical generation while ensuring reversibility. The most related to our method is MolGrow [Kuznetsov and Polykovskiy, 2021], which starts with a single node and splits every node into two recursively to generate the entire molecule. Our MolHF differs from MolGrow in two main aspects: 1) MolGrow transforms the atom matrix and bond tensor simultaneously, therefore the transformed bond tensor becomes continuous hidden states which cannot provide structural information for modeling atoms in the coarsened scale. In contrast, MolHF is able to perform graph convolution to capture the hierarchical structural information, resorting to the aforementioned two-step approach. 2) MolGrow is less flexible and inefficient as the number of atoms in the training molecular graphs must be padded to be a power of 2, leading to extra computational costs. In contrast, we propose a graph coarsening operation to merge any number of nodes and compute the coarsened graph structure, making MolHF adaptable to various molecular datasets by setting appropriate coarsening ratios.

3 Problem Statement and Background

3.1 Hierarchical Molecular Graph Generation

The hierarchical molecular graph generation problem seeks to design a network (e.g., the decoder of VAE or generator of GAN) to generate the entire molecular graph in a coarse-to-fine manner. Put differently, the molecular graph is refined with increasingly fine details, conditioned on the coarser structure. Unfortunately, applying hierarchical generative frameworks to derive molecular graphs is non-trivial, since graphs are discrete objects. The primary challenge is generating intermediate discrete coarsened graph structures, which leads to the inability to backpropagate through samples. To sidestep this non-differentiable issue, we resort to normalizing flows [Dinh et al., 2014] which construct a generative model by defining an invertible transformation. From the perspective of normalizing flows, the challenging problem of hierarchical molecular graph generation can be regarded as the reverse process of hierarchical representation learning.

For a molecular graph G = (X, A) with n atoms, we represent it with an atom matrix X ∈ {0, 1}n×d and a bond tensor A ∈ {0, 1}b×n×n , where d and b are the number of different types of atoms and bonds, respectively. The element Xi,: and A:,i,j are one-hot vectors indicating the type of the ith atom and the type of the bond between the ith atom and j th atom, respectively. Our ultimate purpose is to learn the distribution P(G) via a hierarchical probabilistic generative model.

3.2 Normalizing Flows

The flow-based generative models define a parameterized invertible transformation fθ : X → Z from observational data x ∼ PX (x) to latent variables z ∼ PZ(z), where PX (x) is a complex target distribution and PZ(z) is a simple base distribution (e.g., standard isotropic Gaussian distribution). Under the change-of-variables formula, the exact log-likelihood of x is tractable:

where det ∂fθ(x) ∂x is the Jacobian determinant of fθ at x. Thereby, unlike VAEs, the flow-based models can be trained directly with the maximum likelihood estimation. As for generation, a sample is efficiently generated via the inverse transformation x = f -1 θ (z) where z is sampled from PZ(z).

In practice, fθ = f1 ◦ · · · ◦ fL is composed of a sequence of invertible functions whose Jacobian determinant is convenient to calculate. Herein, we consider using the affine coupling layer, actnorm, and invertible permutation.

Affine coupling layer. One of the most widely used invertible functions is the affine coupling layer [Dinh et al., 2016]. Formally, given x ∈ RD and d < D, the output y ∈ RD can be computed as follows:

where ⊙ stands for Hadamard product, s and t are scale and translation functions, respectively, and r is a rescale function that is commonly instantiated as an exponential function.

Since the log-determinant of the triangular Jacobian of this transformation is efficiently computed as ∑j log r(s(x1:d))j , s and t can be arbitrarily complex functions to improve the expressive power [Dinh et al., 2016]. To mitigate numerical instability, we adopt the sigmoid function following the Swish activation function [Ramachandran et al., 2017] as r with a restricted range, whose effect was thoroughly investigated in [Behrmann et al., 2021].

Actnorm. Actnorm [Kingma and Dhariwal, 2018] is proposed to not only alleviate the problems associated with the numerical instability but also sidestep the noise introduced by the batch normalization. The main idea of actnorm is to normalize the activations to zero mean and unit variance before starting training. Specifically, the parameters of actnorm are initialized with the per-channel mean and standard deviation of the first batch of data, and optimized as regular parameters. As for the feature matrix, actnorm is extended to normalize the feature dimension. For consistency, we continue to call this variant actnorm2D [Zang and Wang, 2020].

Invertible permutation. Invertible 1 × 1 convolution with weight matrix W ∈ Rc×c is a learnable permutation proposed to replace the fixed random permutation [Kingma and Dhariwal, 2018]. With 1 × 1 convolution, all channels are allowed to influence one another, and the expressive power of the following affine coupling layers is subsequently improved. Different from 3D tensor (images), the feature matrix can be permuted by directly multiplying W. To reduce the cost of computing log-determinant of W, we formulate it with LU decomposition [Kingma and Dhariwal, 2018].

4 Methods

In this section, we first introduce the overview of our proposed MolHF, then we elucidate the specifics of MolHF. Lastly, we describe the generation process of MolHF.

4.1 Overview of MolHF

Similar to existing one-shot flow-based models such as GraphNVP [Madhawa et al., 2019] and MoFlow [Zang and Wang, 2020], MolHF aims to generate molecular graphs in one pass through the model, but in a coarse-to-fine manner. To exploit the coarsened graph structure for hierarchical generation while ensuring the reversibility of MolHF, we conditionally factorize the distribution P (G) according to the combinatorial structure of molecular graphs, giving

Based on this factorization, MolHF wishes to hierarchically transform the atom matrix X and bond tensor A respectively. As illustrated in Figure 1, MolHF consists of two flows, namely a coarsened graph conditional flow fX|A for atoms and a CC-Glow fA for bonds. We adopt the well-established dequantization technique [Shi et al., 2020; Zang and Wang, 2020] to convert both X and A into continuous data by adding uniform noise, since discrete data do not fit into the change-of-variables formula. Under the maximum likelihood principle, the parameters of MolHF are updated by gradient descent to minimize the negative log-likelihoods logP(G) over all training molecular graphs:

4.2 Coarsened Graph Conditional Flow

Conditioning on the bond tensor A, coarsened graph conditional flow fX|A aims to learn P(X|A) by transforming the atom matrix X into latent variables zX|A = fX|A(X, A) in a hierarchical manner. Inspired by this motivation, we propose a graph coarsening operation and corresponding coarsened graph coupling layers to extend the multi-scale architecture proposed in [Dinh et al., 2016].

Multi-scale architecture. At each scale, we perform three successive operations. Firstly, a graph coarsening operation is applied to a fine-grained molecular graph, except at the finest scale. Secondly, several steps of transformation which contains actnorm2D followed by an invertible permutation, followed by a coarsened graph coupling layer, are utilized to extract the coarse-grained features. Thirdly, half the features are split into latent variables, except at the coarsest scale. This procedure can be defined recursively,

where h0 = X and A0 = A. Note that the correlations between latent variables of different scales are modeled by assuming ziX|A ∼ N (µi(hi+1 ), ∑i (hi+1 )).

Graph coarsening operation. We merge every ki-1 nodes to construct the coarsened graph, where ki-1 is a hyperparameter controlling the degree of coarsening at each scale. The merged node features are obtained with concatenation, while the coarsened graph structure Ai is computed as:

where S ∈ {0, 1}ni-1× ni-1 ki-1 is a predefined hard assignment matrix assigning every ki-1 nodes as a cluster in order, ni-1 is the number of nodes in the coarsened graph at the i -1th scale, and Ai indicates the connectivity between each pair of clusters [Ying et al., 2018]. As a notorious problem in graph generation, it is tough to define the most suitable atom ordering. We use breadth-first search ordering to preserve spatial locality, which is vital for hierarchical modeling.

Coarsened graph coupling layers. We transform the atom matrix via the affine coupling layers. To achieve permutation equivariance at each scale, we split atom features into two parts along the feature (column) dimension and adopt RGCNs [Schlichtkrull et al., 2018] in the coupling networks, leveraging the coarsened graph structure. The convolution operation of R-GCN can be written as follows:

where Hl is the node embeddings in lth layer, Ai is the ith slice of adjacency tensor, Djj = ∑b,i Ab,i,j, and {W}bi=0 are weight matrices. Compared with the graph coupling layers used in MoFlow [Zang and Wang, 2020] that split atom features along the row dimension, ours are more flexible since the number of coupling layers is decoupled from the number of nodes.

4.3 CC-Glow

CC-Glow fA aims to learn P(A) by transforming bond tensor A into latent variables zA = fA(A) in analogy with multi-channel images. We adopt the architecture of multi-scale Glow but propose a novel Criss-Cross Attention (CCA) [Huang et al., 2019] based coupling layer to offer a better relational inductive bias for modeling bonds.

CCA-based coupling layers. For modeling bonds, it is essential to utilize the local information of the incident atoms. However, the atommatrix is unavailable at this stage in ensuring the reversibility of MolHF. Therefore, we decide to assign atoms with features aggregated from relevant bonds (attached to the atom) through message passing [Gilmer et al., 2017]. Based on the above analysis, we achieve this goal by introducing CCA to extract features from the bond tensor in horizontal and vertical directions as the contextual messages of two atoms connected by a bond, respectively. See Appendix A.1 for more implementation details. We build the coupling network with three layers, where the first and last layers are 3 × 3 convolutions, while the center layer is CCA.

4.4 Two-stepMolecular Graph Generation

Benefiting from the reversibility of MolHF, molecular graph generation is simply executed by drawing z = (zX|A, zA) from PZ(z) and calling the inverse transformation of MolHF. In our two-step generative manner, the bond tensor ˆ A is generated by the inverse transformation of fA at first, then the atom matrix ˆX is generated by the inverse transformation of fX|A, conditioned on the pre-generated ˆ A. Formally, this process can be summarized as,

where the composite operation onehot · argmax maps the dequantized continuous data back to the discrete one-hot data. Lastly, a post-hoc validity correction mechanism [Zang and Wang, 2020] is incorporated to ensure that the generated molecular graphs are valid. This generation process can exploit the hierarchical structural information captured by the inverse transformation of MolHF.

5 Experiments

Following previous works [Jin et al., 2018; Zang and Wang, 2020], we demonstrate the effectiveness of our MolHF in various aspects: generation and reconstruction, property optimization, and hierarchical visualization.

Datasets. We evaluate our method on the ZINC250K dataset [Irwin et al., 2012] for a fair comparison. In addition, we also use the Polymer dataset [St. John et al., 2019] to verify that our method is capable of scaling to larger molecules. The ZINC250K dataset consists of 249,455 molecules with up to 38 atoms of 9 different types, and the Polymer dataset consists of 86,973 larger polymer-like molecules with up to 122 atoms of 7 different types. For the preprocessing strategy, we mostly follow [Zang and Wang, 2020], the difference is that we pad all the molecules in the ZINC250K and Polymer with 2 and 6 extra virtual atoms respectively, to ensure that each molecular graph can be coarsened over multiple scales.

Baselines. We compare MolHF with the following baselines. JT-VAE [Jin et al., 2018] and HierVAE [Jin et al., 2020] are VAE-based models. GCPN [You et al., 2018] is an autoregressive model. GraphAF [Shi et al., 2020] and GraphDF [Luo et al., 2021] are autoregressive flow-based models. MoFlow [Zang and Wang, 2020], GraphCNF [Lippe and Gavves, 2021], and MolGrow [Kuznetsov and Polykovskiy, 2021] are one-shot flowbased models.

Implementation. For the generation and reconstruction, the latent variables are assumed to follow a prior isotropic Gaussian distribution N(0, σ2 I), where σ is a learnable parameter denoting the standard deviation. For property optimization and hierarchical visualization, the same model trained on the ZINC250K dataset is used for all experiments. Further implementation details can be found in Appendix C.

5.1 Generation and Reconstruction

Setup. We evaluate the ability of MolHF to generate and reconstruct molecular graphs via the widely used metrics: Validity and Validity w/o correction are the percentages of chemically valid molecules among all the generated graphs with and without validity correction, respectively. Uniqueness is the percentage of unique molecules among all the generated valid molecules. Novelty is the percentage of generated valid molecules that do not appear in the training set. Reconstruction is the percentage of molecules in the training set that can be reconstructed from latent variables. As for the experiments on the Polymer dataset, to verify that MolHF does generate more large molecules, we introduce Validity w/ filter to calculate the percentage of valid large molecules with more than 38 atoms (the maximum number of atoms in the ZINC250K dataset) among all the generated graphs without validity correction. In practice, higher-quality samples are drawn from the reduced-temperature models [Kingma and Dhariwal, 2018]. To control the trade-off between quality and diversity, we sample from the distribution N(0, (tσ)2 I) with the temperature t < 1, similar to MoFlow and MolGrow. To be comparable with the aforementioned baselines, we sample 10,000 molecular graphs and report the mean and standard deviation of results over 5 runs.

Results. Table 1 shows that MolHF outperforms the stateof-the-art baselines on all the metrics for the ZINC250K dataset. Thanks to the validity correction mechanism, MolHF achieves 100% Validity. We argue that Validity w/o correction can better assess how well models capture the data distribution, since the corrected graphs can be completely different from the original model outputs. Compared with one-shot flow-based models MoFlow and GraphCNF, MolHF attains an increase of 13.1% and 11.5% in Validity w/o correction respectively, implying the superiority of modeling the molecular graphs with hierarchical architecture. Compared with another hierarchical flow-based model MolGrow, MolHF outperforms it by a large margin (37.1%) with respect to Validity w/o correction, indicating that exploiting the coarsened graph structure is crucial for enhancing the performance of hierarchical generative models. Owing to the invertible nature, MolHF is able to reconstruct all molecules, whereas JT-VAE and GCPN can not. Furthermore, experiments on the Polymer dataset demonstrate that MolHF can be applied to model larger molecules that are impractical with previous flow-based models, as shown in Table 2. We find that the performance of all flow-based baselines degrades substantially when modeling larger molecules. MoFlow even cannot fit the data distribution to generate valid molecules. While GraphAF and GraphDF can generate a few small molecules in the early stage of the autoregressive generation process, they become fallible as the current molecules get larger so that only a tiny fraction of the valid molecules are large. GraphCNF, the best flow-based baseline, only generates 11.8% large molecules, which is still unsatisfactory. By contrast, MolHF can generate approximately 43% valid molecules, most of which are large. These results imply that the hierarchical architecture exploiting the coarsened graph structures enables MolHF to model the more complex distribution by better capturing local and global information of large molecules. Compared with the substructure-based model HierVAE, MolHF generates more unique molecules and reconstructs all molecules. We illustrate several valid large molecules (with rich topological complexity) generated by MolHF in Appendix B.1.

Empirical running time. MolHF is significantlymore efficient in aspects of training and sampling. For the ZINC250K dataset, the reported running time of MoFlow and MolGrow is 22 hours and 2 days respectively, while MolHF converges faster and can get the results within 12 hours. For the Polymer dataset, due to the huge memory cost for parallel training, autoregressive GraphAF and GraphDF cost roughly 50 minutes/epoch to train, while MolHF is 3× faster. Besides, GraphAF and GraphDF cost 2 hours and 20 minutes to randomly sample 10,000 molecules (without validity correction) respectively, whereas MolHF only takes 50 seconds, 144× and 24× faster. The efficiency is compared on the same computing facilities using 1 Tesla V100 GPU.

5.2 Property Optimization

Setup. The unconstrained property optimization task aims at discovering novel molecules with the best QED, which measures the druglikeness of the molecules. The constrained property optimization task aims at modifying the given molecules to improve PlogP, which denotes logP score penalized by ring size and synthetic accessibility, with the constraint that the Tanimoto similarity of Morgan fingerprint between the modified and original molecules is above a threshold δ [Jin et al., 2018]. We formulate property optimization as a latent space optimization (LSO) problem [G´omezBombarelli et al., 2018]. Formally, we first train a surrogate network h which takes as input the latent variables z encoded with the trained MolHF to predict the property values, then perform gradient ascent z ← z + α · ∇zh iteratively to navigate the latent space. We decode the latest latent variables with MolHF to obtain the optimized molecule. For constrained property optimization, we start from the latent variables of 800 molecules with the lowest PlogP [Jin et al., 2018]. Note that two different sets of molecules selected from the test set [Jin et al., 2018] or the entire set [Shi et al., 2020] are used in prior methods. We thus report results separately.

Results. We show the best scores discovered in the unconstrained optimization task in Table 3. Like MoFlow, MolHF finds more novel molecules with the best QED score (0.948) appearing in the dataset. Further, we illustrate the histogram of the optimized property scores in Figure 2 to verify that MolHF finds molecules with higher QED scores than MoFlow on average (0.70 vs. 0.56). As shown in Table 4 and Table 5, MolHF achieves higher property improvement than GraphDF under 3 of 4 similarity constraints and outperforms GraphAF and MolGrow under 2 of 4 similarity constraints, while maintaining higher similarities. Compared with the autoregressive GCPN, GraphDF, and GraphAF which finetune the generator with reinforcement learning and optimize each molecule hundreds of times, the optimization strategy of MolHF is more conveniently implementable and efficient. This nonetheless results in the issues of success rate and robustness. Fortunately, these issues can be alleviated with the recent advanced LSO approaches [Tripp et al., 2020; Notin et al., 2021], potentially leading to better performance.

5.3 Hierarchical Visualization

To investigate the generative mechanism of MolHF, we visualize the generated molecules in terms of the semantics of the learned latent variables and the generated substructures.

Semantics of the learned latent variables. We infer the latent variables and resample the latent variables of the finest scale, increasing the coarsest scale affected by this resampling [Dinh et al., 2016]. We decode the resampled latent variables and illustrate several cases in Figure 3, where the leftmost column represents the original molecules, and the subsequent columns (from left to right) are obtained by resampling the latent variables of coarser scale, resampling more and more as we go right. As we can see, latent variables of finer scale affect the local structure (atom and bond), while those of coarser scale affect the global structure (skeletal structure).

Generated substructures. In Figure 4, we illustrate some of the most frequent substructures generated by MolHF. These learned substructures are chemically meaningful, demonstrating that MolHF does capture the rich valid structure and semantic information of molecular graphs. By comparison, we find that these substructures (such as 2butylene, methylethylamine, and acetamide, to name a few) often appear in the fragment vocabulary [Jin et al., 2018; Jin et al., 2020; Xie et al., 2021; Kong et al., 2022] manually extracted based on frequency of occurrence. Thus, we argue that MolHF successfully extracts common substructures from the data without any need for domain knowledge, and matches the prior distribution of fragments to some extent. Meanwhile, due to the hierarchical generation, some of the substructures generated by MolHF are composed of smaller structural motifs. For instance, the substituted benzene rings can be further divided according to the retrosynthesis-based BRICS algorithm [Degen et al., 2008]. To link coarsened graph structure with substructures, we decompose several molecules with the same coarsened graph structure in Appendix B.2, where atoms are colored according to the substructures to which they belong. We observe that MolHF learns to hierarchically generate molecular graphs in a coarseto-fine manner by first generating a coarsened graph structure (skeletal structure) and then refining the substructures.

6 Discussion and Conclusion

In this paper, we present MolHF, a hierarchical flow-based model to generate molecular graphs in a coarse-to-fine manner. As demonstrated with extensive experiments, MolHF achieves state-of-the-art performance in random generation and property optimization, while being capable of modeling larger molecules that are impractical with prior flow-based methods. We hope that our MolHF can spur further research to unleash the great potential of hierarchical generative models in the graph generation field.

Limitations and future work. The strategy used in the property optimization task is inapplicable for real-world drug discovery, as it is expensive to evaluate the properties of molecules with wet-lab experiments [Zhu et al., 2023]. We believe that developing novel sample-efficient optimization approaches is a fruitful avenue. In the future, one might extend MolHF to generate other types of graphs with hierarchical structures.

Acknowledgments

This research was partially supported by the National Key R&D Program of China under grant No.2019YFC0118802, Key R&D Program of Zhejiang Province under grant No.2020C03010, National Natural Science Foundation of China under grant No.62176231. Tingjun Hou's research was supported in part by the National Natural Science Foundation of China under grant No.22220102001.

Contribution Statement

Yiheng Zhu and Zhenqiu Ouyang contributed equally to this work.