Junction Tree Variational Autoencoder for Molecular Graph Generation

02/12/2018 ∙ by Wengong Jin, et al. ∙ 0

We seek to automate the design of molecules based on specific chemical properties. In computational terms, this task involves continuous embedding and generation of molecular graphs. Our primary contribution is the direct realization of molecular graphs, a task previously approached by generating linear SMILES strings instead of graphs. Our junction tree variational autoencoder generates molecular graphs in two phases, by first generating a tree-structured scaffold over chemical substructures, and then combining them into a molecule with a graph message passing network. This approach allows us to incrementally expand molecules while maintaining chemical validity at every step. We evaluate our model on multiple tasks ranging from molecular generation to optimization. Across these tasks, our model outperforms previous state-of-the-art baselines by a significant margin.

READ FULL TEXT VIEW PDF
POST COMMENT

Comments

There are no comments yet.

Authors

page 1

page 2

page 3

page 4

Code Repositories

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

The key challenge of drug discovery is to find target molecules with desired chemical properties. Currently, this task takes years of development and exploration by expert chemists and pharmacologists. Our ultimate goal is to automate this process. From a computational perspective, we decompose the challenge into two complementary subtasks: learning to represent molecules in a continuous manner that facilitates the prediction and optimization of their properties (encoding); and learning to map an optimized continuous representation back into a molecular graph with improved properties (decoding). While deep learning has been extensively investigated for molecular graph encoding (Duvenaud et al., 2015; Kearnes et al., 2016; Gilmer et al., 2017), the harder combinatorial task of molecular graph generation from latent representation remains under-explored.

Prior work on drug design formulated the graph generation task as a string generation problem (Gómez-Bombarelli et al., 2016; Kusner et al., 2017) in an attempt to side-step direct generation of graphs. Specifically, these models start by generating SMILES (Weininger, 1988), a linear string notation used in chemistry to describe molecular structures. SMILES strings can be translated into graphs via deterministic mappings (e.g., using RDKit (Landrum, 2006)). However, this design has two critical limitations. First, the SMILES representation is not designed to capture molecular similarity. For instance, two molecules with similar chemical structures may be encoded into markedly different SMILES strings (e.g., Figure 1). This prevents generative models like variational autoencoders from learning smooth molecular embeddings. Second, essential chemical properties such as molecule validity are easier to express on graphs rather than linear SMILES representations. We hypothesize that operating directly on graphs improves generative modeling of valid chemical structures.

Figure 1: Two almost identical molecules with markedly different canonical SMILES in RDKit. The edit distance between two strings is 22 (50.5% of the whole sequence).

Our primary contribution is a new generative model of molecular graphs. While one could imagine solving the problem in a standard manner – generating graphs node by node – the approach is not ideal for molecules. This is because creating molecules atom by atom would force the model to generate chemically invalid intermediaries (see, e.g., Figure 2), delaying validation until a complete graph is generated. Instead, we propose to generate molecular graphs in two phases by exploiting valid subgraphs as components. The overall generative approach, cast as a junction tree variational autoencoder, first generates a tree structured object (a junction tree) whose role is to represent the scaffold of subgraph components and their coarse relative arrangements. The components are valid chemical substructures automatically extracted from the training set using tree decomposition and are used as building blocks. In the second phase, the subgraphs (nodes in the tree) are assembled together into a coherent molecular graph.

We evaluate our model on multiple tasks ranging from molecular generation to optimization of a given molecule according to desired properties. As baselines, we utilize state-of-the-art SMILES-based generation approaches (Kusner et al., 2017; Dai et al., 2018). We demonstrate that our model produces 100% valid molecules when sampled from a prior distribution, outperforming the top performing baseline by a significant margin. In addition, we show that our model excels in discovering molecules with desired properties, yielding a 30% relative gain over the baselines.

Figure 2: Comparison of two graph generation schemes: Structure by structure approach is preferred as it avoids invalid intermediate states (marked in red) encountered in node by node approach.

2 Junction Tree Variational Autoencoder

Our approach extends the variational autoencoder (Kingma & Welling, 2013) to molecular graphs by introducing a suitable encoder and a matching decoder. Deviating from previous work (Gómez-Bombarelli et al., 2016; Kusner et al., 2017)

, we interpret each molecule as having been built from subgraphs chosen out of a vocabulary of valid components. These components are used as building blocks both when encoding a molecule into a vector representation as well as when decoding latent vectors back into valid molecular graphs. The key advantage of this view is that the decoder can realize a valid molecule piece by piece by utilizing the collection of valid components and how they interact, rather than trying to build the molecule atom by atom through chemically invalid intermediaries (Figure 

2). An aromatic bond, for example, is chemically invalid on its own unless the entire aromatic ring is present. It would be therefore challenging to learn to build rings atom by atom rather than by introducing rings as part of the basic vocabulary.

Our vocabulary of components, such as rings, bonds and individual atoms, is chosen to be large enough so that a given molecule can be covered by overlapping components or clusters of atoms. The clusters serve the role analogous to cliques in graphical models, as they are expressive enough that a molecule can be covered by overlapping clusters without forming cluster cycles. In this sense, the clusters serve as cliques in a (non-optimal) triangulation of the molecular graph. We form a junction tree of such clusters and use it as the tree representation of the molecule. Since our choice of cliques is constrained a priori, we cannot guarantee that a junction tree exists with such clusters for an arbitrary molecule. However, our clusters are built on the basis of the molecules in the training set to ensure that a corresponding junction tree can be found. Empirically, our clusters cover most of the molecules in the test set.

The original molecular graph and its associated junction tree offer two complementary representations of a molecule. We therefore encode the molecule into a two-part latent representation where encodes the tree structure and what the clusters are in the tree without fully capturing how exactly the clusters are mutually connected. encodes the graph to capture the fine-grained connectivity. Both parts are created by tree and graph encoders and . The latent representation is then decoded back into a molecular graph in two stages. As illustrated in Figure 3, we first reproduce the junction tree using a tree decoder based on the information in . Second, we predict the fine grain connectivity between the clusters in the junction tree using a graph decoder to realize the full molecular graph. The junction tree approach allows us to maintain chemical feasibility during generation.

Figure 3: Overview of our molecule generation paradigm: A molecular graph is first decomposed into its junction tree , where each colored node in the tree represents a substructure in the molecule. We then encode both the tree and graph into their latent embeddings and . To decode the molecule, we first reconstruct junction tree from , and then assemble nodes in the tree back to original molecule, guided by .

Notation A molecular graph is defined as where is the set of atoms (vertices) and the set of bonds (edges). Let be the neighbor of

. We denote sigmoid function as

and ReLU function as

. We use for nodes in the tree and for nodes in the graph.

2.1 Junction Tree

A tree decomposition maps a graph into a junction tree by contracting certain vertices into a single node so that becomes cycle-free. Formally, given a graph , a junction tree is a connected labeled tree whose node set is and edge set is . Each node or cluster is an induced subgraph of , satisfying the following constraints:

  1. The union of all clusters equals . That is, and .

  2. Running intersection: For all clusters and , if is on the path from to .

Viewing induced subgraphs as cluster labels, junction trees are labeled trees with label vocabulary . By our molecule tree decomposition, contains only cycles (rings) and single edges. Thus the vocabulary size is limited ( for a standard dataset with 250K molecules).

Tree Decomposition of Molecules Here we present our tree decomposition algorithm tailored for molecules, which finds its root in chemistry (Rarey & Dixon, 1998). Our cluster vocabulary includes chemical structures such as bonds and rings (Figure 3). Given a graph , we first find all its simple cycles, and its edges not belonging to any cycles. Two simple rings are merged together if they have more than two overlapping atoms, as they constitute a specific structure called bridged compounds (Clayden et al., 2001). Each of those cycles or edges is considered as a cluster. Next, a cluster graph is constructed by adding edges between all intersecting clusters. Finally, we select one of its spanning trees as the junction tree of (Figure 3). As a result of ring merging, any two clusters in the junction tree have at most two atoms in common, facilitating efficient inference in the graph decoding phase. The detailed procedure is described in the supplementary.

2.2 Graph Encoder

We first encode the latent representation of by a graph message passing network (Dai et al., 2016; Gilmer et al., 2017). Each vertex has a feature vector indicating the atom type, valence, and other properties. Similarly, each edge has a feature vector indicating its bond type, and two hidden vectors and denoting the message from to and vice versa. Due to the loopy structure of the graph, messages are exchanged in a loopy belief propagation fashion:

(1)

where is the message computed in -th iteration, initialized with . After steps of iteration, we aggregate those messages as the latent vector of each vertex, which captures its local graphical structure:

(2)

The final graph representation is . The mean

and log variance

of the variational posterior approximation are computed from

with two separate affine layers.

is sampled from a Gaussian .

2.3 Tree Encoder

We similarly encode with a tree message passing network. Each cluster

is represented by a one-hot encoding

representing its label type. Each edge is associated with two message vectors and . We pick an arbitrary leaf node as the root and propagate messages in two phases. In the first bottom-up phase, messages are initiated from the leaf nodes and propagated iteratively towards root. In the top-down phase, messages are propagated from the root to all the leaf nodes. Message is updated as:

(3)

where

is a Gated Recurrent Unit 

(Chung et al., 2014; Li et al., 2015) adapted for tree message passing:

(4)
(5)
(6)
(7)
(8)

The message passing follows the schedule where is computed only when all its precursors have been computed. This architectural design is motivated by the belief propagation algorithm over trees and is thus different from the graph encoder.

After the message passing, we obtain the latent representation of each node by aggregating its inward messages:

(9)

The final tree representation is , which encodes a rooted tree . Unlike the graph encoder, we do not apply node average pooling because it confuses the tree decoder which node to generate first. is sampled in a similar way as in the graph encoder. For simplicity, we abbreviate as from now on.

This tree encoder plays two roles in our framework. First, it is used to compute , which only requires the bottom-up phase of the network. Second, after a tree is decoded from , it is used to compute messages over the entire , to provide essential contexts of every node during graph decoding. This requires both top-down and bottom-up phases. We will elaborate this in section 2.5.

2.4 Tree Decoder

We decode a junction tree from its encoding with a tree structured decoder. The tree is constructed in a top-down fashion by generating one node at a time. As illustrated in Figure 4, our tree decoder traverses the entire tree from the root, and generates nodes in their depth-first order. For every visited node, the decoder first makes a topological prediction: whether this node has children to be generated. When a new child node is created, we predict its label and recurse this process. Recall that cluster labels represent subgraphs in a molecule. The decoder backtracks when a node has no more children to generate.

At each time step, a node receives information from other nodes in the current tree for making those predictions. The information is propagated through message vectors when trees are incrementally constructed. Formally, let be the edges traversed in a depth first traversal over , where as each edge is traversed in both directions. The model visits node at time . Let be the first edges in . The message is updated through previous messages:

(10)

where is the same recurrent unit as in the tree encoder.

Figure 4: Illustration of the tree decoding process. Nodes are labeled in the order in which they are generated. 1) Node 2 expands child node 4 and predicts its label with message . 2) As node 4 is a leaf node, decoder backtracks and computes message . 3) Decoder continues to backtrack as node 2 has no more children. 4) Node 1 expands node 5 and predicts its label.
0:  Latent representation
1:  Initialize: Tree
2:  function SampleTree(
3:     Set all cluster labels that are chemically compatible with node and its current neighbors.
4:     Set

with probability

. Eq.(11)
5:     if  and  then
6:        Create a node and add it to tree .
7:        Sample the label of node from . Eq.(12)
8:        SampleTree()
9:     end if
10:  end function
Algorithm 1 Tree decoding at sampling time

Topological Prediction When the model visits node , it makes a binary prediction on whether it still has children to be generated. We compute this probability by combining , node features and inward messages via a one hidden layer network followed by a sigmoid function:

(11)

Label Prediction When a child node is generated from its parent , we predict its node label with

(12)

where is a distribution over label vocabulary . When is a root node, its parent is a virtual node and .

Learning The tree decoder aims to maximize the likelihood . Let and be the ground truth topological and label values, the decoder minimizes the following cross entropy loss:111The node ordering is not unique as the order within sibling nodes is ambiguous. In this paper we train our model with one ordering and leave this issue for future work.

(13)

Similar to sequence generation, during training we perform teacher forcing: after topological and label prediction at each step, we replace them with their ground truth so that the model makes predictions given correct histories.

Decoding & Feasibility Check Algorithm 1 shows how a tree is sampled from . The tree is constructed recursively guided by topological predictions without any external guidance used in training. To ensure the sampled tree could be realized into a valid molecule, we define set to be cluster labels that are chemically compatible with node and its current neighbors. When a child node is generated from node , we sample its label from with a renormalized distribution over by masking out invalid labels.

2.5 Graph Decoder

The final step of our model is to reproduce a molecular graph that underlies the predicted junction tree

. Note that this step is not deterministic since there are potentially many molecules that correspond to the same junction tree. The underlying degree of freedom pertains to how neighboring clusters

and are attached to each other as subgraphs. Our goal here is to assemble the subgraphs (nodes in the tree) together into the correct molecular graph.

Figure 5: Decode a molecule from a junction tree. 1) Ground truth molecule . 2) Predicted junction tree . 3) We enumerate different combinations between red cluster and its neighbors. Crossed arrows indicate combinations that lead to chemically infeasible molecules. Note that if we discard tree structure during enumeration (i.e., ignoring subtree A), the last two candidates will collapse into the same molecule. 4) Rank subgraphs at each node. The final graph is decoded by putting together all the predicted subgraphs (dashed box).

Let be the set of graphs whose junction tree is . Decoding graph from is a structured prediction:

(14)

where is a scoring function over candidate graphs. We only consider scoring functions that decompose across the clusters and their neighbors. In other words, each term in the scoring function depends only on how a cluster is attached to its neighboring clusters , in the tree . The problem of finding the highest scoring graph – the assembly task – could be cast as a graphical model inference task in a model induced by the junction tree. However, for efficiency reasons, we will assemble the molecular graph one neighborhood at a time, following the order in which the tree itself was decoded. In other words, we start by sampling the assembly of the root and its neighbors according to their scores. Then we proceed to assemble the neighbors and their associated clusters (removing the degrees of freedom set by the root assembly), and so on.

It remains to be specified how each neighborhood realization is scored. Let be the subgraph resulting from a particular merging of cluster in the tree with its neighbors , . We score as a candidate subgraph by first deriving a vector representation and then using as the subgraph score. To this end, let specify atoms in the candidate subgraph and let if and if . The indices are used to mark the position of the atoms in the junction tree, and to retrieve messages summarizing the subtree under along the edge obtained by running the tree encoding algorithm. The neural messages pertaining to the atoms and bonds in subgraph are obtained and aggregated into , similarly to the encoding step, but with different (learned) parameters:

(15)

The major difference from Eq. (1) is that we augment the model with tree messages derived by running the tree encoder over the predicted tree . provides a tree dependent positional context for bond (illustrated as subtree A in Figure 5).

Learning The graph decoder parameters are learned to maximize the log-likelihood of predicting correct subgraphs of the ground true graph at each tree node:

(16)

where is the set of possible candidate subgraphs at tree node . During training, we again apply teacher forcing, i.e. we feed the graph decoder with ground truth trees as input.

Complexity By our tree decomposition, any two clusters share at most two atoms, so we only need to merge at most two atoms or one bond. By pruning chemically invalid subgraphs and merging isomorphic graphs, on average when tested on a standard ZINC drug dataset. The computational complexity of JT-VAE is therefore linear in the number of clusters, scaling nicely to large graphs.

3 Experiments

Figure 6: Left: Random molecules sampled from prior distribution . Right: Visualization of the local neighborhood of a molecule in the center. Three molecules highlighted in red dashed box have the same tree structure as the center molecule, but with different graph structure as their clusters are combined differently. The same phenomenon emerges in another group of molecules (blue dashed box).

Our evaluation efforts measure various aspects of molecular generation. The first two evaluations follow previously proposed tasks (Kusner et al., 2017). We also introduce a third task — constrained molecule optimization.

  • [leftmargin=*]

  • Molecule reconstruction and validity We test the VAE models on the task of reconstructing input molecules from their latent representations, and decoding valid molecules when sampling from prior distribution. (Section 3.1)

  • Bayesian optimization Moving beyond generating valid molecules, we test how the model can produce novel molecules with desired properties. To this end, we perform Bayesian optimization in the latent space to search molecules with specified properties. (Section 3.2)

  • Constrained molecule optimization The task is to modify given molecules to improve specified properties, while constraining the degree of deviation from the original molecule. This is a more realistic scenario in drug discovery, where development of new drugs usually starts with known molecules such as existing drugs (Besnard et al., 2012). Since it is a new task, we cannot compare to any existing baselines. (Section 3.3)

Below we describe the data, baselines and model configuration that are shared across the tasks. Additional setup details are provided in the task-specific sections.

Data We use the ZINC molecule dataset from Kusner et al. (2017) for our experiments. It contains about 250K drug molecules extracted from the ZINC database (Sterling & Irwin, 2015). We follow the same training/testing split as the previous work.

Baselines We compare our approach with SMILES-based baselines: 1) Character VAE (CVAE) (Gómez-Bombarelli et al., 2016) which generates SMILES strings character by character; 2) Grammar VAE (GVAE) (Kusner et al., 2017) that generates SMILES following syntactic constraints given by a context-free grammar; 3) Syntax-directed VAE (SD-VAE) (Dai et al., 2018) that incorporates both syntactic and semantic constraints of SMILES via attribute grammar. For molecule generation task, we also compare with GraphVAE (Simonovsky & Komodakis, 2018) that directly generates atom labels and adjacency matrices of graphs.

Model Configuration To be comparable with the above baselines, we set the latent space dimension as 56, i.e., the tree and graph representation and have 28 dimensions each. Full training details and model configurations are provided in the appendix.

3.1 Molecule Reconstruction and Validity

Setup

The first task is to reconstruct and sample molecules from latent space. Since both encoding and decoding process are stochastic, we estimate reconstruction accuracy by Monte Carlo method used in

(Kusner et al., 2017): Each molecule is encoded 10 times and each encoding is decoded 10 times. We report the portion of the 100 decoded molecules that are identical to the input molecule. For a fair comparison, we define two molecules as identical if they have the same SMILES string. At testing time, we convert all generated graphs to SMILES using RDKit.

To compute validity, we sample 1000 latent vectors from the prior distribution , and decode each of these vectors 100 times. We report the percentage of decoded molecules that are chemically valid (checked by RDKit).

Results Table 1 shows that JT-VAE outperforms previous models in molecule reconstruction, and always produces valid molecules when sampled from prior distribution. These sampled molecules have non-trivial structures such as simple chains (Figure 6). We further sampled 5000 molecules from prior and found they are all distinct from the training set. Thus our model is not a simple memorization.

Method Reconstruction Validity
CVAE 44.6% 0.7%
GVAE 53.7% 7.2%
SD-VAE222The SD-VAE result is copied from Table 1 in Dai et al. (2018). 76.2% 43.5%
GraphVAE - 13.5%
JT-VAE 76.7% 100.0%
Table 1: Reconstruction accuracy and prior validity results. Baseline results are copied from Kusner et al. (2017); Dai et al. (2018); Simonovsky & Komodakis (2018).

Analysis We qualitatively examine the latent space of JT-VAE by visualizing the neighborhood of molecules. Given a molecule, we follow the method in Kusner et al. (2017) to construct a grid visualization of its neighborhood. For comparison, we select the same molecule visualized in Dai et al. (2018). Figure 6 shows the local neighborhood of this molecule. Compared to the figure in Dai et al. (2018), our neighborhood does not contain molecules with huge rings (with more than 7 atoms), which rarely occur in the dataset. We also highlight two groups of closely resembling molecules that have identical tree structures but vary only in how clusters are attached together. This demonstrates the smoothness of learned molecular embeddings.

3.2 Bayesian Optimization

Setup The second task is to produce novel molecules with desired properties. Following (Kusner et al., 2017), our target chemical property is octanol-water partition coefficients (logP) penalized by the synthetic accessibility (SA) score and number of long cycles.333 where counts the number of rings that have more than six atoms. To perform Bayesian optimization (BO), we first train a VAE and associate each molecule with a latent vector, given by the mean of the variational encoding distribution. After the VAE is learned, we train a sparse Gaussian process (SGP) to predict

given its latent representation. Then we perform five iterations of batched BO using the expected improvement heuristic.

For comparison, we report 1) the predictive performance of SGP trained on latent encodings learned by different VAEs, measured by log-likelihood (LL) and root mean square error (RMSE) with 10-fold cross validation. 2) The top-3 molecules found by BO under different models.

Results As shown in Table 2, JT-VAE finds molecules with significantly better scores than previous methods. Figure 7 lists the top-3 best molecules found by JT-VAE. In fact, JT-VAE finds over 50 molecules with scores over 3.50 (the second best molecule proposed by SD-VAE). Moreover, the SGP yields better predictive performance when trained on JT-VAE embeddings (Table 3).

Method 1st 2nd 3rd
CVAE 1.98 1.42 1.19
GVAE 2.94 2.89 2.80
SD-VAE 4.04 3.50 2.96
JT-VAE 5.30 4.93 4.49
Table 2: Best molecule property scores found by each method. Baseline results are from Kusner et al. (2017); Dai et al. (2018).
Figure 7: Best three molecules and their property scores found by JT-VAE using Bayesian optimization.
Method LL RMSE
CVAE
GVAE
SD-VAE
JT-VAE
Table 3: Predictive performance of sparse Gaussian Processes trained on different VAEs. Baseline results are copied from Kusner et al. (2017) and Dai et al. (2018).
Improvement Similarity Success
0.0 97.5%
0.2 97.1%
0.4 83.6%
0.6 46.4%
Table 4:

Constrained optimization result of JT-VAE: mean and standard deviation of property improvement, molecular similarity and success rate under constraints

with varied .
Figure 8: A molecule modification that yields an improvement of 4.0 with molecular similarity 0.617 (modified part is in red).

3.3 Constrained Optimization

Setup The third task is to perform molecule optimization in a constrained scenario. Given a molecule , the task is to find a different molecule that has the highest property value with the molecular similarity for some threshold . We use Tanimoto similarity with Morgan fingerprint (Rogers & Hahn, 2010) as the similarity metric, and penalized logP coefficient as our target chemical property. For this task, we jointly train a property predictor (parameterized by a feed-forward network) with JT-VAE to predict from the latent embedding of . To optimize a molecule , we start from its latent representation, and apply gradient ascent in the latent space to improve the predicted score , similar to (Mueller et al., 2017). After applying gradient steps, molecules are decoded from resulting latent trajectories, and we report the molecule with the highest that satisfies the similarity constraint. A modification succeeds if one of the decoded molecules satisfies the constraint and is distinct from the original.

To provide the greatest challenge, we selected 800 molecules with the lowest property score from the test set. We report the success rate (how often a modification succeeds), and among success cases the average improvement and molecular similarity between the original and modified molecules and .

Results Our results are summarized in Table 4. The unconstrained scenario () has the best average improvement, but often proposes dissimilar molecules. When we tighten the constraint to , about 80% of the time our model finds similar molecules, with an average improvement . This also demonstrates the smoothness of the learned latent space. Figure 8 illustrates an effective modification resulting in a similar molecule with great improvement.

4 Related Work

Molecule Generation Previous work on molecule generation mostly operates on SMILES strings. Gómez-Bombarelli et al. (2016); Segler et al. (2017) built generative models of SMILES strings with recurrent decoders. Unfortunately, these models could generate invalid SMILES that do not result in any molecules. To remedy this issue, Kusner et al. (2017); Dai et al. (2018)

complemented the decoder with syntactic and semantic constraints of SMILES by context free and attribute grammars, but these grammars do not fully capture chemical validity. Other techniques such as active learning 

(Janz et al., 2017)

and reinforcement learning 

(Guimaraes et al., 2017) encourage the model to generate valid SMILES through additional training signal. Very recently, Simonovsky & Komodakis (2018) proposed to generate molecular graphs by predicting their adjacency matrices, and Li et al. (2018) generated molecules node by node. In comparison, our method enforces chemical validity and is more efficient due to the coarse-to-fine generation.

Graph-structured Encoders

The neural network formulation on graphs was first proposed by

Gori et al. (2005); Scarselli et al. (2009), and later enhanced by Li et al. (2015) with gated recurrent units. For recurrent architectures over graphs, Lei et al. (2017) designed Weisfeiler-Lehman kernel network inspired by graph kernels. Dai et al. (2016) considered a different architecture where graphs were viewed as latent variable graphical models, and derived their model from message passing algorithms. Our tree and graph encoder are closely related to this graphical model perspective, and to neural message passing networks (Gilmer et al., 2017). For convolutional architectures, Duvenaud et al. (2015) introduced a convolution-like propagation on molecular graphs, which was generalized to other domains by Niepert et al. (2016). Bruna et al. (2013); Henaff et al. (2015) developed graph convolution in spectral domain via graph Laplacian. For applications, graph neural networks are used in semi-supervised classification (Kipf & Welling, 2016)

, computer vision 

(Monti et al., 2016), and chemical domains (Kearnes et al., 2016; Schütt et al., 2017; Jin et al., 2017).

Tree-structured Models Our tree encoder is related to recursive neural networks and tree-LSTM (Socher et al., 2013; Tai et al., 2015; Zhu et al., 2015). These models encode tree structures where nodes in the tree are bottom-up transformed into vector representations. In contrast, our model propagates information both bottom-up and top-down.

On the decoding side, tree generation naturally arises in natural language parsing (Dyer et al., 2016; Kiperwasser & Goldberg, 2016). Different from our approach, natural language parsers have access to input words and only predict the topology of the tree. For general purpose tree generation, Vinyals et al. (2015); Aharoni & Goldberg (2017) applied recurrent networks to generate linearized version of trees, but their architectures were entirely sequence-based. Dong & Lapata (2016); Alvarez-Melis & Jaakkola (2016) proposed tree-based architectures that construct trees top-down from the root. Our model is most closely related to Alvarez-Melis & Jaakkola (2016) that disentangles topological prediction from label prediction, but we generate nodes in a depth-first order and have additional steps that propagate information bottom-up. This forward-backward propagation also appears in Parisotto et al. (2016), but their model is node based whereas ours is based on message passing.

5 Conclusion

In this paper we present a junction tree variational autoencoder for generating molecular graphs. Our method significantly outperforms previous work in molecule generation and optimization. For future work, we attempt to generalize our method for general low-treewidth graphs.

Acknowledgement

We thank Jonas Mueller, Chengtao Li, Tao Lei and MIT NLP Group for their helpful comments. This work was supported by the DARPA Make-It program under contract ARO W911NF-16-2-0023.

References

Appendix A Tree Decomposition

Algorithm 2 presents our tree decomposition of molecules. and contain non-ring bonds and simple rings respectively. Simple rings are extracted via RDKit’s GetSymmSSSR function. We then merge rings that share three or more atoms as they form bridged compounds. We note that the junction tree of a molecule is not unique when its cluster graph contains cycles. This introduces additional uncertainty for our probabilistic modeling. To reduce such variation, for any of the three (or more) intersecting bonds, we add their intersecting atom as a cluster and remove the cycle connecting them in the cluster graph. Finally, we construct a junction tree as the maximum spanning tree of a cluster graph . Note that we assign an large weight over edges involving clusters in to ensure no edges in any cycles will be selected into the junction tree.

   the set of bonds that do not belong to any rings.
   the set of simple rings of .
  for  in  do
     Merge rings into one ring if they share more than two atoms (bridged rings).
  end for
   atoms being the intersection of three or more clusters in .
  
  . Set if or , and otherwise.
  Return The maximum spanning tree over cluster graph .
Algorithm 2 Tree decomposition of molecule
Figure 9: Illustration of tree decomposition and sample of cluster label vocabulary.

Appendix B Stereochemistry

Though usually presented as two-dimensional graphs, molecules are three-dimensional objects, i.e. molecules are defined not only by its atom types and bond connections, but also the spatial configuration between atoms (chiral atoms and cis-trans isomerism). Stereoisomers are molecules that have the same 2D structure, but differ in the 3D orientations of their atoms in space. We note that stereochemical feasibility could not be simply encoded as context free or attribute grammars.

Empirically, we found it more efficient to predict the stereochemical configuration separately from the molecule generation. Specifically, the JT-VAE first generates the 2D structure of a molecule , following the same procedure described in section 2. Then we generate all its stereoisomers using RDKit’s EnumerateStereoisomers function, which identifies atoms that could be chiral. For each isomer , we encode its graph representation

with the graph encoder and compute their cosine similarity

(note that is stochastic). We reconstruct the final 3D structure by picking the stereoisomer . Since on average only few atoms could have stereochemical variations, this post ranking process is very efficient. Combining this with tree and graph generation, the molecule reconstruction loss becomes

(17)

Appendix C Training Details

By applying tree decomposition over 240K molecules in ZINC dataset, we collected our vocabulary set of size

. The hidden state dimension is 450 for all modules in JT-VAE and the latent bottleneck dimension is 56. For the graph encoder, the initial atom features include its atom type, degree, its formal charge and its chiral configuration. Bond feature is a concatenation of its bond type, whether the bond is in a ring, and its cis-trans configuration. For our tree encoder, we represent each cluster with a neural embedding vector, similar to word embedding for words. The tree and graph decoder use the same feature setting as encoders. The graph encoder and decoder runs three iterations of neural message passing. For fair comparison to SMILES based method, we minimized feature engineering. We use PyTorch to implement all neural components and RDKit to process molecules.

Appendix D More Experimental Results

Sampled Molecules Note that a degenerate model could also achieve 100% prior validity by keep generating simple structures like chains. To prove that our model does not converge to such trivial solutions, we randomly sample and plot 250 molecules from prior distribution . As shown in Figure 10, our sampled molecules present rich variety and structural complexity. This demonstrates the soundness of the prior validity improvement of our model.

Neighborhood Visualization Given a molecule, we follow Kusner et al. (2017) to construct a grid visualization of its neighborhood. Specifically, we encode a molecule into the latent space and generate two random orthogonal unit vectors as two axis of a grid. Moving in combinations of these directions yields a set of latent vectors and we decode them into corresponding molecules. In Figure 11 and 12, we visualize the local neighborhood of two molecules presented in Dai et al. (2018). Figure 11 visualizes the same molecule in Figure 6, but with wider neighborhood ranges.

Bayesian Optimization We directly used open sourced implementation in Kusner et al. (2017) for Bayesian optimization (BO). Specifically, we train a sparse Gaussian process with 500 inducing points to predict properties of molecules. Five iterations of batch BO with expected improvement heuristic is used to propose new latent vectors. In each iteration, 50 latent vectors are proposed, from which molecules are decoded and added to the training set for next iteration. We perform 10 independent runs and aggregate results. In Figure 13, we present the top 50 molecules found among 10 runs using JT-VAE. Following Kusner et al.’s implementation, the scores reported are normalized to zero mean and unit variance by the mean and variance computed from training set.

Constrained Optimization For this task, a property predictor is trained jointly with VAE to predict from the latent embedding of . is a feed-forward network with one hidden layer of dimension 450 followed by activation. To optimize a molecule , we start with its mean encoding and apply 80 gradient ascent steps: with . 80 molecules are decoded from latent vectors and their property is calculated. Molecular similarity is calculated via Morgan fingerprint of radius 2 with Tanimoto similarity. For each molecule , we report the best modified molecule with for some threshold . In Figure 14, we present three groups of modification examples with . For each group, we present top three pairs that leads to best improvement as well as one pair decreased property (). This is caused by inaccurate property prediction. From Figure 14, we can see that tighter similarity constraint forces the model to preserve the original structure.

Figure 10: 250 molecules sampled from prior distribution .
Figure 11: Neighborhood visualization of molecule C[C@H]1CC(Nc2cncc(-c3nncn3C)c2)C[C@H](C)C1.
Figure 12: Neighborhood visualization of molecule COc1cc(OC)cc([C@H]2CC[NH+](CCC(F)(F)F)C2)c1.
Figure 13: Top 50 molecules found by Bayesian optimization using JT-VAE.
Figure 14: Row 1-3: Molecule modification results with similarity constraint . For each group, we plot the top three pairs that leads to actual property improvement, and one pair with decreased property. We can see that tighter similarity constraint forces the model to preserve the original structure.