Year: 2021

Abstract:
We developed Distilled Graph Attention Policy Networks (DGAPNs), a reinforcement learning model to generate novel graph-structured chemical representations that optimize user-defined objectives by efficiently navigating a physically constrained domain. The framework is examined on the task of generating molecules that are designed to bind, noncovalently, to functional sites of SARS-CoV-2 proteins. We present a spatial Graph Attention (sGAT) mechanism that leverages self-attention over both node and edge attributes as well as encoding the spatial structure — this capability is of considerable interest in synthetic biology and drug discovery. An attentional policy network is introduced to learn the decision rules for a dynamic, fragment-based chemical environment, and state-of-the-art policy gradient techniques are employed to train the network with stability. Exploration is driven by the stochasticity of the action space design and the innovation reward bonuses learned and proposed by random network distillation. In experiments, our framework achieved outstanding results compared to state-of-the-art algorithms, while reducing the complexity of paths to chemical synthesis.

heading: 1 INTRODUCTION
text: This work aims to address the challenge of establishing an automated process for the design of objects with connected components, such as molecules, that optimize specific properties. Achieving this goal is particularly desirable in drug development and materials science, where manual discovery remains a time-consuming and expensive process (Hughes et al., 2011; Schneider et al., 2020). However, there are two major difficulties that have long impeded rapid progress. Firstly, the chemical space is discrete and massive (Polishchuk et al., 2013), presenting a complicated environment for an Artificial Intelligence (AI) approach to efficiently and effectively explore. Secondly, it is not trivial to compress such connected objects into feature representations that preserve most of the information, while also being highly computable for Deep Learning (DL) methods to exploit.
We introduce Distilled Graph Attention Policy Networks (DGAPN), a framework that advances prior work in addressing both of these challenges. We present a Reinforcement Learning (RL) architecture that is efficiently encouraged to take innovative actions with an environment that is able to construct a dynamic and chemically valid fragment-based action space. We also propose a hybrid Graph Neural Network (GNN) that comprehensively encodes graph objects’ attributes and spatial structures in addition to adjacency structures. The following paragraphs discuss how we addressed limitations of prior work and its relevance to antiviral drug discovery. For more descriptions of key prior methodologies that we used as benchmarks in this paper, see Section 5.
Graph Representation Learning Despite their spatial efficiency, string representation of molecules acquired by the simplified molecular-input line-entry system (SMILES) (Weininger, 1988) suffers from significant information loss and poor robustness (Liu et al., 2017). Graph representations have become predominant and preferable for their ability to efficiently encode an object’s scaffold structure and attributes. Graph representations are particularly ideal for RL since intermediate representations can be decoded and evaluated for reward assignments. While GNNs such as Graph Convolutional Networks (GCNs) (Kipf & Welling, 2016) and Graph Attention Networks (GATs) (Veličković et al., 2017) have demonstrated impressive performance on many DL tasks, further exploitation into richer information contained in graph-structured data is needed to faithfully
represent the complexity of chemical space (Morris et al., 2019; Wang et al., 2019; Chen et al., 2020). In this work, we made improvements to previous studies on attributes encoding and structural encoding. For structural encoding, previous studies have covered adjacency distance encoding (Li et al., 2020), spatial cutoff (Pei et al., 2020) and coordinates encoding (Schütt et al., 2017; Danel et al., 2020). Our work presents an alternative approach to spatial structure encoding similar to Gilmer et al. (2017) which do not rely on node coordinates, but different in embedding and updating scheme. Distinct from Danel et al. (2020) and Chen & Chen (2021), we extended attentional embedding to be edge-featured, while still node-centric for message passing efficiency.
Reinforcement Learning A variety of graph generative models have been used in prior work, predominantly Variational Autoencoders (VAEs) (Simonovsky & Komodakis, 2018; Samanta et al., 2020; Liu et al., 2018; Ma et al., 2018; Jin et al., 2018) and Generative Adversarial Networks (GANs) (De Cao & Kipf, 2018). While some of these have a recurrent structure (Li et al., 2018; You et al., 2018b), RL and other search algorithms that interact dynamically with the environment excel in sequential generation due to their ability to resist overfitting on training data. Both policy learning (You et al., 2018a) and value function learning (Zhou et al., 2019) have been adopted for molecule generation: however, they generate molecules node-by-node and edge-by-edge. In comparison, an action space consisting of molecular fragments, i.e., a collection of chemically valid components and realizable synthesis paths, is favorable since different atom types and bonds are defined by the local molecular environment. Fragment-by-fragment sequential generation has been used in VAE (Jin et al., 2018) and search algorithms (Jin et al., 2020; Xie et al., 2021), but has not been utilized in a deep graph RL framework. In this work, we designed our environment with the Chemically Reasonable Mutations (CReM) (Polishchuk, 2020) library to realize a valid fragment-based action space. Furthermore, we enhanced exploration by employing a simple and efficient technique, adapting Random Network Distillation (RND) (Burda et al., 2018) to GNNs and proposing surrogate innovation rewards for intermediate states.
Antiviral Drug Discovery — A Timely Challenge The severity of the COVID-19 pandemic highlighted the major role of computational workflows to characterize the viral machinery and identify druggable targets for the rapid development of novel antivirals. Particularly, the synergistic use of DL methods and structural knowledge via molecular docking simulations is at the cutting edge of molecular biology — consolidating such integrative protocols to accelerate drug discovery is of paramount importance (Yang et al., 2021; Jeon & Kim, 2020; Thomas et al., 2021). Here we experimentally examined our architecture on the task of discovering novel inhibitors targeting the SARS-CoV-2 non-structural protein endoribonuclease (NSP15), which is critical for viral evasion of host defense systems (Pillon et al., 2021). Structural information about the putative protein-ligand complexes was integrated into this framework with AutoDock-GPU (Santos-Martins et al., 2021), which leverages the GPU resources from leadership-class computing facilities, including the Summit supercomputer, for high-throughput molecular docking (Legrand et al., 2020). We show that our results outperformed state-of-the-art generation models in finding molecules with high affinity to the target and reasonable synthetic accessibility.

heading: 2 PROBLEM FORMULATIONS
text: Our goal is to establish a set of decision rules to generate graph-structured data that maximizes compound objectives under certain constraints. Similar to prior formulations, the generating process is defined as a time homogeneous Markov Decision Process (MDP). We give a formal definition of this process in Appendix A. Under this setting, the action policies and state transition dynamics at step t can be factorized according to the Markov property:
P (at|s0, a0, s1, a1, . . . , st) = P (at|st) := π(at|st) (1) P (st+1|s0, a0, s1, a1, . . . , st, at) = P (st+1|st, at) := ρ(st+1|st, at) (2)
where {st, at}t are state-action sequences. A reward function r(s, a) is used to assess an action a taken at a given state s. The process terminates at an optional stopping time T and sT is then proposed as the final product of the current generating cycle. We aim to estimate the optimal policy π in terms of various objectives to be constructed later in the experiment section.

heading: 3 PROPOSED METHOD

heading: 3.1 ENVIRONMENT SETTINGS
text: In the case of molecular graphs, single-atom or single-bond additions are often not realizable by known biochemical reactions. Rather than employing abstract architectures such as GANs to suggest synthetic accessibility, we use the chemical library CReM (Polishchuk, 2020) to construct our environment such that all next possible molecules can be obtained by one step of interchanging chemical fragments with the current molecule. This explicit approach is considerably more reliable and interpretable compared to DL approaches. A detailed description of the CReM library can be found in Appendix B.1.
At each time step t, we use CReM to sample a set of valid molecules vt+1 as the candidates for the next state st+1 based on st. Under this setting, the transition dynamics are deterministic, the underlying set A of the action space can be defined as equal to S of the state space, and action at is induced by the direct selection of st+1. With an abuse of notation, we let r(st+1) := r(st, at).

heading: 3.2 SPATIAL GRAPH ATTENTION
text: We introduce a graph embedding mechanism called Spatial Graph Attention (sGAT) in an attempt to faithfully extract feature vectors ht ∈ Rdh representing graph-structured objects such as molecules. Two different types of information graphs constructed from a connected object are heterogeneous and thus handled differently in forward passes as described in the following sections. See Figure 1 for an overview.

heading: 3.2.1 ATTENTION ON ATTRIBUTION GRAPHS
text: The attribution graph of a molecule with n atoms and e bonds is given by the triple (A,N ,E), where A ∈ {0, 1}n×n is the node adjacency matrix, N is the node attribution matrix of dimension n× dn and E is the edge attribution matrix of dimension e× de. Each entry aij ofA is 1 if a bond exists between atom i and j, and 0 otherwise. Each row vector ni of N is a concatenation of the properties of atom i, including its atomic number, mass, etc., with the categorical properties being one-hot encoded. E is formed similar to N , but with bond attributes. We denote a row vector of E as eij if it corresponds to the bond between atom i and j.
We proceed to define a multi-head forward propagation that handles these rich graph information: let hnk ∈ R1×dhn denote a given representation for nk, heij ∈ R1×dhe denote a representation for eij , then the m-th head attention αmij from node j to node i (i 6= j) is given by
αmij = softmax j ( ⋃ k: aik=1 { σ([hniWn,m ‖ heikWe,m ‖ hnkWn,m] · attmT ) }) (3)
where softmax j is the softmax score of node j; ‖ is column concatenation; σ is some non-linear activation;Wn,m ∈ Rdhn×dwn ,We,m ∈ Rdhe×dwe are them-th head weight matrices for nodes and edges respectively; attm ∈ R1×(2dwn+dwe ) is the m-th head attention weight. The representations after a feed-forward operation are consequently given as follow:
h′ni = aggr1≤m≤nm σ  ∑
j: aij=1
αmij · hnj + hni Wn,m  (4)
h′eij = aggr1≤m≤nm { σ ([ hniWn,m ‖ heijWe,m ‖ hnjWn,m ] ·Wh,m )} (5)
where Wh,m ∈ R(2dwn+dwe )×dwe ; nm is the total number of attention heads and aggr denotes an aggregation method, most commonly mean , sum , or concat (Hamilton et al., 2017). We note that we have not found significant difference across these methods and have used mean for all aggregations in our experiments. In principle, a single-head operation on nodes is essentially graph convolution with the adjacency matrix Â = Ã + I where Ã is attention-regularized according to (3). This approach sufficiently embeds edge attributes while still being a node-centric convolution mechanism, for which efficient frameworks like Pytorch-Geometric (Fey & Lenssen, 2019) have been well established.

heading: 3.2.2 SPATIAL CONVOLUTION
text: In addition to attributions and logical adjacency, one might also wish to exploit the spatial structure of an graph object. In the case of molecular docking, spatial structure informs the molecular volume and the spatial distribution of interaction sites — shape and chemical complementarity to the receptor binding site is essential for an effective association.
Let G = ( dij −1)
i,j≤n be the inverse distance matrix where dij is the Euclidean distance between
node i and j for ∀i 6= j, and dii−1 := 0. G can then be seen as an adjacency matrix with weighted “edge”s indicating nodes’ spatial relations, and the forward propagation is thus given by
H ′′n = σ (( D̃− 1 2 G̃D̃− 1 2 + I ) HnWn ) (6)
where G̃ is optionally sparsified and attention-regularized from G to be described below; D̃ = diag1≤i≤n {∑n j=1 G̃ij } ; Hn is the row concatenation of {hni}1≤i≤n; Wn ∈ Rdhn×dwn is the
weight matrix. In reality,G inducesO(n) of convolution operations on each node and can drastically increase training time when the number of nodes is high. Therefore, one might want to derive G̃ by enforcing a cut-off around each node’s neighborhood (Pei et al., 2020), or preserving an O(n) number of largest entries in G and dropping out the rest. In our case, although the average number of nodes is low enough for the gather and scatter operations (GS) of Pytorch-Geometric to experience no noticeable difference in runtime as node degrees scale up (Fey & Lenssen, 2019), the latter approach of sparsification was still carried out because we have discovered that proper cutoffs improved the validation loss in our supervised learning experiments. If one perceives the relations between chemical properties and spatial information as more abstract, G should be regularized by attention as described in (3), in which case the spatial convolution is principally fully-connected graph attention with the Euclidean distance as a one-dimensional edge attribution.

heading: 3.3 GRAPH ATTENTION POLICY NETWORK
text: In this section we introduce Graph Attention Policy Network (GAPN) that is tailored to environments that possess a dynamic range of actions. Note that ρ(·|st, at) is a degenerate distribution for deterministic transition dynamics and the future trajectory τ ∼ p(st+1, st+2, . . . |st) is strictly equal in distribution to a ∼ π(at, at+1, . . . |st), hence simplified as the latter in the following sections. To learn the policy more efficiently, we let st and vt share a few mutual embedding layers, and further tested pre-training the first ng layers with supervised learning on an “expert” dataset in the later experiments. Layers inherited from pre-training are not updated during the training of RL. See Figure 2 for an overview of the architecture.

heading: 3.3.1 ACTION SELECTION
text: At each time step t, we sample the next state st+1 from a categorical distribution constructed by applying a retrieval-system-inspired attention mechanism (Vaswani et al., 2017):
st+1 ∼ OHC softmax  ⋃ g∈gt+1 {Lfinal(EQ(gt) ‖ EK(g)}  · vt+1 (7) where OHC{p1, . . . , pnv} is a one-hot categorical distribution with nv categories; gt, gt+1 are the embeddings for st and vt+1 acquired by the shared encoder; EQ, EK are two sGAT+MLP graph encoders with output feature dimension dk; Lfinal : Rb×2dk → Rb is the final feed-forward layer. Essentially, each candidate state is predicted a probability based on its ‘attention’ to the query state. The next state is then sampled categorically according to these probabilities.
There could be a number of ways to determine stopping time T . For instance, an intuitive approach would be to append st to vt+1 and terminate the process if st is selected as st+1. In our experiments, we simply pick T to be constant, i.e. we perform a fixed number of modifications for an input. This design encourages the process to not take meaningless long routes or get stuck in a cycle, and enables episodic docking evaluations in parallelization (further described in Section 3.5). Note that constant trajectory length is feasible because the maximum limit of time steps can be set significantly lower for fragment-based action space compared to node-by-node and edge-by-edge action spaces.

heading: 3.3.2 ACTOR-CRITIC ALGORITHM
text: For the purpose of obeying causal logic and reducing variance, advantages are predominantly used instead of raw rewards in policy iterations. Applying the actor-critic algorithm on discounted rewardto-go yields the following Q-function and advantage function:
Qπ(st, at) = Eπ [ T∑ t′=t γt ′−t · r(st′ , at′) ∣∣∣∣∣st, at ]
(8)
Aπ(st, at) = Q π(st, at)− Eπ [Qπ(st, at)|st] (9)
For a more detailed description of actor-critic algorithm in RL, see Grondman et al. (2012).

heading: 3.3.3 PROXIMAL POLICY OPTIMIZATION
text: We use Proximal Policy Optimization (PPO) (Schulman et al., 2017), a state-of-the-art policy gradient technique, to train our network. PPO holds a leash on policy updates whose necessity is elaborated in trust region policy optimization (TRPO) (Schulman et al., 2015), yet much simplified. It also enables multiple epochs of minibatch updates within one iteration. The objective function is
given as follow:
J∗(θ) = max θ Eµ,πoldθ [ T∑ t=1 min { rt(θ)A πoldθ (st, at), clip (rt(θ))A πoldθ (st, at) }] (10)
where rt(θ) = πnewθ (at ∣∣st)/πoldθ (at∣∣st) and clip (x) = min {max {1− , x} , 1 + }. During policy iterations, πnew is updated each epoch and πold is cloned from πnew each iteration.

heading: 3.4 EXPLORATION WITH RANDOM NETWORK DISTILLATION
text: We seek to employ a simple and efficient exploration technique that can be naturally incorporated into our architecture to enhance the curiosity of our policy. We perform Random Network Distillation (RND) (Burda et al., 2018) on graphs or pre-trained feature graphs to fulfill this need. Two random functions f̂ψ, f∗ that map input graphs to feature vectors in Rdr are initialized with neural networks, and f̂ψ is trained to match the output of f∗:
ψ∗ = arg min ψ Es′∼p̂next‖f̂ψ(s′)− f∗(s′)‖ (11)
where p̂next is the empirical distribution of all the previously selected next states, i.e. the states that have been explored. We record the running errors in a buffer and construct the surrogate innovation reward as:
ri(s ′) = clipη (( ‖f̂ψ(s′)− f∗(s′)‖ −mb )/√ vb ) (12)
where mb and vb are the first and second central moment inferred from the running buffer, clipη(x) = min {max {−η, x} , η}.

heading: 3.5 PARALLELIZATION AND SYNCHRONIZED EVALUATION
text: Interacting with the environment and obtaining rewards through external software programs are the two major performance bottlenecks in ours as well as RL in general. An advantage of our environment settings, as stated in Section 3.3.1, is that a constant trajectory length is feasible. Moreover, the costs for environmental interactions are about the same for different input states. To take advantage of this, we parallelize environments on CPU subprocesses and execute batched operations on one GPU process, which enables synchronized and sparse docking evaluations that reduces the number of calls to the docking program.

heading: 4 EXPERIMENTS

heading: 4.1 SETUP
text: Objectives We evaluated our model against five state-of-the-art models (detailed in Section 5) with the objective of discovering novel inhibitors targeting SARS-CoV-2 NSP15. Molecular docking scores are computed by docking programs that use the three-dimensional structure of the protein to predict the most stable bound conformations of the molecules of interest, targeting a pre-defined functional site. For more details on molecular docking and our GPU implementation of an automated docking tool used in the experiments, see Appendix B.2. In addition, we evaluated our model in the context of optimizing QED and penalized LogP values, two tasks commonly presented in machine learning literature for molecular design. The results for this can be found in Appendix C.
Dataset For the models/settings that do require a dataset, we used a set of SMILES IDs taken from more than six million compounds from the MCULE molecular library - a publicly available dataset of purchasable molecules (Kiss et al., 2012) and their docking scores for the NSP15 target.
Hyperparameters Based on a parameter sweep, we set number of GNN layers to be 3, MLP layers to be 3, with 3 of the GNN layers and 0 of the MLP layers shared between query and key. Number of layers in RND is set to 1; all numbers of hidden neurons 256; learning rate for actor 5−4, for critic 1−4, for RND 2−3; maximum time steps per episode 12, update time steps (i.e. batch
size) 200. Number of epochs per iteration and clipping parameter for PPO are 50 and 0.2. Output dimensions and clipping parameter η for RND are 8 and 5. For more details on hyperparameters and other settings, see our code in supplemental material.

heading: 4.2 RESULTS

heading: 4.2.1 SINGLE-OBJECTIVE OPTIMIZATION
text: The raw docking score is a negative value that represents higher estimated binding affinity when the score is lower. We use the negative docking score as the main reward rm and assign it to the final state sT as the single objective. For DGAPN, we also assign innovation reward to each intermediate state, and the total reward for a trajectory τ is given by r(τ ) = rm(sT ) + ι · ∑T t=1 ri(st) where ι is the relative important of innovation rewards, for which we chose 0.1 and incorporated them with a 100 episode delay and 1,000 episode cutoff. We sampled 1,000 molecules from each method and showed the evaluation results in Table 1. MARS, GCPN and MolDQN are evaluated by the latest molecules seen in training; REINVENT is evaluated in inference mode after training; JT-VAE is evaluated using Bayesian Optimization on the embedded spaces. We note that we have a separate approach for the evaluation of our model that is able to achieve a −9.3 mean and −13.89 best docking score (see the ablation study below), but in order to maintain consistency with GCPN and MolDQN, here we likewise only evaluated the latest molecules found in training.
In the table, validity is checked by examining atoms’ valency and consistency of bonds in aromatic rings, diversity is defined as the average pairwise Tanimoto distance between the Morgan fingerprints (Rogers & Hahn, 2010). Further, we provide additional summary metrics to help gain some perspectives of the generated molecules: QED (Bickerton et al., 2012) is an indicator of drug-likeness, SA (Ertl & Schuffenhauer, 2009) is the synthetic accessibility of molecules, FCD (Preuer et al., 2018) is the Fréchet ChemNet Distance. QED is better when the score is higher, while SA and FCD are the opposite. The details of these metrics can be found in Appendix D. On this task, DGAPN significantly outperformed state-of-the-art models in terms of the top scores and average score, obtaining a high statistical significance over the second best model (MolDQN) with a p-value of 1.67× 10−37 under Welch’s t-test (Welch, 1947). As anticipated, the molecules generated by fragment-based algorithms (JTVAE, MARS, DGAPN) have significantly better SAs. Yet we note that the summary metrics are not of particular interest in single-objective optimization, and obtaining good summary metrics does not always indicate useful results. For example, during the model tuning process, we found out that worse convergence often tend to result in good FCD and diversity score.
Ablation study We performed some ablation studies to better examine the efficacy of each component of our model. Firstly, we segregate spatial graph attention from the RL framework and examine its effect solely in a supervised learning setting. We performed 40 runs with independent samples of 100,000 molecules from the NSP15 dataset, and presented the result in Figure 3. Secondly, we run single-objective optimization on GAPN with (GAPN w/ pre-T) and without layers pre-training, with (DGAPN) and without innovation rewards, as well as the greedy approach on CReM environment. 1,000 molecules are sampled from each setting and the results are shown in Table 2. We note that it is not exactly fair to compare greedy to other approaches since it has access to more information, i.e. the reward for each candidates, when making the decisions, yet our model still beat it in evaluation mode where we take arg max policy instead of sampling policy, expand the number of candidates from 15-20 to 128 and expand the maximum time steps per episode from 12 to 20.
Dock Score
mean 1st 2nd Diversity QED SA
GAPN -7.5 -10.80 -10.75 0.86 0.39 3.14 GAPN w/ pre-T -7.1 -10.24 -10.19 0.86 0.47 2.89 DGAPN -8.3 -12.78 -12.12 0.84 0.36 3.22
CReM Greedy -8.87 -10.60 -10.50 0.85 0.28 2.99 DGAPN eval -9.3 -13.89 -13.33 0.84 0.26 3.24
Table 2: Dock scores and other metrics under different training and evaluation settings
In these experiments, spatial convolution exhibited a strong impact in molecular graph representation learning. From the results by greedy approach, we can see that the environment and the stochasticity design of the action space alone are powerful for the efficacy and exploration of our policies. While the innovation reward helped the discovering of molecules with higher docking scores, it also worsen the two realism metrics, QED and SA. As of layers pre-training, the influence is the exact opposite. We further investigated the balancing of this trade-off in Section 4.2.3. A sample of the molecules generated by DGAPN in evaluation can be found in our supplemental material.

heading: 4.2.2 CONSTRAINED OPTIMIZATION
text: The goal of constrained optimization is to find molecules that have large improvement over a given molecule from the dataset while maintaing a certain level of similarity:
rm′(sT ) = rm(sT )− λ ·max{0, δ − SIM {s0, sT }} (13) where λ is a scaling coefficient, for which we chose 100; SIM {·, ·} is the Tanimoto similarity between Morgan fingerprints. We used a subset of 100 molecules from our dataset as the starting molecules, chose the two most recent and best performing benchmark models in single-objective optimization to compete against, and evaluated 100 molecules generated from theirs and ours. The results are shown in Table 3.
From the results, it seems that MARS is not capable of performing optimizations with similarity constraint. Compared to MolDQN, DGAPN gave better improvements and similarities across all levels of δ, although MolDQN was able to produce molecules with more stable similarity scores.

heading: 4.2.3 MULTI-OBJECTIVE OPTIMIZATION
text: We investigate the balancing between objective and realism by performing multi-objective optimization, and thus provide another approach to generate useful molecules in practice. We weight rm with a function of the two metrics, f(QED,SA). The choice of f and the new main reward are given by
rm′(sT ) = ω · rm(sT ) + (1− ω) · µ1 · [ QED(sT ) + ( 1− SA(sT )− 1
µ2
)] (14)
where µ1 and µ2 are scaling coefficients, for which we chose 8 and 9 respectively. A summary of the results obtained from different ω are shown in Figure 4.
A trade-off between docking reward and QED/SA was identified. We acknowledge that optimizing docking alone does not guarantee finding practically useful molecules, yet our goal is to generate promising chemicals with room for rational hit optimization. We also note that commonly used alternative main objectives such as pLogP and QED are themselves unreliable or undiscerning as discussed in Appendix C. Hence, for methodological study purposes, we believe that docking scores provide a more useful and realistic test bed for algorithm development.

heading: 5 RELATED WORK
text: The REINVENT (Olivecrona et al., 2017) architecture consists of two recurrent neural network (RNN) architectures, generating molecules as tokenized SMILE strings. The “Prior network” is trained with maximum likelihood estimation on a set of canonical SMILE strings, while the “Agent network” is trained with policy gradient and rewarded using a combination of task scores and Prior network estimations. The Junction Tree Variational Autoencoder (JTVAE, Jin et al. (2018)) trains two encoder/decoder networks in building a fixed-dimension latent space representation of molecules, where one network captures junction tree structure of molecules and the other is responsible for fine grain connectivity. Novel molecules with desired properties are then generated using Bayesian optimization on the latent space. Graph Convolutional Policy Network (GCPN, You et al. (2018a)) is a policy gradient RL architecture for de novo molecular generation. The network defines domain-specific modifications on molecular graphs so that chemical validity is maintained at each episode. Additionally, the model optimizes for realism with adversarial training and expert pre-training using trajectories generated from known molecules in the ZINC library. Molecule Deep Q-Networks (MolDQN, Zhou et al. (2019)) is a Q-learning model using Morgan fingerprint as representations of molecules. To achieve molecular validity, chemical modifications are directly defined for each episode. To enhance exploration of chemical space, MolDQN learns H independent Q-functions, each of which is trained on separate sub-samples of the training data. Markov Molecular Sampling (MARS, Xie et al. (2021)) generates molecules by employing an iterative method of editing fragments within a molecular graph, producing high-quality candidates through Markov chain Monte Carlo sampling (MCMC). MARS then uses the MCMC samples in training a GNN to represent and select candidate edits, further improving sampling efficiency.

heading: 6 CONCLUSIONS
text: In this work, we introduced a spatial graph attention mechanism and a curiosity-driven policy network to discover novel molecules optimized for targeted objectives. We identified candidate antiviral compounds designed to inhibit the SARS-CoV-2 protein NSP15, leveraging extensive molecular docking simulations. Our framework advances the state-of-the-art algorithms in the optimization of molecules with antiviral potential as measured by molecular docking scores, while maintaining reasonable synthetic accessibility. We note that a valuable extension of our work would be to focus on lead-optimization — the refinement of molecules already known to bind the protein of interest through position-constrained modification. Such knowledge-based and iterative refinements may help to work around limitations of the accuracy of molecular docking predictions.

heading: A MEASURE THEORY CONSTRUCTION OF MARKOV DECISION PROCESS
text: Let (S,S) and (A,A) be two measurable spaces called the state space and action space; functions Π : S × A → R and T : S × A × S → R are said to be a policy and a transition probability respectively if
1. For each s ∈ S, E → Π(s, E) is a probability measure on (A,A); for each (s, a) ∈ S×A, F → T (s, a, F ) is a probability measure on (S,S).
2. For each E ∈ A, s → Π(s, E) is a measurable function from (S,S) → (R,B); for each F ∈ S, (s, a)→ T (s, a, F ) is a measurable function from (S ×A,S ⊗A)→ (R,B).
We say a sequence of random variable duples (St, At) defined on the two measurable spaces is a Markov decision chain if
P (At ∈ E | σ(S0, A0, S1, A1, . . . , St)) = Π(St, E) (15) P (St+1 ∈ F | σ(S0, A0, S1, A1, . . . , St, At)) = T (St, At, F ) (16)
A function r : S × A → R is said to be the reward function w.r.t. the Markov decision chain if r(st, Et) = EΠ,T [R(st+1) | St = st, At ∈ Et] whereR : S → R is its underlying reward function. With an abuse of notation, we define π(a|s) := Π(s, {a}), ρ(s′|s, a) := T (s, a, {s′}) and let r(s, a) denote r(s, {a}).

heading: B REINFORCEMENT LEARNING ENVIRONMENT AND REWARD EVALUATION DETAILS
text: B.1 ENVIRONMENT - CREM
Chemically Reasonable Mutations (CReM) is an open-source fragment-based framework for chemical structure modification. The use of libraries of chemical fragments allows for a direct control of the chemical validity of molecular substructures and to consider the chemical context of coupled fragments (e.g., resonance effects).
Compared to atom-based approaches, CReM explores less of chemical space but guarantees chemical validity for each modification, because only fragments that are in the same chemical context are interchangeable. Compared to reaction-based frameworks, CReM enables a larger exploration of chemical space but may explore chemical modifications that are less synthetically feasible. Fragments are generated from the ChEMBL database (Gaulton et al., 2012) and for each fragment, the chemical context is encoded for several context radius sizes in a SMILES string and stored along with the fragment in a separate database. For each query molecule, mutations are enumerated by matching the context of its fragments with those that are found in the CReM fragment-context database (Polishchuk, 2020).
In this work, we use grow function on a single carbon to generate initial choices if a warm-start dataset is not provided, and mutate function to enumerate possible modifications with the default context radius size of 3 to find replacements.
B.2 EVALUATION - AUTODOCK-GPU
Docking programs use the three-dimensional structure of the protein (i.e., the receptor) to predict the most stable bound conformations of the small molecules (i.e., its putative ligands) of interest, often targeting a pre-defined functional site, such as the catalytic site. An optimization algorithm within a scoring function is employed towards finding the ligand conformations that likely correspond to binding free energy minima. The scoring function is conformation-dependent and typically comprises physics-based empirical or semi-empirical potentials that describe pair-wise atomic terms, such as dispersion, hydrogen bonding, electrostatics, and desolvation (Huang et al., 2010; Huey et al., 2007). AutoDock is a computational simulated docking program that uses a Lamarckian genetic algorithm to predict native-like conformations of protein-ligand complexes and a semiempirical scoring function to estimate the corresponding binding affinities. Lower values of docking scores (i.e., more negative) indicate stronger predicted interactions (Santos-Martins et al., 2021). The opposite value of the lowest estimated binding affinity energy obtained for each molecule forms the reward.
AutoDock-GPU is an extension of AutoDock to leverage the highly-parallel architecture of GPUs and was implemented in our framework. Within AutoDock-GPU, ADADELTA, a gradient-based method, is used for local refinement (Zeiler, 2012). The structural information of the receptor (here, the NSP15 protein) used by AutoDock-GPU is processed prior to running the framework. In this preparatory step, AutoDockTools (Morris et al., 2009b) was used to define the search space for docking on NSP15 (PDB ID 6W01; Figure 5) and to generate the PDBQT file of the receptor, which contains the atomic coordinates, partial charges, and AutoDock atom types. AutoGrid4 (Morris et al., 2009a) was used to pre-calculate grid maps of interaction energy at the binding site for the different atom types defined in CReM.
In evaluation, after applying an initial filter within RDKit to check whether a given SMILES is chemically valid (i.e., hybridization, ring membership etc.), a 3D conformer of the molecule is generated using AllChem.EmbedMolecule. SMILES that do not correspond to valid compounds are discarded. Next, the molecular geometry is energy minimized within RDKit using the generalized force filed MMFF94. This also provides DGAPN/sGAT the “stable” 3D conformer of the compound, and is used as input for molecular docking via AutoDock-GPU.

heading: C MORE RESULTS ON QED AND PENALIZED LOGP
text: Although QED and penalized LogP are the most popular objectives to benchmark ML algorithms for molecule generation, these benchmarks are questionable for both scientific study and practical use as Xie et al. (2021) pointed out. Most methods can obtain QED scores close or equal to the highest possible of 0.948, making the metric hard to distinguish different methods. As for pLogP, if we simply construct a large molecule with no ring, such as the molecule from SMILES ‘CCCCC...CCCCC’ (139 carbons), it will give us a pLogP score of 50.31 which beats all state-of-the-art models in Table 4. Needless to say, we will achieve a even higher pLogP by continuously adding carbons, which was exactly how REINVENT performed in our experiment. We note that we were able to push our results to around 18 solely by doubling the maximum time step per episode reported in Section 4, yet not so interested in continuing to push the performance on this somewhat meaningless metric by endlessly increasing one hyperparameter.
The results from REINVENT were produced in our own experiments, while others were directly pulled out from the original results reported in the literature.

heading: D QED, SA, AND FCD
text: D.1 QUANTITATIVE ESTIMATE OF DRUGLIKENESS
(QED) is defined as
QED = exp
( 1
n n∑ i=1 ln di
) ,
where di are eight widely used molecular properties. Specifically, they are molecular weight (MW), octanol-water partition coefficient (ALOGP), number of hydrogen bond donors (HBD), number of hydrogen bond acceptors (HBA), molecular polar surface area (PSA), number of rotatable bonds (ROTB), the number of aromatic rings (AROM), and number of structural alerts. For each di,
di(x) = ai + bi
1 + exp ( −x−ci+ di 2
ei
) · 1− 1
1 + exp ( −x−ci+ di 2
fi
)  ,
each ai, . . . , fi are given by a supplementary table in Bickerton et al. (2012).
D.2 SYNTHETIC ACCESSIBILITY
(SA) is defined as

heading: SA = fragmentScore− complexityPenalty
text: The fragment score is calculated as a sum of contributions from fragments of 934,046 PubChem already-synthesized chemicals. The complexity penalty is computed from a combination of ringComplexityScore, stereoComplexityScore, macroCyclePenalty, and the sizePenalty:
ringComplexityScore = log(nRingBridgeAtoms + 1) + log(nSprioAtoms + 1)
stereoComplexityScore = log(nStereoCenters + 1)
macroCyclePenalty = log(nMacroCycles + 1)
sizePenalty = nAtoms1.005 − nAtoms
D.3 FRÉCHET CHEMNET DISTANCE
(FCD) is defined as the distance between a Gaussian distribution pw(.) of real-world molecules with mean and covariance (mw,Cw), and a Gaussian distribution pw(.) of molecules from a generative model with mean and covariance (m,C). FCD is given by
d2((m,C), (mw,Cw)) = ‖m−mw‖22 + Tr (C +Cw − 2(CCw) 1/2).

Review summaries:
The authors propose a new fragment-based method for molecule generation. The model, called DGAPN, uses chemical fragments extracted from a public compound library with their chemical context (atom neighborhoods to which these fragments are attached). This way, all modifications made by the model should produce synthetically accessible compounds with high probability. Additionally, the authors introduce a new class of graph neural networks to predict chemical properties that employ an attention mechanism over atoms and chemical bonds. These models are used to guide the generative MDP trained with reinforcement learning. Inhibition of a SARS-CoV-2 antiviral target, NSP15, is used as an example task for the proposed model. The experimental section shows the results of single- and multi-objective compound generation, in which DGAPN obtains better docking scores of the generated molecules than other methods.

With the goal of generating molecules that bind to functional sites of SARS-Cov-2 protein, the paper proposed a reinforcement learning model with a fragment-based framework for chemical structure modification. In the network part, spatial graph attention and spatial convolution are utilized to extract more structural information from the input graph into the representation of nodes or graphs. Based on the actor-critic algorithm, the reinforcement learning part is designed to find the state with the best docking score computed by docking programs and some novel approaches like PPO and RND are used to train the model more effectively. In the experiments, the model shows great performance while reducing the complexity of chemical synthesis meanwhile.

This paper proposes a method to generate molecular graphs with optimized properties. Molecular graphs are constructed by the iterative addition of molecular fragments in a deep reinforcement learning framework. The method is benchmarked against a set of baselines in a task to generate molecules that maximizes the docking score to a protein (NSP15) from the SARS-CoV-2 virus, and shows good performance.