Symmetry-Aware Actor-Critic for 3D Molecular Design

Introduction

The search for molecular structures with desirable properties is a challenging task with important applications in de novo drug design and materials discovery (Schneider et al. 2019). There exist a plethora of machine learning approaches to accelerate this search, including generative models based on variational autoencoders (VAEs) (Gómez-Bombarelli et al. 2018), recurrent neural networks (RNNs) (Segler et al. 2018), and generative adversarial networks (GANs) (De Cao & Kipf 2018). However, the reliance on a sufficiently large dataset for exploring unknown regions of chemical space is a severe limitation of such supervised models. Recent RL-based methods (e.g., Olivecrona et al. 2017, Jørgensen et al. 2019, Simm et al. 2020) mitigate the need for an existing dataset of molecules as they only require access to a reward function.

Most approaches rely on graph representations of molecules, where atoms and bonds are represented by nodes and edges, respectively. This is a strongly simplified model designed for the description of single organic molecules. It is unsuitable for encoding metals and molecular clusters as it lacks information about the relative position of atoms in 3D space. Further, geometric constraints on the design process cannot be included, e.g. those given by the active site of an enzyme. A more general representation closer to the physical system is one in which a molecule is described by its atoms’ positions in Cartesian coordinates. However, it would be very inefficient to naively learn a model based on this representation. That is because molecular properties such as the energy are invariant (i.e. unchanged) under symmetry operations like translation or rotation of all atomic positions. A model without the right inductive bias would thus have to learn those symmetries from scratch.

In this work, we develop a novel RL approach for designing molecules in Cartesian coordinates that explicitly encodes these symmetry operations. The agent builds molecules by consecutively placing atoms such that if the generated structure is rotated or translated, the agent’s action is rotated and translated accordingly; this way, the reward remains the same (see Fig. 1 (a)). We achieve this through a rotationally covariant state representation based on spherical harmonics, which we integrate into a novel actor-critic network architecture with an auto-regressive policy that maintains the desired covariance. Building in this inductive bias enables us to generate molecular structures with more complex coordination geometry than the class of molecules that were attainable with previous approaches. Finally, we perform experiments on several 3D molecular design tasks, where we find that our approach significantly improves the generalization capabilities of the RL agent and the quality of the generated molecules.

In summary, our contributions are as follows:

we propose the first approach for 3D molecular design that exploits symmetries of the design process by leveraging a rotationally covariant state representation;

we integrate this state representation into an actor-critic neural network architecture with a rotationally covariant auto-regressive policy, where the orientation of the atoms to be placed is modeled through a flexible distribution based on spherical harmonics;

we demonstrate the benefits of our approach on several 3D molecular design tasks, including a newly proposed task that showcases the generalization capabilities of our agent.

Background

The reward function $r(s_{t},a_{t})=-\Delta E(s_{t},a_{t})$ is given by the negative energy difference between the resulting structure described by $\mathcal{C}_{t+1}$, and the sum of energies of the current structure $\mathcal{C}_{t}$ and a new atom of element $e_{t}$ placed at the origin, i.e. $\Delta E(s_{t},a_{t})=E(\mathcal{C}_{t+1})-\left[E(\mathcal{C}_{t})+E(\{(e,\bm{0})\})\right]$. Intuitively, the reward encourages the agent to build stable, low-energy structures. We evaluate the energy using the fast semi-empirical Parametrized Method 6 (PM6) (Stewart 2007) as implemented in Sparrow (Husch et al. 2018; Bosia et al. 2020); see Appendix A for details.

2 Rotationally Covariant Neural Networks

Covariant Policy for Molecular Design

An efficient RL agent needs to exploit the symmetries of the molecular design process. Therefore, we require a policy $\pi(a|s)$ with actions $a=(e,x)$ that is covariant under translation and rotation with respect to the position $x$, i.e., $x$ should rotate (or translate) accordingly if the atoms on the canvas $\mathcal{C}$ are rotated (or translated). In contrast, the policy needs to be invariant to the element $e$, i.e. the chosen element remains unchanged under such transformations (see Fig. 1 (a)). Since learning such a policy is difficult when working directly in global Cartesian coordinates, we instead follow Simm et al. 2020 and use an action representation that is local with respect to an already placed focal atom. If the next atom is placed relative to the focal atom, covariance under translation of $x$ is automatically achieved and only the rotational covariance remains to be dealt with.

A novel actor-critic neural network architecture that implements this policy is illustrated in Fig. 3. In the following, we discuss its state embedding, actor, and critic networks in more detail.

2 Actor

Focal Atom and Element The distribution $p(f|s)$ over the focal atom $f$ is modeled as categorical, $f\sim\text{Cat}(f;h_{f})$, where $h_{f}$ are the logits for each atom in $\mathcal{C}$ predicted by a multi-layer perceptron (MLP). Likewise, the distribution over the element $e$ is given by $p(e|f,s)=\text{Cat}(e;h_{e})$ with $h_{e}=\text{MLP}_{e}(s^{\text{inv}}_{f})$, where $s^{\text{inv}}_{f}$ is the invariant representation for the focal atom. Since the number of possible focal atoms $f$ increases and the set of available elements $e$ decreases during a rollout, we mask out invalid focal atoms $f\notin\{1,\dots,|\mathcal{C}_{t}|\}$ and elements $e\notin\mathcal{B}_{t}$ by setting their probabilities to zero and re-normalizing the categorical distributions. The agent does not make use of chemical concepts like bond connectivity to aid the choice of the focal atom.

3 Critic

The critic needs to compute a value $V$ for the state $s$ that is invariant under translation and rotation. Given $s^{\text{inv}}$, we apply a permutation-invariant set encoding (Zaheer et al. 2017) of the atoms, i.e.

Related Work

Reinforcement Learning for Molecular Design There exists a large variety of RL-based approaches for molecular design using either string- or graph-based representations of molecules (Olivecrona et al. 2017; Guimaraes et al. 2018; Putin et al. 2018; Neil et al. 2018; Popova et al. 2018; You et al. 2018; Zhou et al. 2019). However, the choice of representation limits the molecules that can be generated to a (small) region of chemical space for which the representation is applicable, i.e., single organic molecules. Such representations also prohibit the use of reward functions based on quantum-mechanical properties; instead, heuristics are often used. Lastly, geometric constraints on the design process cannot be imposed as the representation does not include any 3D information.

Molecular Design in Cartesian Coordinates Another downside of string- and graph-based approaches is their neglect of information encoded in the interatomic distances. In light of this, Gebauer et al. 2018; Gebauer et al. 2019 proposed a supervised generative neural network for sequentially placing atoms in Cartesian coordinates. While the model respects local symmetries by construction, atoms are placed on a 3D grid. Similar to other supervised approaches, one further requires a dataset that covers the particular class of molecules to be generated. Hammer and coworkers (Jørgensen et al. 2019; Meldgaard et al. 2020) employed a Deep Q-Network (Mnih et al. 2015) to build planar compounds and crystalline surfaces by placing atoms on a grid. Recently, Simm et al. 2020 presented an RL formulation for molecular design in continuous 3D space. The agent models the position of the next atom to be placed in internal coordinates—i.e. the distance, angle, and dihedral angle with respect to already existing atoms—which are invariant under translation and rotation. By mapping from internal to Cartesian coordinates, they then obtain a policy that is covariant under these symmetry operations. However, as shown in Fig 4, the angle and dihedral angle are only defined with respect to two reference points, which are chosen to be the two closest points to a focal atom. In highly symmetric states, e.g. as commonly encountered in materials, this representation fails to distinguish different configurations as one cannot uniquely select the two closest atoms as reference points anymore. In contrast, we do not rely on such reference points as the agent directly samples the orientation from a spherical distribution.

Covariant Neural Networks in Chemical Science Prior work employed rotationally covariant neural networks to predict translation- and rotation-invariant physical properties (Thomas et al. 2018; Kondor et al. 2018; Weiler et al. 2018; Anderson et al. 2019; Miller et al. 2020; Finzi et al. 2020; Fuchs et al. 2020), e.g. scalars such as the electronic energy. In contrast, we propose a translation-invariant and rotation-covariant neural network architecture for generating molecules. For a more general treatment of covariance (or equivariance) in RL, see van der Pol et al. 2020.

Experiments

We perform experiments to answer the following questions: (1) is the agent able to learn how to build highly symmetric molecules in Cartesian coordinates from scratch, (2) can we increase the validity, diversity, and stability of generated molecules, and (3) does our approach lead to improved generalization? We address (1) and (2) by evaluating the agent on a diverse range of tasks from the MolGym benchmark suite (Simm et al. 2020), and (3) on a newly proposed stochastic-bag task (see Section 5.1) where bags are sampled from a distribution over bags. In Appendix G, we show with an additional experiment that the agent can learn to place water molecules around a given solute to form a solvation shell.

We compare our approach (Covariant) against the RL agent proposed by Simm et al. 2020, which iteratively builds molecules on a 3D canvas by working in internal coordinates (Internal). As an additional baseline, we consider a classical, optimization-based agent (Opt) with access to a black-box function that yields the energy $E(\mathcal{C})$ and the atomic forces $F(\mathcal{C})$ for a given canvas. For the calculation of $E(\mathcal{C})$ and $F(\mathcal{C})$ we employ PM6; the same method as in the reward function. The agent constructs molecules by alternating between randomly placing an atom and optimizing the structure. Moreover, the agent applies several heuristics inspired by fundamental chemical concepts to guide the placement of atoms. To make the comparisons fair, we grant Opt a comparable computational budget in terms of the total number of energy computations. Finally, for some experiments, the best possible performance based on quantum-chemical calculations can be reported. See Appendices E and F for more details on the baselines and experimental settings, and Appendix H for an additional runtime comparison between the agents.

2 Results

Results are listed in Table 1. We observe that Covariant significantly outperforms the other agents on most experiments both in terms of validity and diversity. The difference is particularly large for the more challenging stochastic-bag tasks, where Covariant does similarly well as on the single-bag experiments. This finding is confirmed in Fig. 6 (a), showing that the exact stoichiometry does not need to be known a priori for the agent to build valid molecules. Moreover, the structures generated by Covariant are overall slightly more stable compared to Internal. In contrast, Opt often fails to build valid structures. Inspection of the generated structures reveals that for larger bags the agent tends to build multi-molecular clusters, which are considered invalid in this experiment. The stability for Opt is omitted as all of its valid structures are stable by definition.

Compared to graph-based approaches (e.g., Jin et al. 2017; Bradshaw et al. 2019a; Li et al. 2018b; Li et al. 2018a; Liu et al. 2018; De Cao & Kipf 2018; Bradshaw et al. 2019b), the average validity and diversity achieved by Covariant are still relatively low. This can partly be explained by the fact that state-of-the-art graph-based approaches have the strict rules of chemical bonding in organic molecules encoded into their models. But as a result, they are limited to generating single organic molecules and cannot build molecules for which these rules do not apply (e.g., hypervalent iodine compounds such as $\text{IF}{\vphantom{\text{X}}}_{\smash[t]{\text{5}}}$). In terms of stability, the supervised generative model by Gebauer et al. 2019 reported an average RMSD of approximately $0.25$ Å. While their approach and the considered molecules are significantly different from ours, this suggests that the generated structures are more stable compared to Covariant. Nonetheless, the RL approach presented in this work remains particularly attractive if no dataset exists on which such a supervised model can be trained.

Conclusion

We proposed a novel covariant actor-critic architecture based on spherical harmonics for designing highly symmetric molecules in 3D. We showed empirically that exploiting symmetries of the molecular design process improves the quality of the generated molecules and leads to better generalization. In future work, we aim to employ more accurate quantum-chemical methods (e.g., density functional theory) required for building transition metal complexes or structures in which weak intermolecular interactions are important. For that, however, the sample-efficiency of our agent needs to be improved. Finally, we aim to explore reward functions specifically tailored towards drug design.

Acknowledgements

We would like to thank Austin Tripp and Kris Jensen for useful discussions and feedback. Robert Pinsler receives funding from iCASE grant #1950384 with support from Nokia. This work has been performed using resources provided by the Cambridge Tier-2 system operated by the University of Cambridge Research Computing Service funded by EPSRC Tier-2 capital grant EP/P020259/1.

References

Appendix A Reward Calculation

In the reward function, the energy $E$ has to be computed using quantum-chemical methods. For that, we use the fast semi-empirical Parametrized Method 6 (PM6) (Stewart 2007). In particular, we use the implementation in the software package Sparrow (Husch et al. 2018; Bosia et al. 2020). For each calculation, a molecular charge of zero and the lowest possible spin multiplicity are chosen. All calculations are spin-unrestricted.

Limitations of semi-empirical methods are highlighted in, for example, recent work by Husch & Reiher 2018. More accurate methods such as approximate density functionals need to be employed especially for systems containing transition metals.

Further, we enforce that atoms are not placed too close ($<$ 0.6 Å) nor too far away from each other ($>$ 2.0 Å). If the agent places an atom outside these boundaries, the minimum reward of $-0.6$ is awarded and the episode terminates. Further, the environment encourages the agents to build single molecular structures by terminating the episode and return a reward of $-0.6$ if elements forming stable bimolecular compounds (e.g, $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}$) are placed too far away from other atoms on the canvas.

Appendix B Spherical Harmonics

The spherical harmonics are normalized such that:

Appendix C Calculation of Invariant Features

One can obtain scalar invariats from the covariant features $\hat{f}$ (Anderson et al. 2019):

Appendix D Probability Distribution for Orientation

Appendix E Baselines

Below, we detail the algorithm of the Opt agent. At the beginning of each experiment, the agent is given a canvas $\mathcal{C}_{0}$, a bag $\mathcal{B}_{0}$, and a black-box function that can compute the energy $E(\mathcal{C})$ and the atomic forces $F(\mathcal{C})$ for a given canvas. We assume a total charge of zero and a low-spin configuration. At the end of each experiment, we compute the total reward obtained for the final structure on canvas $\mathcal{C}_{T}$ and report the total number of energy and gradient computations.

If the canvas $\mathcal{C}_{t}$ is not empty, randomly choose a focal atom $f$ from the list of available atoms on the canvas. An atom is considered available if its number of neighbors is less than a predefined number that depends on its element (e.g., one for hydrogen and four for carbon). Two atoms on the canvas are neighbors if their Euclidean distance is below 1.5 Å. If there are no available atoms on the canvas, a focal atom is randomly chosen from the list of atoms on the canvas.

Randomly choose an element $e_{t}$ from the bag $\mathcal{B}_{t}$.

Randomly place the atom $a_{t}=(e_{t},x_{t})$ on a sphere with radial distance $d=1.1$ Å around $x_{f}$ to obtain $\mathcal{C}_{t+1,\text{raw}}$. If the canvas is empty, place the atom at the origin.

Optimize only the position of $a_{t}$ using $F$ to obtain $\mathcal{C}_{t+1,\text{opt}}$.

Compute the energy difference $\Delta E(t)=E(\mathcal{C}_{t+1,\text{opt}})-\left[E(\mathcal{C}_{t})+E(\{e_{t},\bm{0}\})\right]$.

If $\Delta E(t)>0$, return $e_{t}$ to the bag and go back to step 1.

Optimize canvas $\mathcal{C}_{t+1,\text{opt}}$ using $F$ to obtain $\mathcal{C}_{t+1}$.

If the bag is not empty, go back to step 1.

In the experiments, the different agents need to be given a comparable computational budget to ensure a meaningful comparison of their performance. This is difficult as they use different computational resources: Opt runs on a CPU whereas Internal and Covariant perform many of their computations on a GPU. However, we found experimentally that the quantum-chemical calculations are the most computationally expensive ones. These calculations are performed in the same way for all approaches. Therefore, we believe that by granting each approach the same number of PM6 calculations we achieve a fair comparison.

E.2 Optimal Return

The optimal return for the single-bag tasks was derived in the following way. First, we obtained molecular structures for the complexes $\text{SOF}{\vphantom{\text{X}}}_{\smash[t]{\text{4}}}$, $\text{IF}{\vphantom{\text{X}}}_{\smash[t]{\text{5}}}$, $\text{SF}{\vphantom{\text{X}}}_{\smash[t]{\text{6}}}$, and $\text{SOF}{\vphantom{\text{X}}}_{\smash[t]{\text{6}}}$. Subsequently, we performed a structure optimization using the PM6 method. Since the undiscounted return is path-independent, we determined the return $R(s)$ by computing the total interaction energy in the canvas $\mathcal{C}$, i.e.

Appendix F Experimental Details

Experiments were run on an Intel Xeon E5-2650 v4 2.2GHz 12-core processor (96GiB RAM) and an Nvidia P100 GPU (16GiB). Our agent is implemented in the deep learning framework PyTorch (Paszke et al. 2019). Data analysis was performed with the Python libraries matplotlib (Hunter 2007) and pandas (McKinney 2010).

F.2 Implementation Details

The model architecture is summarized in Table 2, where the dimensions of $s^{\text{inv}}$ and $s^{\text{inv}}_{f,e}$ are $d^{\text{inv}}=(L_{\text{max}}+2)\cdot\tau\cdot 2$ and $d^{\text{inv}}_{f,e}=(L_{\text{max}}+2)\cdot\tau_{e}\cdot 2$, respectively. If possible, we made similar architectural choices as Simm et al. 2020, e.g. regarding the number of hidden units/layers, activation functions, and initialization schemes. We initialize the biases of each network with $0$ and each weight matrix as a (semi-)orthogonal matrix. After each hidden layer, a ReLU non-linearity is employed. As explained in the main text, both $\text{MLP}_{f}$ and $\text{MLP}_{e}$ use a masked softmax activation function to guarantee that only valid actions are chosen. To model the distance $d$, we employ a Gaussian mixture model consisting of $M=3$ Gaussians. As we treat the standard deviations $\{\sigma_{m}\}_{m=1}^{3}$ as global parameters, the MDN has 6 outputs. Further, we rescale the means $\mu_{m}\in$ to $\mu_{m}\in[d_{\text{min}},d_{\text{max}}]$. If the sampled distance is negative, we clip the value at $0.001$.

Hyperparameters for Cormorant are listed in Table 3. In our experiments, we found it important to use multiple filters $\tau_{e}$ per element (e.g. $4$) and to set $L_{\text{max}}=4$. This gives the model enough flexibility to represent complex spherical distributions while remaining computationally tractable. For more details on Cormorant, see the original work (Anderson et al. 2019). Further hyperparameters used in our experiments are in Table 4. PPO is known to be relatively robust with respect to the choice of hyperparameters, and we found the default values to be sufficient in most cases. Within the actor, the scaling parameter $\beta$ is important to avoid that the spherical distribution approaches a delta distribution. Note that values of $\beta$ can vary significantly across experiments and might require some tuning. Lastly, the right number of samples $S$ for the global mode estimation of the spherical distribution generally depends on the shape of the distribution. In particular, we would expect that more samples are required as the distribution becomes more peaked. Since we avoid pathological behaviors by scaling the distribution with $\beta$, we found $S=1024$ to be sufficient for all our experiments.

Appendix G Additional Results

Next, we assess the ability of Covariant to generate solvation clusters—a type of molecular structure that cannot be built with graph-based approaches. Following Simm et al. 2020, we task the agent to place $5$ water molecules around a formaldehyde molecule that is already on the canvas at the beginning of each episode. In addition, the reward function is augmented with a penalty term for placing atoms far away from the center, i.e. $r(s_{t},a_{t})=-\Delta_{E}-\rho\|x\|_{2}$, where $\rho$ is a hyper-parameter that is set to $0.01$ (see Simm et al. 2020 for details). Therefore, the agent needs to place the water molecules such that hydrogen bonds can be formed between water molecules and between water molecules and the solute.

From Fig. 10, it can be seen that Covariant can solve this task by constructing stable $\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{2}}}\text{O}$ molecules and placing them in the vicinity of the solute. From visual inspection of the generated structures, it can be observed that in many cases Covariant arranges the molecules such that intermolecular bonds can be formed. However, it should be noted that the quantum-chemical method used in the reward function is not very well suited for modeling these interactions. Finally, Fig. 10 shows that while Internal learns faster at the beginning of training, Covariant is slightly outperforming Internal towards the end.

Appendix H Runtime Evaluation

We compared the runtimes between Opt, Internal, and Covariant. For instance, for the single-bag task with the bag $\text{C}{\vphantom{\text{X}}}_{\smash[t]{\text{3}}}\text{H}{\vphantom{\text{X}}}_{\smash[t]{\text{5}}}\text{NO}{\vphantom{\text{X}}}_{\smash[t]{\text{3}}}$, $T=240$ steps of the last rollout took Covariant and Internal on average 12 and 11 seconds (s), respectively. The final offline evaluation took the agents on average 4 and 1s, respectively. This speed difference is mainly due to the relatively slow rejection sampling procedure in Covariant. Each iteration, policy optimization took on average 2 and 6s for the agents Covariant and Internal, respectively. In this case, Internal is slower than Covariant as it performed around twice as many epochs during optimization due to early stopping. Since there is no training for Opt, this agent was overall faster than the others. Further, we note that the largest fraction of time was spent on the quantum-chemical calculations which are the same for all agents. The time the quantum-chemical calculation takes to converge depends not only on the size but also on the geometry of the input structure. The entire experiment took Covariant approximately 4 hours, Internal 5 hours, and Opt 3 hours.