A Generative Model for Molecular Distance Geometry

Introduction

To compute a molecular property for a molecule, one must sample from $p(\mathbf{x})$. The main approach is to start with one conformation and make small changes to it over time, e.g., by using Markov chain Monte Carlo (MCMC) or molecular dynamics (MD). These methods can be used to accurately sample equilibrium states of molecules, but they become computationally expensive for larger ones (Shim & MacKerell 2011; Ballard et al. 2015; De Vivo et al. 2016). Other heuristic approaches exist in which distances between atoms are set to fixed idealized values (Havel 2002; Blaney & Dixon 2007). Several methods based on statistical learning have also recently been developed to tackle the issue of conformation generation. However, they are mainly geared towards studying proteins and their folding dynamics (AlQuraishi 2019). Some of these models are not targeting a distribution over conformations but the most stable folded configuration, e.g. AlphaFold (Senior et al. 2020), while others are not transferable between different molecules (Lemke & Peter 2019; Noé et al. 2019).

This work includes the following key contributions:

We introduce a novel probabilistic model for learning conformational distributions of molecules with graph neural networks.

We create a new, challenging benchmark for conformation generation, which is made publicly available. To the best of our knowledge, this is the first benchmark of this kind.

By combining a conditional variational autoencoder (CVAE) with an Euclidean distance geometry (EDG) algorithm we present a state-of-the-art approach for generating one-shot samples of molecular conformations for unseen molecules that is independent of their size and shape.

We develop a rigorous experimental approach for evaluating and comparing the accuracy of conformation generation methods based on the mean maximum deviation distance metric.

We show how this generative model can be used as a proposal distribution in an importance sampling (IS) scheme to estimate molecular properties.

Method

Our goal is to build a statistical model that generates molecular conformations in a one-shot fashion from a molecule’s graph representation. First, we describe how a molecule’s conformation can be represented by a set of pairwise distances between atoms and why this presentation is advantageous over one in Cartesian coordinates (Section 2.1). Second, we present a generative model in Section 2.2 that will generate sets of atomic distances for a given molecular graph. Third, we explain in Section 2.3 how a set of predicted distances can be transformed into a molecular conformation and why this transformation is necessary. Finally, we detail in Section 2.4 how our generative model can be used as a proposal distribution in an IS scheme to estimate molecular properties.

We assume that, given a molecular graph $G$, one can represent one of its conformations $\mathbf{x}$ by a set of atomic distances $\mathbf{d}=\{d_{k}\}_{k=1}^{N_{e}}$, where $d_{k}=|\mathbf{r}_{r_{k}}-\mathbf{r}_{s_{k}}|$ is the Euclidean distance between the positions of the atoms $r_{k}$ and $s_{k}$ in this conformation. As the set of edges between the bonded atoms ($E_{\text{bond}}$) alone would not suffice to describe a conformation, we expand the traditional graph representation of a molecule by adding auxiliary edges to obtain an extended graph $\mathcal{G}$. Auxiliary edges between atoms that are second neighbors in the original graph $G$ fix angles between atoms, and those between third neighbors fix dihedral angles (denoted $E_{\text{angle}}$ and $E_{\text{dihedral}}$, respectively). In this work, $E_{\text{angle}}$ are added between nodes in $\mathcal{G}$ which are second neighbors in $G$. After all $E_{\text{angle}}$ have been added, additional edges are added to $\mathcal{G}$ from a node $v$ to a randomly chosen third neighbor of $v$ in $G$ if $v$ has less then three neighbors in $\mathcal{G}$. Therefore, a graph $G$ can give rise to multiple different extend graphs $\mathcal{G}$. In Fig. 2, the process of extending the molecular graph and the extraction of $\mathbf{d}$ from $\mathbf{x}$ and $\mathcal{G}$ are illustrated.

A key advantage of a representation in terms of distances is its invariance to rotation and translation; by contrast, Cartesian coordinates depend on the (arbitrary) choice of origin, for example. In addition, it reflects pair-wise physical interactions and their generally local nature. Auxiliary edges can be placed between higher-order neighbors depending on how far the physical interactions dominating the potential energy of the system reach.

We have a set of $N_{G}$ pairs, $\{G_{i},x_{i}\}_{i=1}^{N_{G}}$, consisting of a molecular graph and a conformation. With the protocol described above, we convert each pair into a pair of an extended molecular graph together with a set of distances $\mathbf{d}$ to obtain $\{\mathcal{G}_{i},\mathbf{d}_{i}\}_{i=1}^{N_{G}}$. With this data, we will train a generative model which we detail in the following section.

2 Generative Model

Here, $q_{\phi}(\mathbf{z}|\mathbf{d},\mathcal{G})$ and $p_{\theta}(\mathbf{d}|\mathbf{z},\mathcal{G})$ are Gaussian distributions, the mean and variance of which are modeled by two artificial neural networks. At the center of this model are message-passing neural networks (MPNNs) (Gilmer et al. 2017). In short, an MPNN is a convolutional neural network that allows end-to-end learning of prediction pipelines whose inputs are graphs of arbitrary size and shape. In a convolution, neighboring nodes exchange so-called messages between neighbors to update their attributes. Edges update their attributes with the features of the nodes they are connecting. The MPNN is a well-studied technique that achieves state-of-the-art performance in representation learning for molecules (Kipf & Welling 2017; Duvenaud et al. 2015; Kearnes et al. 2016; Schütt et al. 2017b; Gilmer et al. 2017; Kusner et al. 2017; Bradshaw et al. 2019a).

The sets of parameters in the encoder and decoder, $\phi$ and $\theta$ (i.e., parameters in $F_{\text{enc},v}$, $F_{\text{enc},e}$, $\{\text{MP}^{(t)}_{\text{enc}}\}_{t=1}^{T}$, $R_{\text{enc}}$, $F_{\text{dec},v}$, $F_{\text{dec},e}$, $\{\text{MP}^{(t)}_{\text{dec}}\}_{t=1}^{T}$, $R_{\text{dec}}$), respectively, are optimized by maximizing the evidence lower bound (ELBO):

where the prior $p_{\theta}(\mathbf{z}|\mathcal{G})$ consists of factorized standard Gaussians. The optimal values for the hyperparameters for the network dimensions, number of message passes, batch size, and learning rate of the Adam optimizer (Kingma & Ba 2014) were manually tuned by maximizing the validation performance (ELBO) and are reported in the Appendix.

3 Conformation Generation through Euclidean Distance Geometry

To compute molecular properties, quantum-chemical methods need to be employed which require the input, i.e., the molecule, to be in Cartesian coordinates. Even though quantum-chemical methods require the input to be in Cartesian coordinates, calculated properties, such as the energy, are invariant under translation and rotation. Therefore, we use an EDG algorithm to translate the set of distances $\{d_{k}\}_{k=1}^{N_{e}}$ to a set of atomic coordinates $\{\mathbf{r}_{i}\}_{i=1}^{N_{v}}$. There are additional constraints due to chirality. However, since they are given by $\mathcal{G}$ and are fixed, they are not modeled by our method.

EDG is the mathematical basis for a geometric theory of molecular conformation. In the field of machine learning, Weinberger & Saul 2006 used it for learning image manifolds, Tenenbaum et al. 2000 for image understanding and handwriting recognition, Jain & Saul 2004 for speech and music, and Demaine et al. 2009 for music and musical rhythms. An EDG description of a molecular system consists of a list of lower and upper bounds on the distances between pairs of atoms $\{(d_{k,\text{min}},d_{k,\text{max}})\}_{k=1}^{N_{e}}$. Here, $p_{\theta}(\mathbf{d}|\mathbf{z},\mathcal{G})$ is used to model these bounds, namely, we set the bounds to $\{(\mu_{d_{k}}-\sigma_{d_{k}},\mu_{d_{k}}+\sigma_{d_{k}})\}$, where $\mu_{d_{k}}$ and $\sigma_{d_{k}}$ are the mean and standard deviation for each distance $d_{k}$ given by the CVAE. Then, an EDG algorithm determines a set of Cartesian coordinates $\{\mathbf{r}_{i}\}_{i=1}^{N_{v}}$ so that these bounds are fulfilled (see the Appendix for details). Often there exist multiple solutions for the same set of bounds. As the bounds are generally tight, the solutions are very similar. Therefore, we only generate one set of coordinates per set of bounds. Together with the corresponding chemical elements $\{\epsilon_{i}\}_{i=1}^{N_{v}}$, we obtain a conformation $\mathbf{x}$.

4 Calculation of Molecular Properties

Related Works

The standard approach for generating molecular conformations is to start with one, and make small changes to it over time, e.g., by using MCMC or MD. These methods are considered the gold standard for sampling equilibrium states, but they are computationally expensive, especially if the molecule is large and the Hamiltonian is based on quantum-mechanical principles (Shim & MacKerell 2011; Ballard et al. 2015; De Vivo et al. 2016).

A much faster but more approximate approach for conformation generation is EDG (Havel 2002; Blaney & Dixon 2007; Lagorce et al. 2009; Riniker & Landrum 2015). Lower and upper distance bounds for pairs of atoms in a molecule are fixed values based on ideal bond lengths, bond angles, and torsional angles. These values are often extracted from crystal structure databases (Allen 2002). These methods aim to produce a low-energy conformation, not to generate unbiased samples from the underlying distribution at a certain temperature.

There exist several machine learning approaches as well, however, they are mostly tailored towards studying protein dynamics. For example, Noé et al. 2019 trained Boltzmann generators on the energy function of proteins to provide unbiased, one-shot samples from their equilibrium states. This is achieved by training an invertible neural network to learn a coordinate transformation from a system’s configurations to a latent space representation. Further, Lemke & Peter 2019 proposed a dimensionality reduction algorithm that is based on a neural network autoencoder in combination with a nonlinear distance metric to generate samples for protein structures. Both models learn protein-specific coordinate transformations that cannot be transferred to other molecules.

AlQuraishi 2019 introduced an end-to-end differentiable recurrent geometric network for protein structure learning based on amino acid sequences. Also, Ingraham et al. 2019 proposed a neural energy simulator model for protein structure that makes use of protein sequence information. Recently, Senior et al. 2020 significantly advanced the field of protein-structure prediction with a new model called AlphaFold. In contrast to amino acid sequences, molecular graphs are, in general, not linear but highly branched and often contain cycles. This makes these approaches unsuitable for general molecules.

Finally, Mansimov et al. 2019 presented a conditional deep generative graph neural network to generate molecular conformations given a molecular graph. Their goal is to predict the most likely conformation and not a distribution over conformations. Instead of encoding molecular environments in atomic distances, they work directly in Cartesian coordinates. As a result, the generated conformations showed significant structural differences compared to the ground-truth and required refinement through a force field, which is often employed in MD simulations.

We argue that our model has several advantages over the approaches reviewed above:

It is a fast alternative to resource-intensive approaches based on MCMC or MD.

Our principled representation based on pair-wise distances does not restrict our approach to any particular molecular structure.

Our model is, in principle, transferable to unseen molecules.

The Conf17 Benchmark

The Conf17 benchmark is the first benchmark for molecular conformation sampling. Datasets such as the one published by Kanal et al. 2018 only include conformers, i.e., the stable conformations of a molecule, and not a distribution over conformations. It is based on the ISO17 dataset (Schütt et al. 2017a) which consists of conformations of various molecules with the atomic composition C7H10O2 drawn from the QM9 dataset (Ramakrishnan et al. 2014). These conformations were generated by ab initio molecular dynamics simulations at 500 Kelvin. From the ISO17 dataset, 430692 valid molecular graph-conformation pairs could be extracted and 197 unique molecular graphs could be identified. We split the dataset into training and test sets such that no molecular graph in the training set can be found in the test or vice versa. Training and test splits consist of 176 and 30 unique molecular graphs, respectively (see Appendix A for details).

In Fig. 4, A, the structural formulae of a random selection of molecules from this benchmark are shown. Most molecules feature highly-strained, complex 3D structures such as rings which are typical of drug-like molecules. It is thus the structural complexity of the molecules, not their number of degrees of freedom, that makes this benchmark challenging. In Fig. 4, B–D, the frequency of distances (in Å) in the conformations are shown for each edge type. It can be seen that the marginal distributions of the edge distances are multimodal and highly context-dependent.

Experiments

We assess the performance of our method, named Graph Distance Geometry (GraphDG), by comparing it with two state-of-the-art methods for molecular conformation generation: RDKit (Riniker & Landrum 2015), a classical EDG approach, and DL4Chem (Mansimov et al. 2019), a machine learning approach. We trained GraphDG and DL4Chem on three different training and test splits of the Conf17 benchmark using Adam (Kingma & Ba 2014). We generated $100$ conformations with each method for molecular graphs in a test set.

We assessed the accuracy of the distance distributions of RDKit, DL4Chem, and GraphDG by calculating the maximum mean discrepancy (MMD) (Gretton et al. 2012) to the ground-truth distribution. In particular, we compute the MMD using a Gaussian kernel, where we set the standard deviation to be the median distance between distances $\mathbf{d}$ in the aggregate sample. For this, we determined the distances in the conformations from the ground-truth and those generated by RDKit, DL4Chem, and GraphDG. For each train-test split and each $G$ in a test set, we compute the MMD of the joint distribution of distances between C and O atoms $p(\{d_{k}\}|G)$ (H atoms are usually ignored), the MMDs of pair-wise distances $p(d_{i},d_{j}|G)$, and the MMDs between the marginals of individual distances $p(d_{i}|G)$. We aggregate the results of three train-test splits, and, finally, compute the median MMDs and average rankings. The results are summarized in Table 1. It can be seen that the samples from GraphDG are significantly closer to the ground-truth distribution than the other methods. RDKit is slightly worse than GraphDG while DL4Chem seems to struggle with the complexity of the molecules and the small number of graphs in the training set.

In Fig. 5, we showcase the accuracy of our model by plotting the marginal distributions $p(d_{i}|G)$ for distances between C and O atoms, given a molecular graph from a test set. It can be seen that RDKit consistently underestimates the marginal variances. This is because this method aims to predict the most stable conformation, i.e., the distribution’s mode. In contrast, DL4Chem often fails to predict the correct mean. For this molecule, GraphDG is the most accurate, predicting the right mean and variance in most cases. Additional figures can be found in the Appendix, where we also show plots for the marginal distributions $p(d_{i},d_{j}|G)$.

2 Generation of Conformations

We passed the distances from our generative model to an EDG algorithm to obtain conformations. For 99.9% of the sets of distances, all triangle inequalities held. For 83% of the molecular graphs, the algorithm succeeded which is 7 pp higher than the success rate we observed for RDKit. For each molecular graph in a test set, we generated 50 conformations with each method. This took DL4Chem, RDKit, and GraphDG on average around hundreds of milliseconds per molecule. All simulations were carried out on a computer equipped with an i7-3820 CPU and a GeForce GTX 1080 Ti GPU. In contrast, a single conformation in the ISO17 dataset takes around a minute to compute.

To assess the approximations made in the IS scheme, we studied the overlap between $p(\mathbf{d}|z,\mathcal{G})$ for a given $\mathcal{G}$ and different samples of $z$. We found experimentally that for 50 samples the overlap between the distributions is small. This finding can be explained by the high dimensionality of $\mathbf{d}$ which is on average $\approx 60$.

In Fig. 6, an overlay of these conformations of six molecules generated by the different methods is shown. It can be seen that RDKit’s conformations show too little variance, while DL4Chem’s structures are mostly invalid, which is due in part to its failure to predict the correct interatomic angles. Our method slightly overestimates the structural variance (see, for example, Fig. 6, top row, second column), but produces conformations that are the closest to the ground-truth.

3 Calculation of Molecular Properties

We estimate expected molecular properties for molecular graphs from the test set with $N=50$ conformational samples each. Due to their poor quality, we could not compute properties $\mathcal{O}(\mathbf{x})$, including the energy $E(\mathbf{x})$, for conformations generated with DL4Chem, and thus, this method is excluded from this analysis. In Table 2, it can be seen that RDKit and GraphDG perform similarly well (computational details can be found in the Appendix). However, both methods are still highly inaccurate for $E_{\text{elec}}$ (in practice, an accuracy of less than 5 kJ/mol is required). Close inspection of the conformations shows that, even though GraphDG predicts the most accurate distances overall, the variances of certain strongly constrained distances (e.g., triple bonds) are overestimated so that the energies of the conformations increase drastically.

Limitations

The first limitation of this work is that the CVAE can sample invalid sets of distances for which there exists no 3D structure. Second, the Conf17 benchmark covers only a small portion of chemical space. Finally, a large set of auxiliary edges would be required to capture long-range correlations (e.g., in proteins). Future work will address these points.

Conclusions

We presented GraphDG, a transferable, generative model that allows sampling from a distribution over molecular conformations. We developed a principled learning representation of conformations that is based on distances between atoms. Then, we proposed a challenging benchmark for comparing molecular conformation generators. With this benchmark, we show experimentally that conformations generated by GraphDG are closer to the ground-truth than those generated by other methods. Finally, we employ our model as a proposal distribution in an IS integration scheme to estimate molecular properties. While orbital energies and the dipole moments were predicted well, a larger and more diverse dataset will be necessary for meaningful estimates of electronic energies. Further, methods have to be devised to estimate how many conformations need to be generated to ensure all important conformations have been sampled. Finally, our model could be trained on conformational distributions at different temperatures in a transfer learning-type setting.

Acknowledgments

We would like to thank the anonymous reviewers for their valuable feedback. We further thank Robert Perharz and Hannes Harbrecht for useful discussions and feedback. GNCS acknowledges funding through an Early Postdoc.Mobility fellowship by the Swiss National Science Foundation (P2EZP2_181616).

References

Appendix A Conf17 Benchmark

The ISO17 dataset (Schütt et al. 2017a) was processed in the following way. First, conformations in which could not be parsed by the tool XYZ2Mol (Jensen 2019) were discarded. Second, the molecular graphs were augmented by adding auxiliary edges for reasons described in the main text. This can lead to an over-specification of the system’s geometry, however, this did not pose a problem in our experiments.

A.2 Input Features

Below we list the node and edge features in the Conf17 benchmark.

Appendix B Model Architecture

The source code of the model (including pre-processing scripts) is available online https://github.com/gncs/graphdg. In Table 5, the model architecture are summarized. After each hidden layer, a ReLU non-linearity is used. In Table 6, all hyperparameters are listed.

Appendix C Computational Details

All quantum-chemical calculations were carried out with the PySCF program package (version 1.5) (Sun et al. 2018) employing the exchange-correlation density functional PBE (Perdew et al. 1996), and the def2-SVP (Weigend & Ahlrichs 2005; Weigend 2006) basis set.

Conformations generated by DL4Chem did not succeed as some atoms were too close to each other. Self-consistent field algorithms in quantum-chemical software such as PySCF do not converge for such molecular structures.

With quantum-chemical methods, we calculate several properties that concern the states of the electrons in the conformation. These are the total electronic energy $E_{\text{elec}}$, the energy of the electron in the highest occupied molecular orbital (HOMO in eV) $\epsilon_{\text{HOMO}}$, the energy of the lowest unoccupied molecular orbital (LUMO in eV) $\epsilon_{\text{LUMO}}$, and the norm of the dipole moment $\mu$ (in debye).

C.2 Euclidean Distance Geometry

We refer the reader to Havel 2002 for theory on EDG, algorithms, and chemical applications. In summary, the EDG procedure consists of the following three steps:

Bound smoothing: extrapolating a complete set of lower and upper limits on all the distances from the sparse set of lower and upper bounds.

Embedding: choosing a random distance matrix from within these limits, and computing coordinates that are a certain best-fit to the distances.

Optimization: optimizing these coordinates versus an error function which measures the total violation of the distance (and chirality) constraints.

We use the EDG implementation found in RDKit (Riniker & Landrum 2015) with default settings.

Appendix D Generation of Conformations

Fig. 7 shows an overlay of 50 conformations from the ground-truth, RDKit, DL4Chem, and GraphDG based on two random molecular graphs from a test set.

Appendix E Distributions over Distances

We show the marginal distributions $p(d_{k}|G)$ and $p(d_{i},d_{j}|G)$ of ground-truth and predicted distances (in Å) for additional molecules from a test set.