Year: 2021

Abstract:
Molecule generation, which requires generating valid molecules with desired properties, is a fundamental but challenging task. Recent years have witnessed the rapid development of atom-level auto-regressive models, which usually construct graphs following sequential actions of adding atom-level nodes and edges. However, these atom-level models ignore high-frequency substructures, which not only capture the regularities of atomic combination in molecules but are also often related to desired chemical properties, and therefore may be sub-optimal for generating high-quality molecules. In this paper, we propose a method to automatically discover such common substructures, which we call graph pieces, from given molecular graphs. We also present a graph piece variational autoencoder (GPVAE) for generating molecular graphs based on graph pieces. Experiments show that our GP-VAE models not only achieve better performance than the state-ofthe-art baseline for distribution-learning, property optimization, and constrained property optimization tasks but are also computationally efficient.

text: Molecule generation, which requires generating valid molecules with desired properties, is a fundamental but challenging task. Recent years have witnessed the rapid development of atom-level auto-regressive models, which usually construct graphs following sequential actions of adding atom-level nodes and edges. However, these atom-level models ignore high-frequency substructures, which not only capture the regularities of atomic combination in molecules but are also often related to desired chemical properties, and therefore may be sub-optimal for generating high-quality molecules. In this paper, we propose a method to automatically discover such common substructures, which we call graph pieces, from given molecular graphs. We also present a graph piece variational autoencoder (GPVAE) for generating molecular graphs based on graph pieces. Experiments show that our GP-VAE models not only achieve better performance than the state-ofthe-art baseline for distribution-learning, property optimization, and constrained property optimization tasks but are also computationally efficient.

heading: 1 INTRODUCTION
text: Molecule generation is a task that aims to produce chemically valid molecules with optimized properties. It is important for a variety of applications, such as drug discovery and material science. Graph-based molecule generation models, which are robust to molecule substructures (You et al., 2018; Kwon et al., 2019), have gained increasing attention recently (Jin et al., 2018; Li et al., 2018a; You et al., 2018; Kwon et al., 2019; De Cao & Kipf, 2018; Shi et al., 2020; Jin et al., 2020b).
Graph-based molecule generation models typically decompose molecular graphs into sequential actions of generating atoms and bonds autoregressively (Li et al., 2018b; You et al., 2018; Li et al., 2018a; Jin et al., 2020b). While this decomposition is natural and straightforward, it inevitably ignores the existence of common substructures in molecular graphs, as illustrated in Figure 1. Compared with using atoms for generating molecules, using graph substructures for generating molecules have three potential benefits. First, using substructures can capture the regularities of atomic combination in molecules, and therefore is more capable of generating realistic molecules. Second, using substructures can better capture chemical properties, as there is a correlation between substructures and chemical properties (Murray & Rees, 2009; Jin et al., 2020b). Third, using substructures enables efficient training and inference. It is evident that using substructures to represent molecular graphs can result in much shorter sequences, therefore the training and inference process can be accelerated. As a result, we believe that models using substructures to represent molecular graphs can generate more realistic molecules with better-optimized properties efficiently.

heading: 83.9% 60.5% 41.3% 32.6% 27.5%
text: In this paper, we present an iterative algorithm to automatically discover common substructures in molecules, which we call graph pieces. Initially, graph pieces correspond to single atoms that appear in graphs of a given dataset. Then for each iteration, we count the occurrence of neighboring pieces in graphs and merge the most frequent neighboring pieces into a new graph piece. Since substructures can be seen as small molecules, we use SMILES (Weininger, 1988), a text-based representation for molecules, to efficiently judge whether two graph pieces are identical. To effectively utilize these substructures, we also propose a graph piece variational autoencoder (GP-VAE). Our model consists of a graph neural network (GNN, Scarselli et al., 2008) encoder and a two-step decoder. The two-step decoder first generates graph pieces auto-regressively and then predicts atom-level bonds between graph pieces in parallel. As a result, our GP-VAE decouples the alternated generation of nodes and edges, achieving significant computational efficiency for both training and generation.
We conduct extensive experiments on ZINC250K (Irwin et al., 2012) and QM9 (Blum & Reymond, 2009; Rupp et al., 2012) datasets. Results demonstrate that our GP-VAE models outperform state-of-the-art models on distribution-learning, property optimization, and constrained property optimization tasks, and are about six times faster than the fastest baseline.

heading: 2 RELATED WORK
text: Molecule Generation Based on different representations for molecules, molecule generation models can be divided into two categories: text-based and graph-based. Text-based models (GómezBombarelli et al., 2018; Kusner et al., 2017; Bjerrum & Threlfall, 2017), which usually adopt the Simplified Molecular-Input Line-entry System (SMILES) (Weininger, 1988) representation, are simple and efficient methods for generating molecules. However, these models are not robust because a single perturbation in the text molecule representation can result in significant changes in molecule structure (You et al., 2018; Kwon et al., 2019). Graph-based models (De Cao & Kipf, 2018; Shi et al., 2020; Jin et al., 2020b), therefore, have gained increasing attention recently. Li et al. (2018b) proposed a generation model of graphs and demonstrated it performed better than text-based generation models on molecule generation. You et al. (2018) used reinforcement learning to fuse rewards of chemical validity and property scores into each step of generating a molecule. Popova et al. (2019) proposed an MRNN to autoregressively generate nodes and edges based on the generated graph. Shi et al. (2020) proposed a flow-based autoregressive model and use reinforcement learning for the goal-directed molecular graph generation. However, these models use atomlevel graph representation, which results in very long sequences, and therefore the training process is typically time-consuming. Our method is graph-based and uses substructure-level representation for graphs, which not only captures chemical properties but also is computationally efficient.
Substructure-level Graph Representation Jin et al. (2018) proposed to generate molecules in the form of junction trees where each node is a ring or edge. Jin et al. (2020a) decomposed molecules into substructures by breaking all the bridge bonds. It used a complex hierarchical model for polymer generation and graph-to-graph translation. Jin et al. (2020b) proposed to extract the smallest substructure which maintains the original chemical property. The extracted substructures usually include most atoms of the original molecules, which are too coarse-grained and exert limitations on the search space. The models proposed by Jin et al. (2020a) and Jin et al. (2020b) are not suitable for the experiments in this paper since they either need graph-to-graph supervised data or are incompatible with continuous properties. There exists various methods to discover frequent subgraphs (Inokuchi et al., 2000; Yan & Han, 2002; Nijssen & Kok, 2004). However, they have difficulty decomposing graphs into frequent subgraphs (Jazayeri & Yang, 2021) since they mainly aim to discover frequent subgraphs as additional features for downstream network analysis. Therefore they can hardly be applied to substructure-level molecular graph representation. Different from Jin et al. (2018; 2020a) which use manual rules to extract substructures, we automatically extract common substructures which better capture the regularities in molecules for substructure-level decomposition.

heading: 3 APPROACH
text: We first give the definition of graph pieces and algorithms for graph piece extraction in Section 3.1. Then we describe the encoder of our graph piece variational autoencoder (GP-VAE) model in Section 3.2. Finally, we descibe the two-step decoder of GP-VAE in Section 3.3.
Algorithm 1: Graph Piece Extraction Input: A set of graphs D and the desired number N of graph pieces to learn. Result: A set of graph pieces S and the counter F of graph pieces.
1 begin 2 S ← {GraphToSMILES(〈{a}, ∅〉)}; . Initially, S corresponds to all atoms a that appear in D. 3 N ′ ← max(N, |S|); 4 while |S| < N ′ do 5 F ← EmptyMap(); . Initialize a counter. 6 foreach G in D do 7 forall 〈Pi,Pj , Ẽij〉 in G do 8 P ← Merge(〈Pi,Pj , Ẽij〉); . Merge neighboring graph pieces into a new graph piece. 9 s← GraphToSMILES(P); . Convert a graph to SMILES representation.
10 F [s] = F [s] + 1; . Update the counter, the default value for a new s is 0. 11 end 12 end 13 s = TopElem(counter); . Find the most frequent merged graph piece. 14 P ← SMILESToGraph(s); . Convert the SMILES string to graph representation. 15 S ← S ∪ {s}; D′ ← {}; 16 foreach G in D do 17 G′ ← MergeSubGraph(G,P); . Update the graph representation if possible. 18 D′ ← D′ ∪ {G′}; 19 end 20 D ← D′ 21 end 22 end
3.1 GRAPH PIECE
A molecule can be represented as a graph G = 〈V, E〉, where V is a set of nodes that correspond to atoms and E is a set of edges that correspond to chemical bonds. Instead of using atoms, we use substructures, which we call graph pieces, as building blocks. We define a graph pieceP as a subgraph 〈Ṽ, Ẽ〉 that appears in a graph G, where Ṽ ⊆ V and Ẽ ⊆ E . It should be noted that either a single atom or a whole graph is a valid graph piece. Given a set of graph pieces S, suppose the graph G can be decomposed into n graph pieces in S, then G can be alternatively represented as 〈{Pi}, {Ẽij}〉, where Ẽij denotes the set of edges between two neighboring graph pieces Pi and Pj . The decomposition of a graph G into graph pieces satisfies the following constraints: (1) the union of all atoms in the graph pieces equals to all atoms in the molecule, namely ⋃
i Vi = V; (2) there is no intersection between any two graph pieces, namely ∀i 6= j, Vi ∩Vj = ∅, and Ẽi ∩ Ẽj = ∅; (3) the union of all connections within and between graph pieces equals to all bonds in the molecule, namely ⋃ i,j(Ẽi ∪ Ẽij) = E , where i range from 1 to n and j range from i+ 1 to n. Figure 2 shows an decomposed molecule.
The algorithm for extracting graph pieces from a given set of graphs D is given in Algorithm 1. Our algorithm draws inspiration from Byte Pair Encoding (Gage, 1994, BPE). Initially, a graph G inD is decomposed into atom-level graph pieces and the vocabulary S of graph pieces is composed of all unique atom-level graph pieces that appear in D. Given the number N of graph pieces to learn, at each iteration, our algorithm enumerates all neighboring graph pieces and edges that connect the two graph pieces in G, namely 〈Pi,Pj , Ẽij〉. As 〈Pi,Pj , Ẽij〉 is also a valid subgraph, we merge it into a graph piece and count its occurrence. We find the most frequent merged graph pieceP and add it into the vocabulary S. After that, we also update graphs G in D that contain P by merging 〈Pi,Pj , Ẽij〉 into P . The algorithm terminates when the vocabulary size exceeds the predefined number N . Note that we use SMILES (Weininger, 1988) to represent a graph piece in our algorithm1, therefore we ensure the uniqueness of a graph piece. A running example of our graph piece extraction algorithm is illustrated in Figure 3. At test time, we first decompose a molecular graph into atom-level graph pieces, then apply the learned operations to merge the graph pieces into larger ones. This process ensures there is a piece-level decomposition for an arbitrary molecule. We provide the pseudo code for the piece-level decomposition in Appendix A for better understanding. We provide the complexity analysis for both algorithms in Appendix B.
1We use RDKit (www.rdkit.org) to perform the conversion between molecular graph and SMILES.

heading: 3.2 GRAPH ENCODER
text: As shown in Figure 4, we use a GNN to encode a molecular graph G represented by graph pieces into a latent variable z. Each node v has a feature vector xv which indicates its atomic type, the type, and generation order of the graph piece it is in. Each edge has a feature vector indicating its bond type. Modern GNNs follow a neighborhood aggregation strategy (Xu et al., 2018), which iteratively update the representations of nodes with AGGREGATE and COMBINE operations. The representation of a node v in the k-th iteration is calculated as follows:
a(k)v = AGGREGATE (k)({h(k−1)u : u ∈ N (v)}), (1)
h(k)v = COMBINE (k)(h(k−1)v ,a (k) v ), (2)
where AGGREGATE and COMBINE vary in different GNNs. N (v) denotes neighboring nodes of v. We set h(0)v = xv and use GIN with edges feature (Hu et al., 2019) as the backbone GNN network, which implements AGGREGATE and COMBINE as follows:
a(k)v = ∑
u∈N (v)
ReLU(h(k−1)u + euv), (3)
h(k)v = hΘ((1 + ε)h (k−1) v + a (k) v ), (4)
where hΘ is a neural network and ε is a constant. We implement hΘ as a 2-layer multilayer perceptron (Gardner & Dorling, 1998, MLP) with ReLU activation and set ε = 0. Since the k-th representation captures k-hop feature of nodes (Xu et al., 2018), we obtain the final representations of nodes as hv = [h (1) v , . . . ,h (t) v ] so that it contains 1-hop to t-hop contextual information. Then
we compute the representation of the graph G through summation hG = ∑ v∈V hv . We use hG to obtain the mean µG and log variance σG of variational posterior approximation q(z|G) through two separate linear layers and use the reparameterization trick (Kingma & Welling, 2013) to sample from the distribution in the training process.

heading: 3.3 TWO-STEP DECODER
text: With a molecule G represented as 〈{Pi}, {Ẽij}〉, our model generates {Pi} and {Ẽij} in two consecutive phases. The two phases move from coarse to fine granularity.
Piece-level Sequence Generation Given a latent variable z, our model first uses an autoregressive sequence generation model P (Pi|P<i, z) to decode a sequence of graph pieces [P1, . . . ,Pn]2. During training, we insert two special tokens “<bos>” and “<eos>” into the begin and the end of the graph piece sequence to indicate the begin and the end, respectively. During testing, the generation stops when a “<eos>” is generated. The sequence model can be RNNs, such as LSTM (Hochreiter & Schmidhuber, 1997) and GRU (Cho et al., 2014). In this work, we use a single layer of GRU and project the latent variable z to the initial state of GRU. The training objective of this stage is to minimize the log-likelihood LP of the ground truth graph piece sequence:
LP = n∑ i=1 − logP (Pi|P<i, z). (5)
2The ordering of graph pieces is not unique. We train with one order and leave this problem for future work.
Bond Completion We generate {Ẽij} non-autoregressively. For the architecture, we use a GNN with the same structure as the graph encoder in Section 3.2 but with different parameters to obtain the representations hv of each atom v in the graph pieces. Given nodes v and u, we predict their connections as follows:
P (euv|z) = Hθ([hv;hu; z]), (6)
where Hθ is a neural network. In this work we adopt a 3-layer MLP with ReLU activation. Apart from the types of chemical bonds, we also add a special type “<none>” to the edge vocabulary which indicates there is no connection between two nodes. During training, we predict both P (euv) and P (evu) to let Hθ learn the undirected nature of chemical bonds. We use negative sampling (Goldberg & Levy, 2014) to balance the ratio of “<none>” and chemical bonds. Since only about 2% pairs of nodes has inter-piece connections, negative sampling significantly improves the computational efficiency and scalability. The training objective of this stage is to minimize the loglikelihood LẼ of the ground truth inter-piece connections:
LẼ = ∑
u∈Pi,v∈Pj ,i6=j
− logP (euv|z). (7)
The reconstruction loss of our GP-VAE model is:
Lrec = LS + LẼ . (8)
To decode {Ẽij} in inference time, the decoder first assigns all possible inter-piece connections with a bond type and a corresponding confidence level. We try to add bonds that have a confidence level higher than δth = 0.5 to the molecule in order of confidence level from high to low. For each attempt, we perform a valency check and a cycle check to reject the connections that will cause violation of valency or form unstable rings which are too small or too large. Since this procedure may form unconnected graphs, we find the maximal connected component as the result of the generation. We present the pseudo code for inference in Appendix C.
We jointly train a 2-layer MLP from z to predict the scores of target properties using the MSE loss. Denote the loss of the predictor as Lprop, the loss function of our GP-VAE is:
L = αLrec + (1− α)Lprop + βDKL, (9)
where α and β balance the reconstruction loss, the property loss, and the KL divergence DKL between the distribution of z and the prior distribution N (0, I).

heading: 4 EXPERIMENTS

heading: 4.1 SETUP
text: Evaluation Tasks We first report the empirical results for the distribution-learning tasks of the GuacaMol benchmark (Brown et al., 2019) to evaluate the ability of the model to generate realistic and diverse molecules. Then we validate our model on two goal-directed tasks. Property Optimization requires generating molecules with optimized properties. Constrained Property Optimization concentrates on improving the properties of given molecules with a restricted degree of modification.
Dataset We use the ZINC250K (Irwin et al., 2012) dataset for training, which contains 250,000 drug-like molecules up to 38 atoms. For GuacaMol benchmark, we also add results on the QM9 (Blum & Reymond, 2009; Rupp et al., 2012) dataset, which has 7,165 molecules up to 23 atoms.
Baselines We compare our graph piece variational autoencoder (GP-VAE) with the following state-of-the-art models. JT-VAE (Jin et al., 2018) is a variational autoencoder that represents molecules as junction trees. It performs Bayesian optimization on the latent variable for property optimization. GCPN (You et al., 2018) combines reinforcement learning and graph representation for goal-directed molecular graph generation. MRNN (Popova et al., 2019) adopts two RNNs to autoregressively generate atoms and bonds respectively. It combines policy gradient optimization to generate molecules with desired properties. GraphAF (Shi et al., 2020) is a flow-based autoregressive model which is first pretrained for likelihood modeling and then fine-tuned with reinforcement learning for property optimization.GA (Nigam et al., 2020) adopts genetic algorithms for property optimization and model the selection of the subsequent population with a neural network.
Implementation Details We choose N = 300 for property optimization and N = 500 for constrained property optimization. GP-VAE is trained for 6 epochs with a batch size of 32 and a learning rate of 0.001. We set α = 0.1 and initialize β = 0. We adopt a warm-up method that increases β by 0.002 every 1000 steps to a maximum of 0.01. More details can be found in Appendix D.

heading: 4.2 RESULTS
text: GuacaMol Distribution-Learning Benchmarks The distribution-learning benchmarks incorporate five metrics on 10,000 molecules generated by the models. Validity measures whether the generated molecules are chemically valid. Uniqueness penalizes models when they generate the same molecule multiple times. Novelty assesses the ability of the models to generate molecules that are not present in the training set. KL Divergence measures the closeness of the probability distributions of a variety of physicochemical descriptors for the training set and the generated molecules. Fréchet ChemNet Distance (FCD) calculates the closeness of the two sets of molecules with respect to their hidden representations in the ChemNet trained for predicting biological activities. Each metric is normalized to 0 to 1, and a higher value indicates better performance. Table 1 shows the results of distribution-learning benchmarks on QM9 and ZINC250K. Our model achieves competitive results in all five metrics, which indicates our model can generate realistic molecules and does not overfit the training set. We present some molecules sampled from the prior distribution in Appendix I.
Property Optimization This task focuses on generating molecules with optimized Penalized logP (Kusner et al., 2017) and QED (Bickerton et al., 2012). Penalized logP is logP penalized by synthesis accessibility and ring size which has an unbounded range. QED measures the drug-likeness of molecules with a range of [0, 1]. Both properties are calculated by empirical prediction models (Wildman & Crippen, 1999; Bickerton et al., 2012), and are widely used in previous works (Jin et al., 2018; You et al., 2018; Shi et al., 2020). We adopt the scripts of Shi et al. (2020) to calculate property scores so that the results are comparable. We directly perform gradient ascending on the latent variable (Jin et al., 2018; Luo et al., 2018). Hyperparameters can be found in Appendix D. Following previous works (Jin et al., 2018; You et al., 2018; Shi et al., 2020), we generate 10,000 optimized molecules from the latent space and report the top-3 scores found by each model. Results in Table 2 show that our model surpasses the state-of-the-art models consistently.
Constrained Property Optimization This task concentrates on improving the property scores of given molecules with the constraint that the similarity between the original molecule and the modified molecule is above a threshold δ. Following Jin et al. (2018); You et al. (2018); Shi et al. (2020), we choose 800 molecules with lowest Penalized logP in the test set of ZINC250K for optimization and Tanimoto similarity with Morgan fingerprint (Rogers & Hahn, 2010) as the similarity metric.
Similar to the property optimization task, we perform gradient ascending on the latent variable with a max step of 80. We collect all latent variables which have better-predicted scores than the previous iteration and decode each of them 5 times, namely up to 400 molecules. Following Shi et al. (2020), we initialize the generation with sub-graphs sampled from the original molecules. Then we choose the one with the highest property score from the molecules that meet the similarity constraint. Table 3 shows our model can generate molecules with a higher Penalized logP score while satisfying the similarity constraint. Since our model uses graph pieces as building blocks, the degree of modification tends to be greater than atom-level models, leading to a lower success rate. Nevertheless, our model still manages to achieve high success rates close to atom-level models.

heading: 4.3 RUNTIME COST
text: We train JT-VAE, GraphAF and our GP-VAE on a machine with 1 NVIDIA GeForce RTX 2080Ti GPU and 32 CPU cores and use them to generate 10,000 molecules to compare their efficiency of training and inference. As shown in Table 4, our model achieves significant improvements on the efficiency of training and inference due to graph pieces and the two-step generation. With graph pieces as building blocks, the number of steps required to generate a molecules is significantly decreased compared to the atom-level models like GraphAF. Moreover, since the two-step generation approach separates the generation of graph pieces and the connections between them into two stage and formalizes the bond completion stage as a link prediction task, it avoids the complex enumeration of possible combinations adopted by JT-VAE. Therefore, our model achieves tremendous improvement in computational efficiency over both atom-level and substructure-level baselines.

heading: 4.4 ABLATION STUDY
text: We conduct an ablation study to further validate the effects of graph pieces and the two-step generation approach. We first downgrade the vocabulary to contain only single atoms. Then we replace the two-step decoder with a fully auto-regressive decoder. We present the performance on property optimization and constrained property optimization after the modification in Table 5 and Table 6. The direct introduction of the two-step generation approach leads to improvement on property optimization but harms the performance on constrained property optimization. This is because separating the generation of atoms and bonds brings massive loss of bond information for the atom generation process. However, the adoption of graph pieces as building blocks alleviates this negative effect since the graph pieces themselves contain abundant bond information. Therefore, while the twostep generation approach enhances computational efficiency, the state-of-the-art performance of our model is mainly credited to the use of graph pieces.

heading: 5 ANALYSIS

heading: 5.1 GRAPH PIECE STATISTICS
text: We compare the statistical characteristics of the vocabulary of JT-VAE that contains 780 substructures and the vocabulary of graph pieces with a size of 100, 300, 500, and 700. Figure 5 shows the proportion of substructures with different numbers of atoms in the vocabulary and their frequencies of occurrence in the ZINC250K dataset. The substructures in the vocabulary of JT-VAE mainly concentrate on 5 to 8 atoms with a sharp distribution. However, starting from substructures with 3 atoms, the frequency of occurrence is already close to zero. Therefore, the majority of substructures in the vocabulary of JT-VAE are actually not common substructures. On the contrary, the substructures in the vocabulary of graph pieces have a relatively smooth distribution over 4 to 10 atoms. Moreover, these substructures also have a much higher frequency of occurrence compared to those in the vocabulary of JT-VAE. We present samples of graph pieces in Appendix H.

heading: 5.2 PROPER GRANULARITY
text: A largerN in the graph piece extraction process leads to an increase in the number of atoms in graph pieces and a decrease in their frequency of occurrence, as illustrated in Figure 6. These two factors affect model performance in opposite ways. On the one hand, the entropy of the dataset decreases with more coarse-grained decomposition (Martin & England, 2011), which benefits model learning (Bentz & Alikaniotis, 2016). On the other hand, the sparsity problem worsens as the frequency of graph pieces decreases, which hurts model learning (Allison et al., 2006). We propose a quantified method to balance entropy and sparsity. The entropy of the dataset given a set of graph pieces S is defined by the sum of the entropy of each graph piece normalized by the average number of atoms:
HS = − 1
nS ∑ P∈S P (P) logP (P), (10)
where P (P) is the relative frequency of graph piece P in the dataset and nS is the average number of atoms of graph pieces in S. The sparsity of S is defined as the reciprocal of the average frequency of graph pieces normalized by the size of the dataset:
SS = M
fS , (11)
where M is the number of molecules in the dataset and fS is the average frequency of occurrence of graph pieces in the dataset. Then the entropy - sparsity trade-off (T ) can be expressed as:
TS = HS + γSS (12)
where γ balances the impacts of entropy and sparsity since the impacts vary across different tasks. We assume that TS negatively correlates with downstream tasks. Given a task, we first sample several values of N to calculate their values of T and then compute the γ that minimize the Pearson correlation coefficient between T and the corresponding performance on the task. With the proper γ, Pearson correlation coefficients for the three downstream tasks in this paper are -0.987, -0.999, and
-0.707, indicating strong negative correlations. For example, Figure 6 shows the curve of entropy - sparsity trade-off with a maximum of 3,000 iteration steps for the property optimization task. From the curve, we choose N = 300 for the property optimization task.

heading: 5.3 GRAPH PIECE - PROPERTY CORRELATION
text: To analyze the graph piece-property correlation and whether our model can discover and utilize the correlation, we present the normalized distribution of generated graph pieces and Pearson correlation coefficient between the graph pieces and the score of Penalized logP (PlogP) in Figure 7.
The curve of the Pearson correlation coeffcient indicates that some graph pieces positively correlate with PlogP and some negatively correlate with it. Compared with the flat distribution under the nonoptimization setting, the generated distribution shifts towards the graph pieces positively correlated with PlogP under the PlogP-optimization setting. The generation of graph pieces negatively correlated with PlogP is also suppressed. Therefore, correlations exist between graph pieces and PlogP, and our model can accurately discover and utilize these correlations for PlogP optimization.

heading: 6 CONCLUSION
text: We propose an algorithm to automatically discover the regularity in molecules and extract them as graph pieces. We also propose a graph piece variational autoencoder utilizing the graph pieces and generate molecules in two phases. Our model consistently outperforms state-of-the-art models on distribution-learning, property optimization and constrained property optimization with higher computational efficiency.Our work provides insights into the selection of granularity on molecular graph generation and can inspire future search in this direction.

heading: A PIECE-LEVEL DECOMPOSITION ALGORITHM
text: Algorithm 2 presents the pseudo code for the piece-level decomposition of molecules. The algorithm takes the atom-level molecular graph, the vocabulary of graph pieces and their frequencies of occurrence recorded during the graph piece extraction process as input. Then the algorithm iteratively merge the graph piece pair which has the highest recorded frequency of occurrence in the vocabulary until all graph piece pairs are not in the vocabulary.
Algorithm 2: Piece-Level Decomposition Input: A graph G that decomposed into atom-level pieces, the set S of learned graphs pieces, and the counter F of learned graph pieces. Result: A new representation G′ of G that consists of graph pieces in S.
1 begin 2 G′ ← G; 3 while True do 4 freq ← −1; P ← None; 5 forall 〈Pi,Pj , Ẽij〉 in G′ do 6 P ′ ← Merge(〈Pi,Pj , Ẽij〉); . Merge neighboring graph pieces into a new graph piece. 7 s← GraphToSMILES(P ′); . Convert a graph to SMILES representation. 8 if s in V and F [s] > freq then 9 freq ← F [s];
10 P ← P ′; 11 end 12 end 13 if freq == −1 then 14 break; 15 else 16 G′ ← MergeSubGraph(G′,P); . Update the graph representation. 17 end 18 end 19 end

heading: B COMPLEXITY ANALYSIS
text: Graph Piece Extraction Since the number of graph piece pairs equals the number of inter-piece connections in the piece-level graph, the complexity is O(NMe), where N is the predefined size of vocabulary, M denotes the number of molecules in the dataset, and e denotes the maximal number of inter-piece connections in a single molecule. The number of inter-piece connections decreases rapidly in the first few iterations, therefore the time cost for each iteration decreases rapidly. It cost 6 hours to perform 500 iterations on 250,000 molecules in the ZINC250K dataset with 4 CPU cores.
Piece-Level Decomposition Given an arbitrary molecule, the worst case is that each iteration adds one atom to one existing subgraph until the molecule is finally merged into a single graph piece. In this case the algorithm runs for |V| iterations. Therefore, the complexity is O(|V|) where V includes all the atoms in the molecule.
C INFERENCE ALGORITHM FOR BOND COMPLETION
Algorithm 3 shows the pseudo code of our inference algorithm. We first predict the bonds between all possible pairs of atoms in which the two atoms are in different graph pieces and sort them by the confidence level given by the model from high to low. Then for each bond with a confidence level higher then the predefined threshold δth, which is 0.5 in our experiments, we add it into the molecular graph if it passes the valence check and cycle check. The valence check ensures the given bond will not violate valence rules. The cycle check ensures the given bond will not form unstable rings with nodes less than 5 or more than 6.
Algorithm 3: Inference Algorithm for Bond Completion Input: An incomplete molecular graph G composed of graph pieces where inter-piece bonds
are absent, the predicted bond type for all possible inter-piece connections B and the map to their confidence level C, the threshold for confidence level δth
Result: A valid molecular graph G′ 1 begin 2 G′ ← G; 3 B ← SortByConfidence(B, C); . Sort the bonds by their confidence level from high to low. 4 forall buv in B do 5 if C[buv] < δth then 6 continue; . Discard edges with confidence level lower than the threshold. 7 end 8 if valence check(buv) and cycle check(buv) then 9 G′ ← AddEdge(G′, buv); . Add edges that pass valence and cycle check to G′
10 end 11 end 12 G′ ← MaxConnectedComponent(G′); . Find the maximal connected component in G′ 13 end

heading: D EXPERIMENT DETAILS
text: Model and Training Hyperparameters We present the choice of model parameters in Table 7 and training parameters in Table 8. We represent an atom with three features: atom embedding, piece embedding and position embedding. Atom embedding is a trainable vector of size eatom for each type of atoms. Similarly, piece embedding is a trainable vector of size epiece for each type of pieces. Positions indicate the order of generation of pieces. We jointly train a 2-layer MLP from the latent variable to predict property scores. The training loss is represented as L = α · Lrec + (1 − α) · Lprop + β · DKL where α balances the reconstruction loss and prediction loss. For β, we adopt a warm-up method that increase it by βstage every fixed number of steps to a maximum of βmax. We found a β higher than 0.01 often causes KL vanishing problem and greatly harm the performance. Our model and the baselines are trained on the ZINC250K dataset with the same train / valid / test split as in Kusner et al. (2017).
Property Optimization We use gradient ascending to search in the continuous space of latent variable. For simplicity, we set a target score and optimize the mean square error between the score given by the predictor and the target score just as in the training process. The optimization stops if the mean square error does not drop for 3 iterations or it has been iterated to the maxstep. We normalize the Penalized logP in the training set to [0, 1] according to the statistics of ZINC250K. By setting a target value higher than 1 the model is supposed to find molecules with better property than the molecules in the training set. To acquire the best performance, we perform a grid search with lr ∈ {0.001, 0.01, 0.1, 1, 2}, maxstep ∈ {20, 40, 60, 80, 100} and target ∈ {1, 2, 3, 4}. For optimization of QED, we choose lr = 0.01,maxstep = 100, target = 2. For optimization of Penalized logP, we choose lr = 0.1,maxstep = 100, target = 2.
Constrained Property Optimization We use the same method as property optimization to optimize the latent variable. We also perform a grid search with lr ∈ {0.1, 0.01} and target ∈ {2, 3}. We select lr = 0.1,maxstep = 80 and target = 2. For decoding, we first initialize the generation with a submol sampled from the original molecule by teacher forcing. We follow Shi et al. (2020) to first sample a BFS order of all atoms and then randomly drop out the last m atoms with m up to 5. We collect all latent variables which have better predicted scores than the previous iteration and decode each of them 5 times, namely up to 400 molecules. Then we choose the one with the highest property score from the molecules that meet the similarity constraint. For the baseline GA (Ahn et al., 2020), we adjust the number of iterations to 5 and the size of population to 80, namely traversing up to 400 molecules, for fair comparison.

heading: E DATA EFFICIENCY
text: Since the graph pieces are common subgraphs in the molecular graphs, they should be relatively stable with respect to the scale of training set. To validate this assumption, we choose subsets of different ratios to the training set for training to observe the trend of the coverage of Top 100 graph pieces in the vocabularies as well as the model performance on the average score of the distributionlearning benchmarks. As illustrated in Figure 9, with a subset above 20% of the training set, the constructed vocabulary covers more than 95% of the top 100 graph pieces in the full training set, as well as the model performance on the distribution-learning benchmarks.

heading: F FUSED RINGS GENERATION
text: We conduct an additional experiment to validate the ability of GP-VAE to generate molecules with fused rings (cycles with shared edges), because at first thought it seems difficult for GP-VAE to handle these molecules due to the non-overlapping nature of graph pieces. We train atom-level and piece-level GP-VAEs on all 4,431 structures consisting of fused rings from ZINC250K. Then we sample 1,000 molecules from the latent space to calculate the proportion of molecules with fused rings. The results are 94.5% and 97.2% for the atom-level model and the piece-level model, respectively. The experiment demonstrates that the introduction of graph pieces as building blocks will not hinder the generation of molecules with fused rings.

heading: G DISCUSSION
text: Universal Granularity Adaption The concept and extraction algorithm of graph pieces resemble those of subword units (Sennrich et al., 2015) in machine translation. Though subword units are designed for the out-of-vocabulary problem of machine translation, they also improve the translation quality (Sennrich et al., 2015). In this work, we demonstrate the power of graph pieces and are curious about whether there is a universal way to adapt atom-level models into piece-level counterparts to improve their generation quality. The key challenge is to find an efficient and expressive way to encode inter-piece connections into feature vectors. We leave this for future work.
Searching in Continuous Space In recent years, reinforcement learning (RL) is becoming dominant in the field of optimization of molecular properties (You et al., 2018; Shi et al., 2020). These RL models usually suffer from reward sparsity when applied to multi-objective optimization (Jin et al., 2020b). However, most scenarios that incorporate molecular property optimization have multi-objective constraints (e.g.,drug discovery). In this work, we show that with graph pieces, even simple searching method like gradient ascending can surpass RL methods on single-objective optimization. It is possible that with better searching methods in continuous space our model can achieve competitive results on multi-objective optimization.

heading: H GRAPH PIECE SAMPLES
text: We present 50 graph pieces found by our extraction algorithm in Figure 11.

heading: I MORE MOLECULE SAMPLES
text: We further present 50 molecules sampled from the prior distribution in Figure 12.

Review summaries:
This paper proposes a new molecular graph generation method and empirically shows its effectiveness. The proposed method first decomposes molecular graphs into smaller parts (called graph pieces), followed by training a variational autoencoder so that it can generate the collected graph pieces. The trained model is expected to generate a variety of graphs with desirable properties.