Differentiable Scaffolding Tree for Molecule Optimization
Authors: Tianfan Fu
Wenhao Gao
Cao Xiao
Jacob Yasonik
Connor W. Coley
Jimeng Sun
Venue: ICLR 2022 Poster (2022)

Abstract: The structural design of functional molecules, also called molecular optimization, is an essential chemical science and engineering task with important applications, such as drug discovery. Deep generative models and combinatorial optimization methods achieve initial success but still struggle with directly modeling discrete chemical structures and often heavily rely on brute-force enumeration. The challenge comes from the discrete and non-differentiable nature of molecule structures. To address this, we propose differentiable scaffolding tree (DST) that utilizes a learned knowledge network to convert discrete chemical structures to locally differentiable ones. DST enables a gradient-based optimization on a chemical graph structure by back-propagating the derivatives from the target properties through a graph neural network (GNN). Our empirical studies show the gradient-based molecular optimizations are both effective and sample efficient (in terms of oracle calling number). Furthermore, the learned graph parameters can also provide an explanation that helps domain experts understand the model output. The code repository (including processed data, trained model, demonstration, molecules with the highest property) is available at https://github.com/futianfan/DST.

Published as a conference paper at ICLR 2022
DIFFERENTIABLE SCAFFOLDING TREE FOR MOLECULAR OPTIMIZATION
Tianfan Fu'-, Wenhao G Cao Xiao², .acob
Jacob Yasonik?, Connor W. Coley? & Jimeng Sun"
Georgia Institute of Technology, MMasaacuusetss Institute of Technology 3 'Ampliude. *University of Illinois Urbana-Cham Chhapa-Caampaign.. "Equal Contribution
cfu428gatech.eda edu {whgao, jyasonik, ccoley @mit edu (danicaxiao, jimeng. sun) @ma il com
ABSTRACT
The structural design of functional molecules, also called molecular optimization, is an essential chemical science and engineering task with important applications, such as drug discovery. Deep generative models and combinatorial optimization methods achieve initial success but still struggle with directly modeling discrete chemical structures and often heavily rely on brute-force enumeration. The chal- lenge comes from the discrete and non-differentiable nature of molecule structures. To address this, we propose differentiable scaffolding tree (DST) that utilizes a learned knowledge network to convert discrete chemical structures to locally dif- ferentiable ones. DST enables a gradient- bbsed optimization on chemical graph structure by back back-propagating the derivatives from the target properties through a graph neural network (GNN). Our empirical studies show the gradient- -based molecular optimizations are both effective and sample efficient (in terms of oracle calling number). Furthermore, the learned graph parameters can also provide an explanation that helps domain experts understand the model output. The code repos- itory (including processed data, trained model, demonstration, molecules with the highest property) is available httptts:: /ttp//ggtttt//tt com/ fut an
INTRODUCTION
The structural design of new funetional molecules, also called molecular optimization, is the key to many scientific and engineering challenges, such as finding energy storage materials (Hachmann al.
201T: Janet et al. 2020), small molecule pharmaceutics (Kuntz, 1992 Zhavoronkov al.
2019
and unvironment--rienddy material (Zimmerman et al.) 2020). The objective is to identify novel molecular structures with desirable chemical or physical properties
Bombarell et al. Dai et al.] 2018 Jin et al. 2018: You et al. 2018 Jin et al. 2019 Shi et al. 2020 [hou et al.
2018:
2019
Fu et al. 2020; Jin et al. 2020 Zang & Wang 2020 Xie et al. 2021). Recent advances
deep generative models (DGM) allow learning the distribution of molecules and optimizing the
latent embedding vectors of molecules. Models in this category are exemplified by the variational
autoencoder (VAE) (Gmmez-Bommarelli al.) 2018 Dai et 2018 Jin al.) [2018] 2020) and generative adversarial network (GAN) (De Cao & Kipl 2018). On the other hand, because of the discrete nature of the enormous chemical space, applying combinatorial optimization algorithms with some structure enumeration has been the predominant approach (You et al. 2018 Jensen 2019 Zhou al. 2019 Nigam et al., 2020, Xie et al., 202T). Deep learning models have also been used to guide these combinatorial optimization algorithms. For example, You et al. (2018); Zhou et al. (2019); Jin al. (2020); Gottipati al.) (2020) tried to solve the problem with deep reinforcement learning; Nigam et al. (2020) enhanced genetic algorithm with neural network discriminator; Xie et al. (2021); Fu et al. (2021) approached the problem with Markov Chain Monte Carlo (MCMC) explore the target distribution guided by graph neural networks.
Despite the initial success of these previous attempts, the following limitations remain: (1) deep generative models optimize the molecular structures a learned latent space, which requires the latent space to be smooth and discriminative. Training such models need carefully designed networks and well-distributed datasets. (2) most combinatorial optimization algorithms, featured by evolutionary

---

Published as a conference paper at ICLR 2022
learning methods (Nigam et al. 2020 Jensen 2019:
Xie al. 2021 Fu et al.
2021 exhibit random-
walk behavior, and leverage trial -aal-and-oee -ndo- strategies to explore the discrete chemical space. The recent deep reinforcement learning methods (You et al. 2018 Zhou et al. 2019; Jin et al. 2020 Coutpatti et al. 2020) aim to remove random-walk search using a deep neural network to guide the searching. It is challenging to design the effective reward function into the objective (Jin et al. 2020).
(3) Most existing methods require a significant number of oracle calls (a property evaluator; see Def. T to proceed with an efficient search. However, realistic oracle functions, evaluating with either experiments or high-fidedt high-iideityy computational simulations, are usually expensive. Examples include using biological assays to determine the potency of drug candidates (Wang et al. (0177) or conducting electronic structure calculation to determine photoelectric properties (Long et al. 2011).
We propose differentiable scaffolding tree (DST) to address these challenges, where we define a differentiable scaffolding tree for molecular structure and utilize trained GNN to obtain the local derivative that enables continuous optimization. The main contributions are summarized as:
We propose the differentiable scaffolding tree to define a local derivative of a chemical graph. This concept enables gradient- based optimization of a discrete graph structure.
We present a general molecular optimization strategy utilizing the local derivative defined by the differentiable scaffolding tree. This strategy leverages the property landscape's s geometrie structure and suppresses the random-walk behavior, exploring the chemical space more efficiently. We also incorporate a determinantal point process (DPP) selection strategy to enhance the diversity of generated molecules.
We demonstrate encouraging preliminary results on de novo molecular optimization with multiple computational objective functions. The local derivative shows consistency with chemical intuition, providing interpretability of the chemical stroecure-rroperty relationship. Our method also requires less oracle calls, maintaining good performance in limited oracle settings.
2 RELATED WORK
Existing molecular optimization methods can mainly be categorized as deep generative models and combinatorial optimization methods.
Deep generative models model a distribution of general molecular structure with a deep network
model so that one can generate molecules by sampling from the learned distribution. Typical algorithms include variational autoencoder (VAE), generative adversarial network (GAN), energy-
based models, flow-baser model BombarlBo al., 2018 Jin et al.,
2018 De Cao & Kipf
2018: Segler et al. 2018 [Tin et al. 2009||Hoaa et al. 2019 Madhawa et al. 2019 Shi et al. 2020
Tin et al. 2020 Zang & Wang [200]] Kotsias et al./ 2020 Chen et al. [022] Fu et al., 2020 Liu
al. 202 Bagal et al. (022)) Shen et al. (2021) also leverages inverse learning based on SELFIES
representation. However, its performance is not satisfactory, primarily due to the failure of training an
adequate surrogate oracle. In addition, DGMs can leverage Bayesian optimization in latent spaces to optimize latent vectors and reconstruct to obtain the optimized molecules (Jin et al. 2018). However, such approaches usually require a smooth and discriminative latent space and thus elaborate network architecture design and well-distributed data set. Also, as they learn the reference data distribution, their ability explore diverse chemical space is relatively limited, evidenced by the recent molecular optimization benchmarks (Brown et al. 2019: Huang et al. 2021).
Conbinatorial optimization methods mainly include deep reinforcement learning (DRL) You et al.) 2018 Zhou et al. 2019 JJin et al., 2020) Gottipati et al. 2020) and evolutionary learning meth- ods (Nigam et al. 22200] Jensen 2019 Xie al. [202] Fu et al. [022]) They both formulate molecule optimization as discrete optimization task. Specifically. they modify molecule substructures (or tokens in a string representation (Weininger) 1988)) locally, with an oracle score or a policy/value network to tell if they keep it or not. Due to the discrete nature of the formulation, most of them conduct an undirected search (random-wal behavior), while some recent ones like reinforcement learning try to guide the searching with a deep neural network, aiming to rid the random-walk nature. However, it is challenging to incorporate the learning objective target into the guided search. Those algorithms still require massive numbers of oracle calls, which is computationally inefficient during the inference time (Kroovina et al. 2020). Our method, DST, falls into this category, explicitly leverages the objective function landscape, and conducts efficient grie-trie search. Instead
2

---

Published as a conference paper at ICLR 2022
Abstract
GNN Predict
Gradient Ascent
Assemble and Diversify
EXPAND
S
UNCHANGE
Molecular Graph
Scaffolding Tree
Differentiable Scaffolding Tree
Molecular Graph
Figure 1: Illustration of the overall approach: During inference, we construct the corresponding scaffolding tree and differentiable scaffolding tree (DST) for each molecule. We optimize each DST along its gradient back propagated from the GNN and sample scaffolding trees from the optimized DST. After that, we assemble trees into molecules and diversify them for the next iteration.
of operating on molecular substructure or tokens, we define the search space as a binary and multinomial variables to indicate the existence and identity of nodes, respectively, and make it locally differentiable with learned GNN as a surrogate the oracle. This problem formulation can find its root in conventional computer-aided molecular design algorithms with samcchandd-buug algorithms as solutions (Sinha et al., 1999 Sahinidis Tawarmalani 2000).
METHOD
We first introduce the formulation of molecular optimization and differentiable scaffolding tree (DST) illustrate the pipeline in Figure T then describe the key steps following the order:
in Section 3.1
We leverage GNNs to imitate property oracles, which are targets of molecular optimization (Section 3.2). Oracle GNN is trained once and for all. The training is separately from optimizing DST below.
Oracle GNN construction:
Optimizing differentiable scaffolding tree: We formulate the discrete molecule optimization into a locally differentiable problem with a differentiable scaffolding tree (DST). Then a DST can be optimized by the gradient back-propp pocc-propggtted from oracle GNN (Section 3.3
Molecule Diversification After that, we describe how we design a determinantal point process (DPP) based method to output diverse molecules for iterative learning (Section 3.4
3.1 PROBLEM FORMULATION AND NOTATIONS
3.1.1 Molecular optimization problem Oracles are the objective functions for molecular optimiza- tion problems, e.g., QED quantifying a molecule's drug-likeness (Bickerton et al. 2012).
Definition 1 (Oracle O). Oracle O is a black-box function that evaluates certain chemical or biological properties of a molecule X and returns the ground truth property O(X).
In realistic discovery settings, the oracle acquisition cost is usually negligible. Suppose optimize P molecular properties specified by oracle O,, Op, we can formulate multi-objective molecule optimization problem through scalarization as represented in Eq. (I).
arg maxxeQ F(X; 02, 0p) ((O,(X), (x))),
where X is molecule, Q denotes the set of valid molecules; is the composite objective combining all the oracle scores, e.g., the mean value of P oracle scores.
3.1.2 Scaffolding Tree The basic mathematical description of molecule is molecular graph, which contains atoms as nodes and chemical bonds as edges. However, molecular graphs are not easy to generate explicitly as graphs due the presence of rings, relatively large size, and chemical validity constraints. For case of computation, we convert molecular graph to scaffolding tree as higher-level representation, tree of substructures, following Jin al. (2018) 2019).
Definition (Substructure). Substructures can be either alom single ring. The substructure sei is denoted S (vocabulary set), covering frequent atoms and single rings in drug-like molecules.
Definition 3 (Scaffolding Tree T). A scaffolding tree, Tx, is a spanning tree whose nodes are substructures. is a higher-level representation molecular graph X.
3

---

Published as a conference paper at ICLR 2022
Scaffolding Tree
1
non-leaf
leaf
expansion
0
Molecular Graph
64
0
des
dgg
0
0
9
0
0
Sag
0
0
9
May 10
da
0
Differentiable Scaffolding Tree
Figure 2: Example of differentiable scaffolding tree. We show non-leaf nodes (grey), leaf nodes (yellow), expansion nodes (blue). The dashed nodes and edges are learnable, corresponding to nodes' identity and existence. w and A share the learnable parameters (W3, W1, Waj3, Wals: W711, W
Tx represented by (i) node indicator matrix, (ii) adjacency matrix, and (iii) node weight vector. We distinguish leaf and non-leaf nodes in Tx. Among the K nodes in Tx, there are Kleat leaf nodes (nodes connecting to only one edge) and Kieat non-leaf nodes (otherwise). The sets of leaf nodes and non-leaf nodes are denoted Vieat and Voleaf correspondingly.
N Naculeaa Neaf
Definition 4. Node indicator matrix N is decomposed as N
E (00 {0,1) where
Nuonteaf E (K- corresponds to non-leaf nodes while Nieaf E {0, 13 Knax|S] corre- sponds to leaf nodes. Each row of N one-hot vector, indicating which substructure the node belongs to.
Definition 5. Adjacency matrix is denoted A {0, [0,1]KxK, Ag = indicates the i-th node and j-th node connected while 0 indicates unconnected.
Definition 6. Node weight vector, w = [1, ,1]T E RR, indicates the K nodes in scaffolding free equally weighted.
3.1.3 Differentiable scaffolding tree Similar to a scaffolding tree, a differentiable scaffolding tree (DST) also contains (i) node indicator matrix, (ii) adjacency matrix, and (iii) node weight vector, but with additional expansion nodes. Specifically, while inheriting leaf node set Veaa and non-leaf node set Viniteat from the original scaffolding tree, we add expansion nodes and form expansion node set, Vespand {u,,/v E Vear Voenteur): [Vespnaal Kexpand K, where u, is connected u in the original scaffolding tree. We also define differentiable edge set, A {(t, e')
Vear OR v' E V Yexpandi vf are connected] to incorporate all the edges involving leaf-nonleaf node and leaf/nonleaff- expansion node connections. To make locally differentiable, we modify the tree parameters from two aspects: (A) node identity and (B) node existence. Figure an example to illustrate DST.
(A) We enable optimization on node identity by allowing the corresponding node indicator matrix learnalle:
Definition 7. Differentiable node indicator matrix Ñ takes the form:
Snonleof Nteof Nexpondd
S M
Ñ
R
BR(K++ (K+
[S|
NN==1.
K= Kespand-
Nunntef N nonteaf {0, 1 x(S| are fixed, equal to the part in the original scaffolding tree, each row is a one-hot vector, indicating that we fix all the non-leaf nodes. In contrast, both N Nexpand and N are learnable, we use softmax activation to implicitly encode the constraint Σ,
No IN are the parameters to learn. This constraint guarantee that each row
i.e. Nij ΣS- .x
exp(NN.))
of N is valid substructures' distribution.
(B) We enable optimization on node existence by assigning learnable weights for the leaf and expansion nodes, constructing adjacency matrix and node weight vector based those values:
K depends molecular graph. During optimization (Section 3.3 and 3.4 after molecular structure
changes, K updated.
4

---

Published as a conference paper at ICLR 2022
Definition 8. Differentiable adjacency matrix A IR Kespure) x( Keequd)) takes the form:
1/0,
(i,j) & A Ordisconnected, 1:connected
Ay Azi o W.
(2,)) E Ai Vient J V Vnontea
(3)
o (Wag)
(1,5) E E, Verpands J E Vieaf U Vunleaf
where A is the differentiable edge set defined above, Sigmoid function o-)) imposes the constraint Kopust 0 S Aij 1. W E KK Koyus are the parameters, each leaf node and expansion node has one learnable parameter. For connected and j. when E Visaf= j E Vnonteofs Aij or (W:) measures the existence probability of leaf node i; when E Vespand E ReeaLVunnteeuf Aij = o (Wa1s) measures the conditional probability of the existence of expand node given the original node When is a leaf node, it naturally embeds the inheritance relationship between the leaf node and the corresponding expansion node.
Definition 9. Differentiable node weight vector w E +K takes the form:
E Vmonleaf
i Vieaf
W;-
i E Verpand J E Vieafs
(1,5) E A,
(4)
(Way) σ Wals),
Verpand= J E V Vuonteafs (i, E.,
where all the weights range from 0 to 1. The weight of the expansion node connecting to the leaf node relies the weight of the corresponding leaf node. w and A (Def. 8) shares the learnable parameter W. Figure 2 shows an example to illustrate DST.
3.2 TRAINING ORACLE GRAPH NEURAL NETWORK
This section constructs a differentiable surrogate model to capture the knowledge from any oracle function. We choose graph neural network architecture for its state-of-the-ar performance in modeling structure-property relationships. In particular, we imitate the objective function F with GNN:
ÿ = GNN(N;; = F(X; O2, ,0p) =
(5)
where ④ represents the GNN's parameters. Concretely, we use a graph convolutional net- work (GCN) (Kipf & Welling 2016). Other GNN variants, such Graph Attention Network (GAT) (Velicković et 2017). Graph Isomorphism Network (GIN) (Xu et al. 2018). can also be used our setting. The initial node embeddings H))) NE E & xd stacks basic embeddings of all the nodes in the scaffolding tree, d is the GCN hidden dimension, N is the node indicator matrix
(Def. 4 E E R is the embedding matrix of all the substructures in vocabulary set S, and is randomly initialized. The updating rule GCN for the [-th layer is
H RE + A(H (A (U()) 1=
where L is GCN's depth, A is the adjacency matrix (Def. 5 H Rxd is the nodes' embedding of layer 1, (()) = [b (()),.. .... b E ERKx ()) E Iddxd are bias/weight parameters.
We generalize the GNN from discrete scaffolding tree differentiable one. Based on learnable weights for each node, we leverage the weighted average as the readout function of the last layer's (L-th) node embeddings, followed by multilayer perceptron (MLP) yield the prediction y, i.e..
to 1, ((()) (L) 2 wa Σk We ()).). in discrete scaffolding tree, weights for all the nodes are equal
T = MLP
is the k-th row of H (4), sum, the prediction can be written as
ÿ GNN(X;;) GNN(Tx (), X Q
where θ- (E) U (B(), U are the GNN's parameters. We train the GNN by minimizing the discrepancy between GNN prediction ÿ and the ground truth y.
e. arg ming ExxMED
Clv
OP), GNN(X;
(8)
where is is loss function, e.g., mean squared error; D is the training set. After training, have GNN parameterized by O, to approximate the black-box objective function F (Eq. Worth mentioning that Oracle GNN is trained once and for all. The training is separately from optimizing DST below.
5

---

Published as a conference paper at ICLR 2022
3.3 OPTIMIZING DIFFERENTIABLE SCAFFOLDING TREE
Overview With a little abuse of notations, via introducing DST, we approximate molecule optimiza- tion locally difserentiable problem
X = arg x F
((-11) arg x F(X)
(1) structured combinatorial optimization
(II) iterative local discrete search
(9)
Txain arg (XXX) GNN (Nx A Axt, (x
(TT) local differentiable optinization
where (() is the molecule at f-th iteration, .V S Q is the neighborhood set of ( (Def. 10
Next, we explain the intuition behind these approximation steps. Molecular optimization is generally a discrete optimization task, which is prohibitively expensive due to exhaustive search. The first approximation is to formulate the problem
iterative local discrete search via introducing neighborhood molecule set ((x (9). Second, to enable differentiable learning,
GNN to
imitate black-box objective F' (Section (.2) and further reformulated it into a local differentiable
optimization problem. Then we can optimize DST (Tx = {N , A Axa, Wx }) in continuous domain for (xx using gradient- optimization method.
3.3.1 Local Editing Operations For a leaf node in the scaffolding tree, we can perform three editing operations, (1) SHRINK: delete node v; (2) REPLACE: replace a new substructure over u; (3) EXPAND: add new node u, that connects to node For a nonleaf node v. support (1) EXPAND: add new node connecting to (2) do nothing. If we EXPAND and REPLACE, the new substructures are sampled from the vocabulary S. We define molecule neighborhood set:
Definition 10 (Neighborhood set). Neighborhood set of molecule X. denoted N (X)), is the set of all the possible molecules obtained by imposing one local editing operation to scaffolding tree Tx and assembling the edited frees into molecules.
3.3.2 Optimizing DST Then within the domain of neighborhood molecule set (XX), the objective
function can be represented differentiable function of X's DST (N x, A x, Wx). We address the following optimization problem to get the best scaffolding tree within (XX).
N.. A., w, = arg max (Ñx.A. ..x) GNN( x. Ax, wx +.).
(10)
where the GNN parameters +. (Eq. (8) are fixed. Comparing with Eq. (7. it is differentiable with regard to (N, A, } for all molecules in the neighborhood set ((X). Therefore, we can optimize the DST using gradient-based optimization method, e.g., an Adam optimizer (Kingma & Ba 2015). The
whole DST pipeline leverages iterative local discrete search (Equation Algorithm 1). specifically, in i-th iteration, we optimize the DST of ( i.e., X in Equation (TOp is (().
3.3.3 Sampling from DST Then we sample the new scaffolding tree from the optimized DST. Concretely, (i) for each leaf node v E Vieat and the corresponding expansion node u, Vexpand, we select one of the following step with probabilities (w.p.) as follows,
T DST-SamplerrN...... A., W.)
1. SHRINK: delete leaf node v,
w.p. (W.)es
2. EXPAND: add uy, select substructure at u, based on (N+)u.,
w.p. (w.)o (Wa)) )
3. REPLACE: select substructure at u based on (N=)a.
w.p. (W., )o (1- (w. )..l
(11)
(ii For each nonleaf node v, we expand new node with probability (w. (W))allo ju. If expanding, select substructure at u. based on (Ñ.)... > Otherwise, we "do nothing".
3.3.4 Assemble Each scaffolding tree corresponds to multiple molecules due to the multiple ways substructures can be combined. We enumerate all the possible molecules following Jin et al. (2018) (See Section C.5 in Appendix for details) for the further selection described below.
3.4 MOLECULE DIVERSIFICATION
In the current iteration, we have generated M molecules (X,,... XXM) and need to select C molecules for the next iteration. We expect these molecules to have desirable chemical properties

---

Published as a conference paper at ICLR 2022
(high F score) and simultancously maintain higher structural diversity. To do so, resort to the determinantal point process (DPP) (Kulesza & Taskar 2012), which models the repulsive correlation x between data points. Specifically, for M data points, whose indexes are {1, 2, ]. denotes the similarity matrix between these data points. To create diverse subset (denoted R) with fixed size C, the sampling probability should be proportional to the determinant of the submatrix SR RCxC, iee. P(R) det(SR). where R (1,2, ,,, [R] = C. Combining the objective (F) value and diversity, the composite objective is
DPP (R) = Σ F(X,) log P(R) = log det (R)) + log dct(SR), (12)
arg max RC(122," MIRRI=
PER
where the hyperparamter À > 0 balances the two terms, the diagonal scoring matrix V (iag(eexpAAF(x))),.. ".. )))) VR RCx' is sub-matrix of V indexed by R. When À goes to infinity. it is equivalent to select top-C candidates with the highest F score, same as conven- tional evolutionary learning in (Jensen| 2019 Nigam et al.] (020). Inspired by generalized DPP meth- ods (Chen et al. 2018), we further transform LDPP(R), (R), LDp ( log det VR) ) + log (Sx) =
log det (V SR R = log det
) where V is symmetric positive semi-definin
Then it can be solved by generalized DPP methods in ((C*M) (Chen et al. 2018 (Sec. F The computational complexity of DST is O(TMC²) (Sec. C.8 Algorithm summarizes the pipeline.
Algorithm 1 Differentiable Scaffolding Tree (DST)
1: Input: Iteration T, beam width C, input molecule (1). $ = (x(1).
2: Output: Generated Molecule Set f2
3: Learn GNN (Eq. 8): O = arg ming E(x,))D L(y. 1) # Section 3.2
4: 10 ... do
5: Initialize set O.
6: for (1)) E & do
7: Initialize DST (Ň, A. } for (9) (Eq. 2
8: Optimize DST: N., A., W. arg max (Ñ. GNN A, : +.) (Eq 10
9: Sample from DST: T (t+1) DST-Sampler [N., A., w. ), = 1,2, (Bq
10: Assemble scaffolding tree (++1) into molecules x((+)), r r= rx
11: end for
12: Select & C T, |6| = C' based on Eq. 12
Ω = RUU $
# Section 3.4
13:
end for
4 EXPERIMENT
4.1 EXPERIMENTAL SETUP
Molecular Properties contains QED; LogP; SA: JNK3; GSK333, following (Jin al. 2020 Nigam et al., (2020; Moss et al., [2020] Xie et al. 2021). where QED quantifies likg-lieeness; LogP indicates the water-octanol partition coefficient; SA stands for synthetic accessibility and used prevents the formation of chemically unfeasible molecules; JNK3/GSK33 measure inhibition against c-Jun N-terminal 3/llyog3 synthase kinase beta. For all the scores (including normalized SA), higher is better. We conducted (1) single-objective generation that optimizes JNK3, GSK33 and LogP separately and (2) multi-objective generation that optimizes the mean value of "JNK+GGK33'" and "0ED+SA+INK++GSSS
in the main text. Details are in Section
Dataset: ZINC 250K contains around 250K druglike molecules (Sterling & Irwin 2015). We select the substructures that appear more than 1000 times in ZINC 250K as the vocabulary set S, which contains 82 most frequent substructures. Details are in Section C.T The code repository (including processed data, trained model, demonstration, molecules with the highest property) is available?
Baselines. (1) LigGPT (string-based distribution learning model with Transformer as a decoder) (Ba gal et al. 2021) (2) GCPN (Graph Convolutional Policy Network) (You al. 2018); (3) MolDQN
https://gi thtp:///gttttuttt/ com http://////ftaiaaa///sst
7

---

Published as a conference paper at ICLR 2022
(Molecule Deep Q-Network)
(Zhou et al. 2019;; ; (4) GA+D (Genetic Algorithm with Discriminator
network) (Nigam et al.
2020 (5) MARS (Markov Molecular Sampling) (Xie et al. 2021); (6)
RationaleRL Jin et al. 2020); (8) BOSS (Bayesian Optimization over String Space) (Moss et al. 2020); 99LTM (Long short term memory) (Brown al.] 2019); (10) Graph-GA (graph level genetic algorithm) (Brown ct al. 2019). Among them, LigGPT belongs to the deep generative model, where all the oracle calls can be precomputed; GCPN, MolDQN are deep reinforcement learning methods; GA+D, MARS, Graph-GA are evolutionary learning methods; RationaleRL is a deep generative model fine-tuned with RL techniques. ChemBO and BOSS are Bayesian optimization methods. We also consider DST variant: DST-rand, Instead of optimizing and sample from DST, DST-rand leverages random local search, i.e., randomly selecting basic operations (EXPAND, REPLACE, SHRINK) and substructure from the vocabulary. We also select subset of all the random samples with high surrogate GNN prediction scores to improve efficiency. All the baselines except LigGPT require online oracle calls. Details are in Scction[B]
2020); (7) ChemBO (Chemical Bayesian Optimization) (Korovina et al.
Metries. For each method, we select top-100 molecules with highest property scores for evaluation, and consider the following metrics following Jin et al. (2018); You et al. (2018); Din et al. (2020);
Xie et al.] (2021) (I) Novelty (Nov) (% of the generated molecules that are not in training set); (2)
Diversity (Div) (average pairwise Tanimoto distance between the Morgan fingerprints); (3) Average Property Score (APS) (average score of top- molecules); (4) # of oracle calls: DST needs to
call oracle labeling data for GNN (precomputed) and DST based de novo generation (online), we
show the costs for both steps. For each method in Table | and 2 we set the number of oracle calls so
that the property score nearly converge w.r.t. oracle call's number. Details are in Section C.4
INK34GSK3/8
OED-SA-INK3+CSS3//
Table I: Multi-objective de novo design. #oracle= precommppde oracle call" (to label molecules in existing database) + online oracle call" (during learning).
Matiod
Novi
Divf
APST
eoraclej
Nove
Divt
APS
#oralee
LigGiPT
100%
0.345
0.271
100000
100%
0.9902
0.378
1000-0
GCPN
100%
0.578
0.293
D4200K
100%
0.596
0.450
0+200K
MalDON
100%
0.605
0.348
B+200K
100%
(.597
0+200K
GA-D
100%
0.657
0.608
0450K
97%
0.681
0+50K
RationaleR.
100%
0.700
0.795
25K-67K
90%
0.720
0.675
25K-66K
MARS
100%
0.711
0.789
0450K
100%
0.714
0.662
0+-0KK
ChemBo
0,702
0+50K
0.701
0:50K
ROSS
09%
0.564
0.504
0450K
0.561
0.504
0-50K
LSTM
100%
0.712
0.680
0+30K
100%
0.706
0.672
0-30K
Gmpb-GA
100%
0.634
0.825
5+25K
100%
0.723
0.714
5+25K
DST-ned
0.456
0.765
20K-5K
DS7
100%
0.750
0.827
10K-5K
100%
0.755
20K+5K
4.2 OPTIMIZATION PERFORMANCE
The results of multi-objective and single-objective generation shown in Table I and 2 We find that DGM (LigGPT) and RL-based methods (GCPN and MolDQN) fail in some tasks, which is consistent with the results reported in RationaleRL (Jin et al. 2020 and MARS (Xie et al., 2021). Overall, DST obtains the best results in most tasks. In terms of success rate and diversity. DST outperformed all baselines in most tasks. It also reached the highest scores within T 50 iterations in most optimization tasks (see Table 5 and 6 in Appendix). Especially in optimizing LogP, the model successfully learned to add six-member ring (see Figure in Appendix) each step. which is theoretically the optimal strategy under our setting. Combined with the ablation study compared with random selection (see Figure 12 in Appendix), our results show the local gradient defined by DST is a useful direction indicator, consistent with the concept of gradient. Further, achieving high diversity validates the effect of the Dap-based selection strategy. Although the novelty is not the highest, is still comparable to baseline methods. These results show our gradient-based optimization strategy has a strong optimization ability to provide a diverse set of molecules with high objective functions.
4.3 ORACLE EFFICIENCY
As mentioned above, oracle calls for realistic optimization tasks can be time-consuming and expensive. From Table and 2 we can see that majority of de novo optimization methods require oracle calls online (instead of precomputation), including all of RL/evolutionary algorithm based baselines. DST takes fewer oracle calls compared with baselines. DST can leverage the precomputed oracle calls to label the molecules in an existing database (i.e., ZINC) for training the oracle GNN and dramatically saving the oracle calls during reference. In the three tasks in Table two-thirds of the oracle calls
8

---

Published as a conference paper at ICLR 2022
Table Single-objective de novo molecular generation.
GSK3,9
LegP
Nowt
Divl
APS1
foralee
Nevt
Drrt
APS:
fomcie]
Nevi
Divt
AP81
Poracled
TigeFFT
1047
0.337
0.307
[0K40
100%
0.867
0.255
10.00:0
F0055
0.568
4.55
Do0Ke
GCPN
100%
0.554
0.365
D+200K
100%
0.519
0,400
0+2008
100%
0.532
5.43
0-200K
MoIDON
100%
0.605
0.499
B4200K
100%
0.545
0.39%
0+200%
1005
0.455
6,00
0-300K
GA+D
99%
0,702
0.15
D-50K
0.687
0.678
D+50X
100%
0,221
30.2
O+50K
30K-45K
ChanBB
9886
0&50K
BOSS
98%
0.901
0,471
0+50K
0,432
D+50X
1077
9:84
0+50K
LSTM
100%
0.745
0.645
B+20K
100%
0.715
0.674
0+20X
100%
0,702
10.1
0=30K:
Ginll-AA
100%
0.674
0.865
S-10K
100%
0.785
0.875
5+
1005
0.721
45.0
SaiOK
DS7 rand
100%
0.754
0.413
10K+10K
9016
0.793
0.455
10K-10K
100%
0.715
36.1
[0K+15K
DS?
100%
0.732
0.928
10K45K
100%
0.743
0.869
10K+5K
1007
0,704
47.1
10K+5K
(10K) can be precomputed or collected from other sources. To further verify the oracle efficiency, we explore a special setting molecule optimization where the budget of oracle calls is limited to fixed number (2K, 5K, 10K, 20K, 50K) and compare the optimization performance. For GCPN, MolDQN, GA+D, and MARS, the learning iteration number depends on the budget of oracle calls. RationaleRL (Tin al. 2020) is omitted because it requires intensive oracle calls collect enough reference data, exceeding the oracle budget in this scenario. DST, we use around 80% budget to label the dataset (i.e., training GNN) while the remaining budget to conduct de novo design. Specifically, for 2K, 5K, 10K, 20K, 50K, we use 1.5K, 4K, 8K, 16K, and 40K oracle calls to label the data for learning GNN, respectively. We show the average objective values of top- molecules under different oracle budgets in Figure 3 Our method has significant advantage compared to all the baseline methods in limited budget settings. We conelude that supervised learning a well-studied and much easier task than generative modeling.
JNK3
GSK38
JNK3+6SK3B
DST
RationaleRL
-MARS
BOSS
GA+D
ChemB0
MoIDON
DST-rend
08
GCPN
Figure Oracle efficiency test. Top-100 average score v.s. the number of oracle calls.
4.4 INTERPRETABILIA ANALYSIS
lesf node
gradied 0.001
14:3.10
Positive weight gradient expand node means positive contribution to the objective value increment (at least locally), while negative gradient represents negative one, Due to space limit, we only show the gradient of two expand nodes.
2883:8.91
1508:8.722
Iopand node
ansca:: E401
0.38
9.79
0.12
SHRINK
leaf node
gadlan 0,003
300319.00
0938.0.13
Figure 4: Two steps in optimizing "OED+SAHINK++GSK33
To obtain more insights from the local gradient, We visualize modification steps in Figure 4 The node weights and gradient values interpret the property at the substructure level. This is similar ar to most attribution-ban based interpretable ML methods, e.g., saliency map (Simonyan et al. 2013).
CONCLUSION
This paper proposed differentiable scaffolding tree (DST) to make a molecular graph locally dif- ferentiable, allowing a continuous gradient-baan based optimization. To the best our knowledge, is the first attempt make the molecular optimization problem differentiable at the substructure level, rather than resorting to latent spaces or using RL/evolutionary algorithms. We constructed general molecular optimization strategy based on DST, corroborating thorough empirical studies.

---

Published as a conference paper at ICLR 2022
REPRODUCIBILITY STATEMENT
The code repository is available at http::://gtthub..com/ttttte README file, all the codes for data preprocessing. model architecture, molecules with the highest property. All the data we use are publicly available. Section [ describes an experimental setup for all the baseline methods. We elaborate on the implementation details in Section Concretely, Section C.T describes the dataset we use; Section C.2 describes software and hardware configuration; Section introduces the target molecular properties optimize; Section C.4 describes all the metrics for evaluating the performance; Section describes the setup for model implementation, especially hyperparameter setup.
including
ACKNOWLEDGEMENT
This work supported by NSF award SCH-2014438, PPoSS 2028839, IIS-1838042, NIH award ROI IR0INS107291-01, OSF Healthcare, ONR under grant number N00014-21-1-2195, and the MIT UROP office.
REFERENCES
Viraj Bagal, Rishal Aggarwal, PK Vinod, and U Deva Priyakumar. LigGPT: Molecular generation using a transformer-decoder model. 2021.
G Richard Bickerton, Gaia V Paolini, Jérémy Besnard, Sorel Muresan, and Andrew L. Hopkins. Quantifying the chemical beauty of drugs. Nature chemistry, 4(2):90, 2012.
Rasmus Bro and Age K Smilde. Principal component analysis. Analytical methods, 6(9):2812-288 2831, 2014.
Nathan Brown, Marco Fiscato, Marwin HS Segler, and Alain C. Vaucher. Guacamol: benchmarking models for de molecular design. Journal of chemical information and modeling, 59(3): 1096-1108, 2019.
Krzysztof Chalupka, Christopher KI Williams, and lain Murray. A framework for evaluating approximation methods for gaussian process regression. Journal of Machine Learning Research, 14:333-350, 2013.
Binghong Chen, Tianzhe Wang, Chengtao Li, Hanjun Dai, and Le Song. Molecule optimization by explainable evolution. In International Conference Learning Representations, 2021.
Laming Chen, Guoxin Zhang, and Hanning Zhou. Fast greedy map inference for determinantal point process to improve recommendation diversity. In Neural Information Processing Systems, pp. 5627-5638, 2018.
Sangwoo Cho, Logan Lebanoff, Hassan Foroosh, and Fei Liu. Improving the similarity measure of determinantal point processes for extractive multi-document summarization. In Association for Computational Linguistics, ACL, 2019.
Hanjun Dai, Yingtao Tian, Bo Dai, Steven Skiena, and Le Song. Syntax-directed variational autoencoder for structured data. 2018.
Nicola De Cao and Thomas Kipf. Molgan: An implicit generative model for small molecular graphs. arXiv preprint arXiv: 18885./993... 2018.
Peter Ertl and Ansgar Schuffenhauer. Estimation of synthetic accessibility score of drug-like molecules based on molecular complexity and fragment contributions. Journal of cheminfor- matics, 1(1):8, 2009.
Tianfan Fu, Cao Xiao, and Jimeng Sun. CORE: Automatic molecule optimization using copy and refine strategy. AAAI, 2020.
Tianfan Fu, Cao Xiao, Xinhao Li, Lucas Glass, and Jimeng Sun. MIMOSA: Multi-constrrint molecule sampling for molecule optimization. AAAI, 2021.
10

---

Published as a conference paper at ICLR 2022
Wenhao Gao and Connor W Coley. The synthesizability of molecules proposed by generative models. Journal of chemical information and modeling, 60(12):5714-572 2020.
Rafael GómezzBombarrell, Jennifer N Wei, David Duvenaud, José Miguel Hernández-Lobato
Benjamín Jânchez-Lengelling Dennis Sheberla, Jorge Aguilera- paarragurre,,,,reee, Timothy D Hirzel, Ryan P Adams, and Alán Aspuru-Guzik, Automatic chemical design using a data-driven continuous representation of molecules. ACS central science, 2018.
Sai Krishna Gottipati, Boris Sattarov, Sufeng Niu, Yashaswi Pathak, Haoran Wei, Shengchao Liu, Simon Blackburn, Karam Thomas, Connor Coley, Jian Tang, et al. Learning to navigate the synthetically accessible chemical space using reinforcement learning. In International Conference on Machine Learning, pp. 3668-3679, PMLR, 2020.
Johannes Hachmann, Roberto Olivares-Amaya Sule Atahan-Evrenk, Carlos Amador-Bedolle Roel S Sánchez- Carrera, Aryeh Cold-Parker, Leslie Vogt, Anna M Brockway, and Alán Aspuru-Guzie Guzik. The harvard clean energy project: large-scale computational screening and design of organic photovoltaics on the world community grid. The Journal of Physical Chemistry Letters, 2(17): 2241-22511 2011.
Shion Honda, Hirotaka Akita, Katsuhiko Ishiguro, Toshiki Nakanishi, and Kenta Oono. Graph residual flow for molecular graph generation. arXiv preprint 1909 2019.
Kexin Huang, Tianfan Fu, Wenhao Gao, Yue Zhao, Yusuf Roohani, Jure Leskovec, Connor W Coley, Cao Xiao, Jimeng Sun, and Marinka Zitnik. Therapeutics data commons: machine learning datasets and tasks for therapeutics. Proceedings of Neural Information Processing Systems, NeurIPS Datasets and Benchmarks, 2021.
Jon Paul Janet, Sahasrajit Ramesh, Chenru Duan, and Heather J Kulik. Accurate multiobjective design in a space of millions of transition metal complexes with neural-network-driivee efficient global optimization. ACS central science, 6(4):513-524 2020.
Jan Jensen. A graph-based genetic algorithm and generative model/monte carlo tree search for the exploration of chemical space. Chemical science, 10(12):3567-3572) 2019.
Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Junction tree variational autoencoder for molecular graph generation. ICML, 2018.
Wengong Jin, Kevin Yang, Regina Barzilay, and Tommi Jaakkola. Learning multimodal graph-to- graph translation for molecular optimization. ICLR, 2019.
Wengong Jin, Regina Barzilay, and Tommi Jaakkola. Multi-objective molecule generation using interpretable substructures. In International Conference on Machine Learning, pp. 4859. PMLR, 2020.
Diederik P Kingma and Jimmy Ba. Adam: A method for stochastic optimization. International Conference on Learning Representations (ICLR), 2015.
Thomas N Kipf and Max Welling. Semi-supervised classification with graph convolutional networks. International Conference on Learning Representations (ICLR), 2016.
Ksenia Korovina, Sailun Xu, Kirthevasan Kandasamy, Willie Neiswanger, Barnabas Poczos, Jeff Schneider, and Eric Xing. Chembo: Bayesian optimization of small organic molecules with syn- thesizable recommendations, In International Conference on Artificial Intelligence and Statistics, pp. 3393-3403. PMLR, 2020.
Panagiotis-Christos Kotsias, Josep Arús-Pous, Hongming Chen, Ola Engkvist, Christian Tyrchan, and Esben Jannik Bjerrum. Direct steering of de novo molecular generation with descriptor conditional recurrent neural networks. Nature Machine Intelligence, 255):254-265. 2020.
Alex Kulesza and Ben Taskar. Determinantal point processes for machine learning. arXiv preprint /207.662776 2012.
Irwin D Kuntz. Structure-based strategies for drug design and discovery. Science, 257(3073): 1078-1082, 1992.
11

---

Published as a conference paper at ICLR 2022
Greg Landrum et al. RDKit: Open-source cheminformatics, 2006.
Yibo Li, Liangren Zhang, and Zhenming Liu. Multi-objective de novo drug design with conditional graph generative model. Journal of cheminformatics, 2018.
Meng Liu, Keqiang Yan, Bora Oztekin, and Shuiwang Ji. Graphebm: Molecular graph generation with energy-based models. arXiv preprint arXiv:2 2021.
Mengqiu Long, Ling Tang, Dong Wang, Yuliang Li, and Zhigang Shuai. Electronic structure and carrier mobility in graphdiyne sheet and nanoribbons: theoretical predictions. ACS nano, 5(4): 2593-2600, 2011.
Kaushalya Madhawa, Katushiko Ishiguro, Kosuke Nakago, and Motoki Abe. Graphnyp: An invertible flow model for generating molecular graphs. arXiv preprint
2019.
Henry B Moss, Daniel Beck, Javier Gonzalez, David S Leslie, and Paul Rayson. Boss: Bayesian optimization over string spaces. NeurIPS, 2020.
AkshatKumar Nigam, Pascal Friederich, Mario Krenn, and Alán Aspuru-Guzis. Augmenting genetic algorithms with deep neural networks for exploring the chemical space. 2020.
Mariya Popova, Mykhailo Shvets, Junier Oliva, and Olexandr Isayev. Moleculamnn: Generating realistic molecular graphs with optimized properties. arXiv preprint arXiv:1905.13372 2019.
Nikolaos V Sahinidis and Mohit Tawarmalani. Applications of global optimization to process and
molecular design. Computers & Chemical Engineering, 2409-10)2157-2166) 2000.
Marwin HS Segler, Thierry Kogej, Christian Tyrchan, and Mark P Waller. Generating focused molecule libraries for drug discovery with recurrent neural networks. ACS central science, 4(1): 120-131, 2018.
Cynthia Shen, Mario Krenn, Sagi Eppel, and Alan Aspuru-Guzia Deep molecular dreaming: Inverse machine learning for de-novo molecular design and interpretability with surjective representations. Machine Learning: Science and Technology, 2021.
Chence Shi, Minkai Xu, Zhaocheng Zhu, Weinan Zhang, Ming Zhang, and Jian Tang. GraphAF: flow-based autoregressive model for molecular graph generation. 2020.
Karen Simonyan, Andrea Vedaldi, and Andrew Zisserman. Deep inside convolutional networks:
Visualising image classification models and saliency maps. arXiv preprint vrXiv:1312.6033. 2013.
Manish Sinha, Luke EK Achenie, and Gennadi M Ostrovsky. Environmentally benign solvent design by global optimization. Computers & Chemical Engineering,
1999.
Teague Sterling and John J Irwin. Zinc 15-ligand discovery for everyone. Journal of chemical information and modeling. 55(1 (5(11):2324-233 2337, 2015.
Paraskevi Supsana, Theodoros Liaskopoulos, Stavroula Skoulika, Antonios Kolocouris, Petros G Tsoungas, and George Varvounis. Thermal rearrangement of spiro [maphthalene (naphthopyranofu- razan)] oxides spiro (naphthalene (phenalenofurazzan) oxides. probable furazan oxide triggered tandem isomerisation process. Telrahedron, 61 6125):6131-6133 2005.
Petar Veličković, Guillem Cucurull, Arantxa Casanova, Adriana Romero, Pietro Lio, and Yoshua Bengio. Graph attention networks. arXiv preprint arXiv: 10900, 2017.
Yanli Wang, Stephen H Bryant, Tiejun Cheng, Jiyao Wang, Asta Gindulyte, Benjamin A Shoemaker, Paul A Thiessen, Siqian He, and Jian Zhang. Pubchem bioassay: 2017 update. Nucleic acids research, 45(51))D955-D963, 2017.
David Weininger. Smiles, a chemical language and information system. 1. introduction to methodol- ogy and encoding rules. Journal of chemical information and computer sciences, 28(1):31-36. 1988.
Yutong Xie, Chence Shi, Hao Zhou, Yuwei Yang, Weinan Zhang, Yong Yu, and Lei Li. MARS: Markov molecular sampling for multi-objective drug discovery. 2021.
12

---

Published as a conference paper at ICLR 2022
Keyulu Xu, Weihua Hu, Jure Leskovec, and Stefanie Jegelka. How powerful graph neural networks? arXiv preprint arXiv:1810.00826, 2018.
Jiaxuan You, Bowen Liu, Rex Ying. Vijay Pande, and Jure Leskovec. Graph convolutional policy network for goal-directed molecular graph generation. In NIPS, 2018.
Chengxi Zang and Fei Wang. MoFlow: an invertible flow model for generating molecular graphs. In ACM SIGKDD, pp. 617-626, 2020.
Cheng Zhang, Hedvig Kjellstrom, and Stephan Mandt. Determinantal point processes for mini-batch diversification. Uncertainty in Artificial Intelligence (UAI), 2017.
Alex Zhavoronkov, Yan Ivanenkov, Alex Aliper, Mark S Veselov, Vladimir A Aladinskiy, Anas- tasiya V Aladinskaya, Victor Terentiev, Daniil A Polykovskiy, Maksim D Kuznetsov, Arip Asadulaev, et al. Deep learning enables rapid identification of potent ddrl kinase inhibitors. Nature biotechnology. 3739):1038-1044 2019.
Zhenpeng Zhou, Steven Kearnes, Li Li, Richard Zare, and Patrick Riley. Optimization of molecules via deep reinforcement learning. Scientific reports, 2019.
Julie B Zimmerman, Paul T Anastas, Hanno C Erythropel, and Walter Leitner. Designing for green chemistry future. Science, 367(4476)33774400, 2020.
13

---

Published as a conference paper at ICLR 2022
CONTENTS
Introduction
Related Work
3
Method
Problem Formulation and Notations
3.2
Training Oracle Graph Neural Network
3.3
Optimizing Differentiable Scatfolding Tree
3.4
Molecule Diversification
Experiment
4.1
Experimental Setup
Optimization Performance
4.3
Oracle Efficiency
8
4.4
Interpretability Analysis
5
Conclusion
Complete Mathematical Notations.
16
Baseline Setup
16
Implementation Details
18
Dataset
18
Software/Hardware Configuration
18
Target molecular properties)
18
Evaluation metrics
19
Assembling Molecule from Scaffolding Tree
20
Details on GNN Learning and DST Optimization
20
Results of Different Random Seeds
21
Complexity Analysis
22
Additional Experimental Results
22
D.I
Additional results of de novo molecular generation
22
D.2
De novo molecular optimization on QED (potential limitation of DST
23
D.3
Results Analysis for Distribution Learning Methods (LigGPT)
24
D.4
Ablation study
25
D.5
Chemical space visualization
25
D.6
Additional Interpretability Analysis
26
14

---

Published as a conference paper at ICLR 2022
Theoretical Analysis
26
Convergence Analysis
27
F
Extension of Molecule diversification
27
G
Proof of Theoretical Results
29
G.
Proof of Lemma
29
G.2 Proof of Lemmal
30
G.3 Proof of Lemmal
31
Proof of Theorem
31
APPENDIX TO DIFFERENTIABLE SCAFFOLDING TREE FOR MOLECULAR OPTIMIZATION
The appendix is organized follows. First, we list the complete mathematical notations for ease of exposition. Then, we show additional experimental setup and empirical results, including baseline setup in Section implementation details of our method in Section C additional experimental results in Section D Then, we provide theoretical analysis in Section E extend molecule diversification in Section, and prove the theoretical results in the main paper in Section G
Table 3: Complete Mathematical Notations.
Notations
Descriptions
Oracle function, e.g. evaluator of molecular property (Def
objective function of molecule generation (Eq Number of larget oracles.
P N.
Q
Set of all the valid chemical molecules. Vocabulary set, i... substructure set. A substructure is an atom or ring.
T
Scaffolding tree (Del
K = IT1
number nodes scaffolding tree T.
A;w
Node indicator matrix; adjacency matrix; node weight. Leaf node set in scaffolding tree T.
Vier
Vicntear
Nonleaf node set in scaffolding tree T. Expansion node set in scaffolding tree T. Size of node set. Size of expansion node Kiest Kespans-
Vespast
Kiaf = Vial
Kespand
GNN hidden dimension.
L N
GNN depth.
(E) (B U Mem
Leamable parameter GNN.
E lsika
embedding stackings of all the substructures vocabulary set bias parameters at /-th layer.
BdD Ru
0 c Rexa
weight parameters at (-th layer.
Node embedding 1-th layer GNN
()) NE Raxd
initial node embeddings, stacks basic embeddings of all nodes scaffolding tree. multilayer perceptron
MLP
ReLU
ReLU activate function
GNN prediction. groundtruth Loss function GNN. the truining Neighborhood molecule of X (Der diflerentiable edge set.
N(X)
Ñ: A;
Differentiable node indicator matrix; adjacency matrix; node weight.
del()
Determinant of square matrix
M N
Number of all possible molecules select.
C CEN
Number of selected molccules.
S MMN
Similariyy kernel matrix.
V MXM
Diagonal scoring matrix.
R
subset of [1, ... M). index of select molecules. hyperparameler in Eq. and
A>0 0
balances desirable property and diversity.
15

---

Published as a conference paper at ICLR 2022
A COMPLETE MATHEMATICAL NOTATIONS.
In this Section, we show all the mathematical notations in Table 3 for completeness.
We also illustrate the difference between DST and existing methods in Figure
Existing methods (GCPN, MolecularRNN, GraphAF)
DST (differentiable scaffolding tree)
(2) predict property
(1) generate molecule
GNN generator
(2) evaluate
(1) DST
Pretrained GNN as property predictor
maximize predicted property
objective
(differentiable scaffolding tree)
(3) Gradient BP, Update GNN
(4) Gradient BP, update DST
Figure 5: Left: Most of the existing methods (including GCPN (You et al. 2018), Molecular-
RNN (Popova al. 2019), GraphAF (Shi al. 2020) use GNN as graph generator. Specifically,
(1) generate molecule; (2) evaluate learning objective (loss in deep generative model or reward in reinforcement learning); (3) back-propagate (BP) gradient to update GNN. In sum, the learning objective is differentiable w.r.t. the GNN's parameters.
Right: Regarding DST, given the pretrained GNN as surrogate oracle model (i.e., property predictor), we have several steps: (1) construct differentiable scaffolding tree (DST); (2) predict the property via GNN: (3) maximize the predicted property ý (learning objective); (4) pack-propagate (BP) gradient to update DST. In sum, the learning objective is differentiable w.r.t. DST (also input of GNN). Summary: DST makes the learning objective differentiable w.r. molecule graph structure, while prior works make the learning objective differentiable w.r.t. neural networks' parameters. Our approach directly optimizes molecular graph structures, while prior works indirectly search for the molecule graphs with the help of neural network.
B BASELINE SETUP
In this section, we describe the experimental setting for baseline methods. Most of the settings follow the original papers.
LigGPT (string- -aaieg--bssed distribution learning model with Transformer as a decoder) (Bagal et al.] 2021) is trained for 10 epochs using the Adam optimizer with a learning rate of 6e 4. LigGPT comprises stacked decoder blocks, each of which is composed of masked self-attention layer and fully connected neural network. Each self-attention layer returns vector of size 256 that is taken as input by the fully connected network. The hidden layer of the neural network outputs vector of size 1024 and uses a GELU activation, and the final layer again returns vector of size 256 to be 7 used as input for the next decoder block. LigGPT consists of 8 such decoder blocks. LigGPT has around 6M parameters.
GCPN (Graph Convolutional Policy Network) You al.) 2018 leveraged graph con- volutional network and policy gradient to optimize the reward function that incorporates target molecular properties and adversarial loss. In each step, the allowable action to the current molecule could be either connecting a new substructure or an atom with an existing molecular graph or adding bond to connect existing atoms. GCPN predicts the actions and is trained via proximal policy optimization (PPO) to optimize an accumulative reward, including molecular property objectives and adversarial loss. Both policy network and adversarial network (discriminative training) use the same neural architecture, which is lhree-layer graph convolutional network (GCN) (Kipf Welling 2016) with 64 hidden nodes. Batch normalization is adopted after each layer, and sum-pooling is used as the aggregation function. Adam optimizer is used with le-3 initial learning rate, and batch size is 32.
MolDQN (Molecule Deep Q-Networks) (Zhou et al. 2019), same GCPN, formulate the molecule generation procedure as Markov Decision Process (MDP) and use Deep Q-Network solve it. The reward includes target property and similarity constraints.
16

---

Published as a conference paper at ICLR 2022
Following the original paper, the episode number is 5,000, maximal step in each episode is 40. Each step calls oracle once; thus, 200K oracle calls are needed in one generation process. The discount factor is 0.9. Deep Q-network is multilayer perceptron (MLP) whose hidden dimensions are 1024, 512, 128, 32, respectively. The input of the Q-network is the concatenation of the molecule feature (2048-bit Morgan fingerprint, with radius of 3) and the number of left steps. Adam is used as an optimizer with le-4 as the initial learning rate. Only rings with a size of and 6 are allowed. It leverages e-greedy together with randomized value functions (bootstrapped-DON) as an exploration policy, c is annealed from 1 to 0.01 a piecewise linear way.
GA+D (Genetic Algorithm with Discriminator network) (Nigam et al. 2020) uses deep neural network as a discriminator to enhance exploration in a genetic algorithm. / is an important hyperparameter that weights the importance of the discriminator's loss in the overall fitness function, and we set to 10. The generator runs 100 generations with population size of 100 for de novo molecular optimization and 50 generations with population size of 50 for molecular modification. Following the original paper (Nigam et al. 2020). the architecture of the discriminator is two-layer fully connected neural network with ReLU activation and a sigmoid output layer. The hidden size is 100, while the size of the output layer 1. The input feature is a vector of chemical and geometrical properties characterizing the molecules. We used Adam optimizer with 1e-3 as the initial learning rate.
RationaleRL (Tin et al. 2020) a deep generative model that grows a molecule atom-by- atom from an initial rationale (subgraph). The architecture of the generator is a message- passing network (MPN) followed by MLPs applied in breadth-firs order. The generator is pre-trained on general molecules combined with an encoder and then fine-tuned to maximize the policy gradient reward function. The encoder and decoder MPNs both have hidden dimensions of 400. The dimension of the latent variable is 20. Adam optimizer is used on both pre-training and fine-tuning with initial learning rates of 1e-3, Se-4, respectively. The annealing rate is 0.9. We pre-trained the model with 20 epochs.
MARS (Xie et al. 202T) leverage Markov chain Monte Carlo sampling (MCMC) on molecules with an annealing scheme and an adaptive proposal. The proposal is parameterized by a graph neural network trained on MCMC samples. We follow most of the settings in the original paper. The message passing network has six layers, where the node embedding size is set to 64. Adam is used as optimizer with 3e-4 initial learning rate. To generate basic unit, top-1000 frequent fragments are drawn from ZINC database (Sterling & Irwin 2013) by enumerating single bonds to break. During the annealing process, the temperature T 0.95[4/5] would gradually decrease to O.
ChemBO (chemical Bayesian optimization) (Korovina et al., 2020 leverage Bayesian optimization. It also explores the synthesis graph in sample-efficie
way and produces synthesizable candidates. Following the default setting in the original paper, the number of steps of acquisition optimization is set to 20. The initial pool size is set to 20, while the maximal pool size is 1000. Regarding the kernel, we leveraged the kernel proposed in ChemBO. which is the catimaal-trnssoot-bbaea graphical information explicitly. We did extra experiments to validate that it empirically outperforms the Tanimoto similarity kernel based on the molecular fingerprint. As for the acquisition, we followed the original paper and adopted the ensemble method using the EI (Expected Improvement), UCB (Upper Confidence Bound), and TTEI (Top-Two Expected Improvement) acquisitions: we conducted an extra experiment and found it empirically outperforms a single acquisition in our scenario, which is consistent to the ChemBO paper. GPs do not scale well with big data. To address this issue, use Subset of Date et al. (203)), which uses a subset of all the data points to reduce the size of the kernel matrix and learn the surrogate model efficiently. Following Chalupka et al. (2013), we selected subset of data points using k-means clustering based on molecular fingerprint and sampled data points from each cluster. The size of the subset is set to m 1000. The computational complexity is O(m²).
distance and kernel that accounts for
BOSS (Bayesian Optimization over String Space) (Moss et al. 2020) builds a Gaussian process surrogate model based on Sub-sequence String Kernel (SSK which naturally sup- ports SMILES strings with variable length, and maximizing acquisition function efficiently
based SMILES string
17

---

Published as a conference paper at ICLR 2022
for spaces with syntactical constraints. Following their default setting, we only keep the SMILES string whose length less than 50 (except the LogP task). At each BO step, BOSS samples 100 candidates, querying those that maximize the acquisition function predicted by SSK. It utilized genetic algorithm (GA) to optimize acquisition function efficiently under syntactical constraints. We empirically compared GA based acquisition optimizer with random search based acquisition optimizer (provided in the BOSS code repository) to validate its superiority. GA requires intensive oracle calls and is leveraged in inner-loop maximization in BOSS, so we do not spend too much computational resources in a single GA step. The population size is set to 100, the generation (evolution) number is set to 100. These setups are nearly optimal, and tuning these hyperparameters didn't 't improve the performance.
DST-rand a DST variant, which randomly selects basic operations (EXPAND, REPLACE, SHRINK) and substructure from vocabulary). During each iteration, it randomly selects basic operations (EXPAND, SHRINK, REPLACE) and substructure from vocabulary to generate at most K₃ molecules, select at most top-Ky molecules with highest GNN surrogate predictions. Then we evaluate the K₂ candidates with the real oracle and select C molecules as starting points for the next iteration. For a fair comparison, C' is set to 10, same as DST. To get the near-optimal setup, we try different combinations of (K1, K2). and finally find that K = 1000, K₂ = 100 achieves the best performance in terms of the average of top-100 molecules' score. Other setups of DST-rand follow DST.
C IMPLEMENTATION DETAILS
C.I DATASET
We use ZINC 250K dataset, which contains around 250K druglike molecules extracted from the ZINC database (Sterling & Irwin 2015). We first cleaned the data by removing the molecules containing out-of-vocabulary substructure and having 195K molecules left.
Vocabulary S: set of substructure. The substructure is the basic building block in our method, including frequent atoms and rings. On the other hand, atom-wise molecule generation is tricky due to the existence of rings. To select the substructure set S, we break all the ZINC molecules into substructures (including single rings and single atoms), count their frequencies, and include the substructures whose frequencies are higher than 1000 into vocabulary set S.
The final vocabulary contains 82 substructures, including the frequent atoms like carbon atom, an oxygen atom, nitrogen atom, and frequent rings like benzene ring. The vocabulary set covers over 80% of molecules in ZINC databases. After removing the molecules that contain out-of vocabl-vocayulary substructure, we use random subset of the remaining molecules train the GNNs, depending on the oracle budget.
The vocabulary size big enough for this proof-of-coneep: study. Other works also need to constrain their design space, such as MolDQN only allowing three types of atoms in a generation: "C", "N", "O" (Zhou al. 2019h: JTVAE (Jin et al., 2018). and RationaleRL (Jin al. 2020) only using frequent substructures similar to our setting. On the other hand, we may not want infrequent atoms or substructures because rare substructures ZINC may have some undesired properties such as toxicity, may not be stable, may not be easily synthesizable (Gao & Coley 2020). Also, rare substructures may impede the learning of oracle GNN. Note that users can enlarge the substructure space when they apply our method. We show all the 82 substructures in S in Figure 6
C.2 SOFTWARE/HARDWARE CONFIGURATION
We implemented DST using Pytorch 1.7.0, Python 3.7, RDKit v2020.09.1.0 on an Intel Xeon E5-2690 machine with 256G RAM and 8 NVIDIA Pascal Titan GPUs.
C.3 TARGET MOLECULAR PROPERTIES
Target molecular properties include
18

---

Published as a conference paper at ICLR 2022
QED represents quantitative estimate of drug-likeneess QED score ranges from 0 to 1. It can be evaluated by the RDKit package (https: www..wdki..or..... org/).
LogP represents octanol-water partition coefficient, measuring molecules' solubility, LogP score ranges from -os to +00. Thus, when optimizing LogP individually, we use the GNN model to do regression.
SA (Synthetic Accessibilitty score measures how hard it to synthesize a given molecule, based on a combination of the molecule's fragments contributions (Frtl & Schuffenhauer 2009h. It is evaluated via RDKit (Landrum et al.) 2006). The raw SA score ranges from to 10. A higher SA score means the molecule is hard to be synthesized and is not desirable. In the multiple-objective optimization, we normalize the SA score to [0. so that a higher normalized SA value is easy to synthesize. Following (Gao & Coley 2020), use the normalize function for raw SA score,
1,
SA(X) < P SA(X) 2
normalized-SA
exp (- (SA(X)-H)) )
20
where = 2.230044, a 0.6526308.
JNK3 (c-Jun N-terminal Kinases-3) belongs to the mitogen-activated protein kinase family and are responsive stress stimuli. such as cytokines, ultraviolet irradiation, heat shock, and osmotic shock. Similar to GSK33, JNK3 is also evaluated by well-trained random forest classifiers using ECFP6 fingerprints using ExCAPE-DE dataset (Li al. 2018 Jin
et al. 2020), and the range is also [0, 1].
GSK33 (Glycogen synthase kinase 3 beta) is an enzyme that in humans is encoded by the GSK3,3 gene. Abnormal regulation and expression of GSK333 is associated with an
increased susceptibility towards bipolar disorder. It is evaluated by well-trained" random forest classifiers using ECFP6 fingerprints using ExCAPE-DE dataset (Li al. 2018 i
et al. 2020). GSK38 score of a molecule ranges from 0 to
For QED, LogP, normalized SA, JNK3, and GSK33, higher scores are more desirable under our experimental setting.
C.4 EVALUATION METRICS
We leverage the following evaluation metrics measure the optimization performance:
Novelty is the fraction of the generated molecules that do not appear in the training set.
Diversity of generated molecules is defined as the average pairwise Tanimoto distance between the Morgan fingerprints (You et al. 2018 Jin al. 2020 Xie et al. 2021).
diversity = 1 | 1)
Σ
sim Z1, Z2).
(13)
24, 7,7262,,7177
where Z is the set of generated molecules. sim( Z2) is the Tanimoto similarity between molecule Z1 and Zz. (Tanimoto) Similarity measures the similarity between the input molecule and generated molecules. It defined as
b by
sim [X, Y)
||b lla b
bx is the binary Morgan fingerprint vector for the molecule X. In this paper, it is 2048-bit binary vector.
SR (Success Rate) the percentage of the generated molecules that satisfy the property constraint measured by objective f defined in Equation (1). For single-objective de novo molecular generation, the objective f is the property score, the constraints for JNK3,
"The test AUROC score is 0.86 (Jin et 2020
'The test AUROC score is also 0.86 (Jn al. 2020
19

---

Published as a conference paper at ICLR 2022
GSK33 and LogP are JNK32 0.5, GSK332 0.5 and LogP= 5.0 respectively. For multi- objective de novo molecular generation. the objective f is the average of all the normalized target property scores. Concretely, when optimizing "JNK3+GSK33", both JNK3 and GSK333 ranges from 0 to 1, is average of JNK3 and GSK333 scores; when optimizing first normalized SA 0 to 1. is average of QED, normalized SA, JNK3 and GSK3B scores. The constraint is the f score is greater than 0,4.
"QED-SA+INK3+GSKK"
# of oracle calls during the generation process. DST needs to call oracle in labeling data for GNN and DST based de novo generation, thus we show the costs for both steps.
chemical validities. As we only enumerate valid chemical structures during the recovery from scaffolding trees (Section C.3p, the chemical validities of the molecules produced by DST are always 100%.
N
CI
Br
Figure 6: All the substructures in the vocabulary set S, drawn from ZINC 250K database (Sterling &
Irwin 2015). includes atoms and single rings appearing more than 1000 times in the ZINC250K database.
C.5 ASSEMBLING MOLECULE FROM SCAFFOLDING TREE
Each scaffolding tree corresponds to multiple molecules due to rings' multiple combination ways. For each scaffolding tree, we enumerate all the possible molecules following Jin et al.) (2018) for further selection. We provide two examples in Figure illustrate it. Two examples are related to ring-atom combination and ring-ring combination, respectively. For ring-ring combination, our current setting does not support the spiro compounds (contains rings sharing one atom but no bonds) phenalene- like compounds (contains three rings sharing one atom, and each two of them sharing a bond). These two cases relatively rare chemical structures in the context of drug discovery (Supsana al.) 2005). As we only enumerate valid chemical structures during the recovery from scaffolding trees, the chemical validities are always 100%.
C.6 DETAILS GNN LEARNING AND DST OPTIMIZATION
Both the size of substructure embedding and hidden size of GCN (GNN) in Eq. 6 are d = 100. The depth of GNN L 3. When training GNN, the training epoch number is 5, and we evaluate the loss function on the validation set every 20K data passes. When the validation loss would not decrease, we terminate the training process. When optimizing "JNK3", "GSK3Ø3", "QED". "INK3+GSK3,3" and NOED=SA++NKK3GGSSS3
we use binary cross entropy as loss criterion. When optimizing "LogP", since LogP ranges from -oo to +o,, we leverage GNN to conduct regression tasks, and use means square emror (MSE) as loss criteria C. In the de
generation, in each generation, we keep C = 10 molecules for the next iteration. In most cases in the experiment, the size of the neighborhood set (Definition. [T) is less than 100. We use Adam optimizer with 1e-3 learning rate in training and inference procedure, optimizing the GNN and differentiable scaffolding tree, respectively. When optimizing DST(EEqation|| 10 method processes DST time, conduct 1000 Adam steps
20

---

Published as a conference paper at ICLR 2022
Molecule
EXPAND
CI
Scaffolding Tree
(a) Ring-atom connection. When connecting atom and ring in molecule, atom can be connected to any possible atoms in the ring. the example, there are 4 possible ways to add a Chlorine atom ("CI") as an expansion node to the target ring, a leaf node in the scaffolding tree.
800
Malecule
EXPAND
Scaffolding
(b) Ring-ring connection. When connecting ring and ring, there are two general ways, (1) one to use bond (single, double, triple) to connect the atoms in the two rings. (2) another two rings share two atoms and bond. In the example, there 14 possible ways add Cyclohexane ring (SMILES is "CICCCCC1") connect the target ring, leaf node the scaffolding tree.
Figure 7: Assemble examples.
in every iteration, which is sufficient to converge in almost all the cases. As complete generation algorithm, we optimize batch parallelly and select candidates based on DPP. When we use up oracle budgets, we stop it. All DST results in the tables take at most T 50 iterations.
C.7 RESULTS DIFFERENT RANDOM SEEDS
This section presents the empirical results that use different random seeds for multiple runs. In our pipeline, the random comes from in two steps: (1) Training oracle GNN: data selection/split, training process including data shuffle and GNN's parameter initialization. (2) Inference (Optimizing DST): before optimizing DST, we initialize the learnable parameter randomly, including N, w, A, which also brings randomness. To measure the robustness of the proposed method, we use 5 different random seeds for the whole pipeline and compare the difference of 5 independent trials. The results are reported in Table 4 find that almost all the metrics would not change significantly among various trials, validating the robustness of the proposed method.
Table 4: Results of 5 independent trials using different random seeds. For novelty, diversity, and SR (success rate), report the average value of 5 runs and their standard deviation.
Novdly!
OSK3
0.738:0.047
10K-5K
Log!
100.15 1.0.0%
07161:0002
100.0% 1.00%%
10K-5K
JNKIHGSK38
98.69 1.15
0.721
91.3%
10X-55K
ORD+SA-INK 3408K3
09.2% 9%
0'731200199
79.499 1.2%
20K-1KK
21

---

Published as a conference paper at ICLR 2022
C.8 COMPLEXITY ANALYSIS
We did computational analysis in terms of oracle calls and computational complexity. (1) oracle calls. DST requires O(TM) oracle calls, where is the number of iterations (Alg T M is the number of generated molecules (Equation. 12). scaffolding tree, for small molecule, K is very small. is the number of enumerated candidates in each node. As shown Figure 7 is also upper-bounded (J - 4 14 for the example in Figure 7 (2) computational complexity. The computational complexity is O(TMMC²) (the main bottleneck is DPP method, Algorithm 2). where the size of selected molecules C 10 for all the tasks (Section 3.4 & C. For all the tasks In Table and [
have M < KJ, K is the number of nodes in the
DST can be finished in 12 hours on an Intel Xeon E5-2699
562 machine with 256G RAM and NVIDIA Pascal Titan GPUs. The complexity and runtime are acceptable for molecule optimization.
D ADDITIONAL EXPERIMENTAL RESULTS
This section presents the additional empirical results, including additional results on de novo genera- tion, ablation study, chemical space visualization, and interpretability analysis (case study).
D.1 ADDITIONAL RESULTS de novo MOLECULAR GENERATION
This section presents some additional results of de novo molecular generation for completeness.
First, we present the optimization curve for all the optimization tasks in Figure g We observe that our method can efficiently reach high objective value within 10 iterations in all the optimization tasks. Worth mentioning that when optimizing LogP, the model successfully learned to add six- mix-member ring each step, as shown in Figure 9 and the objective (F) value grows linearly as a function of iteration number, which is theoretically the optimal strategy under our setting. Then, in Figure 10 we show the molecules with the highest objective (F) scores generated by proposed method on optimizing QED and "OED+SA+IN3+GGS337 Then compare method with baseline methods on 3 molecules with the highest objective (F) scores Table and for single-objective and multi-objective generation, respectively.
Average objective value
Average ofjective value
Average abjectiv value
Neration
Alration
Runalion
(a) LogP
(b) JNK3
(c) GSK33
Averrage objeective valure
Drver agge objeettive Live vvelure
Auerage abjectivn
Ton Tenation
Auration
Amation
(d) QED
(e) JNK3+GSK3B
(I) QED+SA+INK3+GSS333
Figure 8: The optimization curves in de novo optimization experiments. The objective value (F) is function of iterations.
22

---

Published as a conference paper at ICLR 2022
Logi -002277777777
Logh 20040
1025-442516 9479601119
Logi 19901/999005
Figure 9: The first eight steps in the de novo optimization procedure of LogP. The model successfully learned to add a six-member ring to each step.
0.60 2.25
77 1.67
0. 72 2.56
0. 75 2.50
0.68 2.01
0.93 0.83
0. 8.85
0.81 0. 83
78 0.79
0.69 0. 83
(a) Molecules with highest average QED, normalized -SA, JNK3 and GSK333 scores, four scores represent QED, raw SA, JNK3, and GSK33 scores, respectively.
0.947
946
0.946
0.946
0.946
(b) Molecules with highest QED.
Figure Sampled molecules with the highest scores.
Table 5: Highest scores generated molecules on single-objective de novo molecular generation. We present the result of DST in the first 50 iterations, but please note the setting of generation varies among the models, and an entirely fair comparison impossible.
JNK3
GSK3/3
LogP 2nd
Method
1st
2nd
3rd
1st
2nd
3rd
1st
3rd
GCPN
0.57
0.56
0.57
0.56
0.56
8.0
7.8
MolDQN
0.64
0.63
0.63
0.54
0.53
0.53
11.8
11.8
11.8
GA+D
0.81
0.80
0.80
0.79
0.79
0.78
20.5
20.4
20.2
RationaleRL
0.90
0.90
0.90
0.93
0.92
0.92
MARS
0.92
0.91
0.90
0.95
0.93
0.92
45.0
44.3
DST
0.97
0.97
0.97
0.95
0.95
0.95
49.1
49.1
49.1
D.2 De novo MOLECULAR OPTIMIZATION ON QED (POTENTIAL LIMITATION OF DST)
As we have touched in Section the optimization on QED is not as satisfactory as other oracles. We compare the performance various methods on single-objective de novo molecular generation for optimizing QED score and show the result in Table 7 Additional baseline methods include JTVAE (junction tree variational autoencoder) [Jin et al. 2018) and GraphAF (Graph Flow-based Autoregressive Model) (Shi et al.] 2020). The main reason behind this result is that our GNN predicts the target property based on a scaffolding tree instead of a molecular graph, as shown in Equation 7 A scaffolding tree omits rings' assembling information,
shown in Figure 2 Compared with other properties like LogP, JNK3, GSK33, drug- likeneiss more sensitive to how substructures connect (Bickerton et al. 2012). This behavior impedes the training of GNN and leads to the failure
23

---

Published as a conference paper at ICLR 2022
Table 6: Highest scores of generated molecules on multi-objectivv de novo molecular generation. The score is the average value of all objectives.
JNK3+GS33
QED+SA+JNK3+GSK33P
Method
top-1
top-2
top-3
top-1
lop-2
top-3
GCPN
0.31
0.31
0.30
0.57
0.56
0.56
MolDQN
0.46
0.45
0.45
0.45
0.45
0.44
GA+D
0.68
0.68
0.67
0.71
0.70
0.70
RationaleRL
0.81
0.81
0.81
0.76
0.76
0.75
MARS
0.78
0.78
0.77
0.72
0.72
0.72
DST
0.89
0.89
0.89
0.83
0.83
0.83
of optimization. We report the learning curve in Figure 11] where we plot the normalized loss the validation set as a function of epoch numbers when learning GNN. For fairness comparison, validation loss is normalized by dividing the validation loss at scratch (i.e., O-th epoch) that all the validation losses are between 0 and 1. For most of the target properties, the normalized loss value on the validation set would decrease significantly, and GNN can learn these properties well, except QED. verifies the failure of training the GNN on optimizing QED. differentiable molecular graph at atom-wise resolution may potentially solve this problem.
Table 7: Comparison of different methods on optimizing QED for single-objective de nove molecular generation. The results for baseline methods are copied from You et al. (2018); Zhou et (2019); Shi et al. (2020); Xie al. (202I). The results of JTVAE copies from
You al. (2018).
Method
top-1
top-2
top-3
TIVAR al. 2018
0.910
GCPN Your sh. 2018
0.948
0.947
0.945
MoIDON
GraphAP Shi er
2020
0.948
0.948
0.947
MARS Xie
0.948
0.948
0,948
DST
0.947
0.946
0.946
normalized validation loss
0.8
0.6
LogP
JNK3
GSK3B
0.4
0.2
epoch
Figure 11: Normalized validation loss-epoch learning curves. For fairness of comparison, validation loss is normalized by dividing the validation loss at scratch (i.e., O-th epoch) so that all the validation losses are between and For most of the target properties, the normalized loss value on the validation set would decrease significantly, and GNN can learn these properties well, except QED. The key reason for the failure of GNN on optimizing QED is the limitation of the expressive power of the scaffolding tree itself. QED is property that is highly dependent how substructures connect (Bickerton et al.) 2012). while our scaffolding tree currently ignores that information. See Section D.2 for more details and analysis.
D.3 RESULTS ANALYSIS FOR DISTRIBUTION LEARNING METHODS (LIGGPT)
As shown in Table 2 and T distribution learning methods (LigGPT) (Bagal et al. 2021) have much weaker optimization ability. DST and all the other baselines fall into the category of dire-didected
24

---

Published as a conference paper at ICLR 2022
molecule generation, a.k.a., molecule optimization, which generates molecules with high scores for a given oracle. In contrast, LigGPT belongs to distribution learning (a different method category), which learns the distribution of the training set. We refer to Brown et al. (2019) for more description of two categories of methods. Consequently, conditioned generation learns from the training set, cannot generate molecules with property largely beyond the training distribution, and can not optimize a property directly, even though they claim to solve the same problem. Problem formulation of distribution learning methods leads inability to generate molecules with property largely beyond the training set distribution, which means they are much weaker in optimization.
D.4 ABLATION STUDY
As described in Section during molecule sampling, we sample the new molecule from the differentiable scaffolding tree (Equation[TT] To verify the effectiveness of strategy, compare with a random-walk sampler, where the topological edition (i.e., expand, shrink or unchanged) and substructure are both selected randomly. We consider the following variants:
"DST DPP Both topology and substructure to fill are sampled from the optimized differentiable scaffolding tree, as shown in Equation (TT). This is what we use this paper.
"random DPP". Changing topology randomly, that is, each leaf node, "expand", "shrink" and "unchange" probabilities are fixed to 0.5, 0, 0.4. Substructure selection is sampled from the substructures' distribution in the optimized differentiable scaffolding tree. Then uses DPP (Section 3.4 to select diverse and desirable molecules for the next iteration.
"DST top-K". Same as "DST + DPP". it uses DST to sample new molecules. The difference is when selecting molecules for the next iteration, it selects the top- K molecules with highest score. It is equivalent to X too in Equation (12).
We show the results in Figure 12 We find that both DST sampling and DPP-based diversification play critical role performance. We check the results for "DST + top-K", during some period, the objective does not grow, we find it is trapped into local minimum, impeding its performance, especially convergence efficiency. 'random+DPP" exhibits the andom-walk behavior and it would reach satisfactory performance. When optimizing LogP, "DST +DPP" and "DST +top-K" achieved similar performance because logP score will prefer larger molecules with more carbon atoms, which less sensitive to the diversity and relatively easier to optimize. Overall, "DST +DPP" is the best strategy compared with other variants.
DST+DFP
DST-DFP
0.8
Objective valuee
DST+lopp- DPnomm+DPP
DST-log-K
20
random=DFP
0.6
0.6
25
Cojective
0.4
Objective
10.
0.2
DST+DFP
0.2
US14top-K
0.0
rancom+DPP
Iterations
iterations
iterations
(a) JNK3
(b) GSK3B
(c) LogP
Figure 12: Ablation study. Objective value (F) as function of iterations. See Section D.4 for more details.
D.5 CHEMICAL SPACE VISUALIZATION
We use principal component analysis (PCA) to visualize the distribution of explored chemical struc- tures in optimizing "JNK3 & GSK3,3'. Specifically, we fit two-dimensional principal component analysis (PCA) (Bro & Smilde 2014 for 2048-bil bit Morgan fingerprint vectors of 20K ZINC molecules, which randomly selected from the ZINC database (Sterling & Irwin 2015). Then to project the fingerprint of the generated molecule from various generations into wo--imensioona vector to observe their trajectories. The results are reported in Figure 13 where the grey points
the PCA
25

---

Published as a conference paper at ICLR 2022
represent the two-dimensional vector of ZINC molecules. We find that method explores different parts of the 2D projection the chemical space and covers a similar chemical space as the ZINC database after 20 to 30 iterations.
(a) After iterations
(b) After iterations
(c) After 10 iterations
(d) After 20 iterations
(c) After 30 iterations
Figure 13: Visualization chemical space covered during optimization. We used PCA to reduce the dimension of Morgan's fingerprint. The gray points the ZINC 250k data set. while colored points are generated molecules after corresponding iterations.
D.6 ADDITIONAL. INTERPRETABILITY ANALYSIS
leaf weight 0.99
expand weight
0.99
0.42
0.41
0.10
0.02
Figure 14: Interpretability analysis when optimizing LogP.
We provide interpretability example in Figure 14 At the leaf node (yellow), from the optimized differentiable scaffolding tree, we find that the leaf weight and expand weight are both 0.99. Thus we decide to EXPAND, the six-member ring is selected and filled in the expansion node (blue). This is consistent with our intuition that the logP score will prefer larger molecules with more carbon atoms.
E THEORETICAL ANALYSIS
In this section, we present some theoretical results of the proposed method. First, in Section E.I we conduct convergence analysis of DST under certain mild assumptions.
26

---

Published as a conference paper at ICLR 2022
E.1 CONVERGENCE ANALYSIS
In this section, we discuss the theoretical properties of DST in the context of de nove molecule design (learning from scratch). We restrict attention special variant of DST, named DST-greedy; the t-th iteration, given scaffolding tree Z(),, DST-greedy pick up only objective value from 2 neighborhood set set (2(()) i.e. Z(++1) arg max maxzeN ZEN(2(t) F(Z(() is exactly solved. We theoretically guarantee the quality of the solution produced by DS Dgg-redd,., First, make some assumptions and explain why these assumptions hold.
molecule with highest
Assumption 1 (Molecule Size Bound), The sizes (i.e., number of substructures) of all the scaffolding trees generated by DST are bound by Nain and Nmav
We focus on small molecule optimization; the target molecular properties would decrease significantly when the molecule size is too large, e.g., QED (enugllikeness) (Bickerton et al. 2012). Thus it is reasonable to bound the size of the scaffolding tree. addition, we use submodulantty and smoothness to characterize the geometry of the objective landscape.
Assumption (Submodularity and Smoothness). Suppose X1, X2. X3 are generated successively by DST-greedy via growing (i.e., EXPAND) a substructure on the corresponding scaffolding tree. We assume that the objective gain (i.e., AF) brought by adding a single substructure would not increase as the molecule grows (EXPAND).
F(Xs) - F S F(X2) F(X1),
(submodularity)
(14)
where X2 EXEAND $1), X3 EXPAND $2). 81,92 substructures add. Submodu- larity plays the role of concavity/corvexity in the discrete regime. On the other hand, we specify the smoothness of the objective function F by assuming
F(X3) F(X2) 2 (F((2)) F(XI)), 0 D<7<I (smoothness)
holds for the X1, X2, X3 described above, whose sizes are smaller than Nmin-
Then we theoretically guarantee the quality of the solution under these assumptions.
Theorem 1. Suppose Assumption Mand" hold; we have the following relative improvement bound with the optimum
1 y (F(X.)-F(XX)
( (1- 7)
(15)
where Z, is the local optimum found by DST-greedy, X, is the ideal optimal molecule, Xo is an empty molecule, starting point of de novo molecule design. In molecule generation setting, molecule is local optimum when its objective value is maximal within its neighbor molecule set, i.e. F(Z.)= 2 F(Z) for V Z E (Z..).
The proof is given in Section
Proof Sketch. We first show that DST-greedy can converge to local optimum within a finite step in Lemma[]] Then we decompose the successive generation path and leverage the geometric information of the objective landscape to analyze the quality of the local optimum.
Lemma (Local optimum). DST-greedy would converge to local optimum within finite steps.
The proof is given in Section
G.3
F EXTENSION OF MOLECULE DIVERSIFICATION
In the current iteration, we have generated molecules (',,... XM) and need to select C molecules for the next iteration. We expect these molecules to have desirable chemical properties (high F score) and simultaneously maintain higher structural diversity.
To quantify diversity, resort the determinantal point process (DPP) (Kulesza & Taskar 2012).. DPP models the repulsive correlation between data points (Kulesza & Taskar 2012)) and has been successfully applied many applications such text summarization (Cho al. 2019). mini-batch
27

---

Published as a conference paper at ICLR 2022
sampling (Zhang al. 2017), and recommendation system (Chen et al. 2018). Generally, have M data points, whose indexes are 2, ... .M. E R denotes the similarity kernel matrix between these data points, jith element of S measures the Tanimoto similarity between i-th and j-th molecules. We want to sample subset (denoted R. of M data, R is subset of {1,2, ... .M3 with fixed size C, it assigns the probability
P(R) & (et(SR), where R S (1,2, ,M, |R| =
(16)
where SR. CR is the sub-matrix of S, det(SR) is the determinant of the matrix SR. For instance, we want to sample a subset of size 2, i.e., R ji. then we have P(R) xx det(Sx)
S 5 SJ 1- S Syy more similarity between i-th and data points lower the probability of their co-oceurrence. DPP thus naturally diversifies the selected subset. DPP can be calculated efficiently using the following method.
Definition 11 (DPP-grest -opeggredyy (Chen et al. 2018). For any symmetric positive semidefinite (PSD) matrix S E RRMx M and fixing the size of R to C. Problem (16) can be solved in a greedy manner by DDP- gre-gre in polynomial time OCC2M) is denoted DeP-greeddd(( ), C?)
We describe the DPP-greedy algorithm in Algorithm 2 for completeness. Each iteration selects one data sample that maximizes the current objective, as described in Step 5 in Algorithm 2
Algorithm 2 DPP-greedy (Chen et al. 2018
1: Input: symmetric positive semi-definite matrix S E RMXxM, number selected data N
C'< M.
2: Output: R S (1.2. ...
.M.............................. .M. R| C.
3: W (1,2. ... .M.
4: for -ri== 1,2, 22..... ... ,.,,Col do
5: j arg max log deg kEW
6: R RU (s}.
7: W = W - -i}.
8: end for
As mentioned, our whole target is to select the molecules with desirable properties while maintaining the diversity between molecules. The objective is formulated as
arg max
(P)
Σ or + log P(S₂) = log det((e))))))))) + log (Sx),
(17)
(1,2, (M. IR-
SE
where the hyperparamter À O balances the two terms, the diagonal matrix V is
= =diag(expp(Anij ":xp))):"....(()"))))))
where U1 F(X,), .. UM = F(XM):
(18)
where n; is the F-score of the i-th molecule (Bq. [.. VR. is a sub-matrix of V indexed by R. For any square matrix M1, M2 of the same shape, have
det(M, = det (M1) det(MM)))))) = (et detM)),
we further transform CDpe (R) as below to construct symmetric matrix,
CDPP = log det VR.) lotd (Sx) = log det SR. = log det (v
(19)
where V = diag((expp) Av Ar=), ),...., exp( Aum ))). Then we present the following lemma for the usage of the DPP gre-dry method.
Lemma 2. Suppose S MM is the (Tanimoto) similarity kernal matrix of the M molecules.
b.b,
i.e. Sū |abl
by is the binary fingerprint vector for the i-th molecule, V is diagonal
matrix defined in Eg. 18 then we have (1) v is positive semidefiite;; (2) V Su V (ViSV
28

---

Published as a conference paper at ICLR 2022
The proof is given in Section
G.1
Thus, Problem 19 can be transformed as
LDPP(R) log det ((VAsV))
arg max R1 MI-RIIC
(20)
which means we can use DPP greedy (Def. to solve Problem
19 and obtain the optimal R.
Discussion. In Eq. (17). we have two terms specify the constraints the molecular property and structural diversity, respectively. When we only consider the first term (A ELER vt), the selection strategy is to select C molecules with the highest F score for the next iteration, same conventional evolutionary learning in Brown et al. (2019) : (2019); Nigam et al. (2020).
On the other hand, if only consider the second term in Eq. [7]. show the effect of selection strategies under certain approximations. molecules X1, X2, with high diversity among them, then leverage DST to optimize these C molecules respectively, and obtain C clusters of new molecules, i.e. ",,,..., ,,,, Ltd. DMG-Sammperr(N (x,)* A (x,) W (xx))i:... "..,,..., Zcve itd.
Suppose we have C
DMG-Sampler(N (xe)' A (Xe): W ((x))) Then we present the following lemma to show that when only considering diversity, under certain assumptions, Problem [17 reduces to multiple chain MCMC methods.
In Eq. I7 À is key hyperparamter, larger À corresponds to more weights on objective function F while smaller À specifies
diversity. When À goes to infinity, i.e., only considering the first term (A EIER Ut), it is equivalent to selecting C' molecule candidates with the highest F score for the next iteration, same conventional evolutionary learning in JJnssn (2019): Nigam al. (2020).
On the other hand, if we only consider the second term, we show the effect of selection strategies under certain approximations. Suppose we have C' molecules X1, X2, Xx with high diversity among them, then we leverage DST to optimize these C' molecules respectively, and obtain C' clusters of new molecules, i.c.
211, Zu, ~ DST-Sampler(N) N (x:) (Nxx)) A (xx) w(x)))
"01,... -oc iid DaT.Saaplee((((x))>>>>> A(x)) ((x)))
Then we present the following lemma to show that when only considering diversity, under certain assumptions, Problem (17) reduces to multiple independent Markov chains.
Lemma 3. Assume (1) the inter-cluster similarity is upper-bounded, i.e., sim Zip= 2n S C1 for any * j; (2) the intra-cluster similarity lower-bounded, i.e., sim( Zaps Figg) 2 1 62 for any
i E 11.2..... MM and # qi when both e1, approach to 0+, the optimal solution to Problem
(Zmp: ?zpo. Zcre..
where Pe = arg max, ((ep)) ) for c= 1,
, ,
The proof is given in Section
G.2
Remark. When the inter-cluster similarity low enough, and intra-cluster similarity high enough, molecule selection strategy reduces to multiple independent Markov chains. However, these assumptions are usually too restrictive for small molecules.
G PROOF OF THEORETICAL RESULTS
In this Section, we provide the proof of all the theoretical results Section and
G.1 PROOF OF LEMMA 2
Proof. (I) Visvi is positive semidefinite.
29

---

Published as a conference paper at ICLR 2022
First, let prove similarity kernel matrix S RBMxM based on molecular Tanimoto similarity is positive semidefinite (PSD), we know that the (i,j)-th element of S measures the Tanimoto similarity between i-th and j-th molecules, i.e. Siy
b bj
b 2 by
where b; E [0, 1] is the P-bit fingerprint vector for the i-th molecule (in this paper, P = 2048). S can be decomposed as
S
BB
where matrix B is the stack of all the normalized (divided by l2 norm, | . (la) fingerprint vector,
B b
b2
bp
E R PxM
|
|b2l22
|
For V x E RM,, we have
x
= = = (Bx) (Bx) 20 0.
Thus, S is PSD.
Then, similarly, for x RM,, we have
(BV v xT ((V)) " B BV (BV 2 0.
where V is diagonal matrix, so V = (V+)T. Thus, Visvi is symmetric and positive semidefi- nite.
(I)) S (VESVž)
Without loss of generalization, we assume R = "., ,c),, where 4 < 12 < ,., V is diagonal.
exp( Any 2
Atie
exp(
where
Vt. FX,,)
is the objective function of 4:-th molecule Xt.- The i, j-th element of V V is
v
Aun.
XUr
exp
SS,, exp
(21)
2
ij
Au On the other hand, the j-h element of V v is exp( of (v
2
Si exp( Ar
Then the i,j-th element
is
((V!VV)) =exp Auz. St,d exp
AU
(22)
2
Combining Equation
21 and 22 we prove VESR V2 = =V
G.2 PROOF OF LEMMA 3
Proof. We consider two cases in the solution R. (A) one molecule for each input molecule Zi,-. Ze. (B) other cases. Our solution belongs to Case (A).
(A) First, we prove for (A), our solution is optimal. We consider the second term in Equation (17).
SR is diagonal dominant. Also, determinant function is continuous function with regard to all
30

---

Published as a conference paper at ICLR 2022
the elements. Thus, det(SRR = I (Sx): goes to 1. Intuitively, all the selected molecules dissimilar to each other, and the diversity is maximized. On the other hand, to maximizing the first term in Equation (T). during each & 21, C), we select molecule with highest F score from (Zk1.... Zute J. That is our solution.
(B) Then we prove all the possible combinations in (B) are worse than our solution. In (B), based on the pigeonhole principle, there is at least one imput molecule Zy that corresponds at least two selected molecules. Without loss of generalization, we denoted them 2N, and Zkg- Since Sæ is diagonally dominant, its determinant can be decomposed as
(et = II (Sx)
k=1
there is at least one S₂ whose shape is greater than 1. Based on the definition of determinant, for matrix A E IR MMM
M II Amt)))) 1=1
det(A) = E sgn(n)
(23)
()Per((M)
where Perm(M) is the set of all permutations of the set 2.. ..,,M,,,, M). sgn(n) denotes the signature of 17, value that +1 whenever the reordering given by 17 can be achieved by successively inter- changing two entries an even number of times, and -1 whenever it can be achieved by an odd number of such interchanges. For exactly half of all r7s, sgn(n) and the other half are equal to -1. For the matrix A whose shape is greater than I and all the elements are equal to 1, the determinant is equal to Enpperm(M) sgn(1) = O.
The determinant function is a continuous function with regard to all the elements. When E2 goes to 04, all the elements of S, approach to I, the determinant goes to Thus, det(SR) also goes to The objective in Equation (17) goes negative infinity. Thus, it is worse than solution.
G.3 PROOF OF LEMMA 1
Proof. For the de novo design, DST-greedy start from scratch (empty molecule). First, we show in this setting, there is no "REPLACE" or "DELETE" by mathematical induction and contradiction. Since start from an empty molecule, at the 1-st step, the action is "EXPAND". Then we show the first / steps are "EXPAND", the (++1)-th step is still "EXPAND". Now we have (t) EXPAND(X st-1) (where 8:-1 is substructure). Suppose the (t + 1)-th step DST's action is "REPLACE", e.g.. (+++)) REPLACE(X st) (where s, is substructure), based definition of DST DST-greedy, have (++ ()) Since DST only REPLACE the leaf node, find that (++1) and (1) are both in neighbor molecule set of (t-1), i..., (t-1(,, which contradict with the fact that (1) arg max maxxxN((x=-1) (()( Similarly, we show that there would not exist "DELETE".
Then based on Assummption we find that DST-greedy converges at most Nmax steps.
G.4 PROOF OF THEOREM 1
Proof. Based on the Proof of Lemma I we find that there is only "EXPAND" action, then we can decompose the generation path as follows. Starting from scratch, i.c., Xu. suppose the path to optimum X, is
X, X1 + X2 + XR, X., where each step, one substructure added. The path produced by DST-greedy is
Zo(X,) Z₁ - Zn
Zkz Z..
Based on the definition of each step, only one substructure is added.
Based on Assumption [ we have Nmin & kj k2 S Nmux- There might some overlap within the first several steps, without loss of generalization, we assume 7 = X and Ze+1 # X where k can be 0, 1, 1........ A,. Based on Assumption we have F(X,) F(Xo) F(X2) F(X,) 2 F(X2)= F(X2)
F(Xx,) - F(XE,11)-
31

---

Published as a conference paper at ICLR 2022
Then, have
k1 (F((1) - F(Xo))
(((X))-F((X))) +
+ (F(Xs) - (F(X)) + (F(Xx)) F(Xx)
=F(Xk,)- -F(X))
Thus, we get
1 (F(x,))
F(XI) - F(Xo) 2
(24)
1 (F((.)
(F(Xx,) (0xx,)-((M))) F((o)) Nia
Since Zo Xo. according to the definition of greedy algorithm, we have F(Z1) 2 F(Xi). Based on Assumption 2
we have
F(ZN F(
2 7 (Z Natio-1 - -F(ZN..- 2 (F(NNNN - F(ZN..-aa)
2 ... 2 New (F((z) - -F(Zo)).
Based on Assumption 1
we have F(Z.) F(Zo) 2 F(ZNam) F(Zo). Then we have
F(Z)) - F(Zo)
FZ (Zoo)
(F(Z Nmn 1)) (F(ZN..-1) F(NN..22)) +.... (F((Z) F(Zo))
(25)
2 (1 + + 72 + 7......m....... Nnn (F(Z)) -F(Zo))
1 - 7 Nuin
(F(Z)) - F(Zo))
Combining Equation
(24) and 25). we have
F(Z.) F(Xo) 1 Nin
N (F(X.) - (F(X.)-F(Xx))
We observe that objective F's improvement relatively lower bounded.
32


Reviews:
[{"review_id": "t-evfeADi7A", "reviewer": "Reviewer_BjiD", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "I thought that this paper was well-written and had a lot of genuinely interesting ideas. It was interesting to read. My main concern is that I think it is possible that the proposed algorithm is overly complicated and doesn't work as well as the authors claim. I will provide more details below.\n\n== Pros ==\n\n- To my knowledge, the differential scaffolding tree (DST) idea is novel. I think it is a very interesting contribution\n- the details of the DST method are very well-described in section 3. I think it would be possible to re-implement their method based on the information given, which is a good sign\n- Choice of experiments seemed good (I'm glad this is not yet another paper maximizing logP). It was nice to see that you included many baselines\n- Good mentioning of relevant related work.\n- The interpretability of the method is really interesting!\n\n== Cons ==\n\n- Looking at algorithm 1 (the main contribution), it is clearly a genetic algorithm. There is a population, and a loop over \"generations\" which includes a \"mutation\" step and a \"new population selection\" step. I found it slightly odd that the authors did not point out this connection, or use any terminology related to genetic algorithms. Given the similarity, it would be nice to see the authors discuss how their GA compares to previously proposed GAs for molecule optimization (such as the GAs from Guacamol, or the various SELFIES GAs by Alan Aspuru-Guzik's group)\n- Thinking of the contribution as a genetic algorithm, it actually isn't clear to me that the whole DST thing is necessary. Fundamentally, the proposed mutation step takes a molecule X (in the form of a scaffold tree), and returns a set of candidate molecules whose scaffold tree differs from the current molecule's scaffold tree by 1 local editing operation (called $\\mathcal{N}(X)$ in the paper). The DST effectively provides a way of sampling from $\\mathcal{N}(X)$, where samples are (probably) predicted to have a high score by the GNN. While this is sensible, I can think of many similar alternatives which are much less complicated (see below). Overall I worry that the proposed method is just an overly complex way of doing something quite simple.\n  1. Just enumerate all of $\\mathcal{N}(X)$. My guess is that for small molecules this set is actually not that large ($\\le 10,000$) although I could be wrong about this.\n  2. Since the above approach may propose many poor molecules (and therefore not be sample efficient), one could enumerate $\\mathcal{N}(X)$ and return a subset which the GNN predicts will have high scores. In fact, any oracle could be used here (e.g. random forest, GP) which might be less prone to overfitting and work with less training data.\n  3. For computational efficiency, the approaches above could be modified to just randomly sample from $\\mathcal{N}(X)$ rather than enumerating it. This is the approach taken by most GAs.\n- Many details of the experimental section are not clear in the main text. As a reader I felt that the experiments offered little value without reading a lot of the appendix. Aware that there is a page limit, I nonetheless think that the following things needed to be in the main text:\n  1. Brief explanation of the objectives. Why should I care about them? Are they novel objectives, or were they proposed in previous work? (seems like the latter)\n  2. Why were your baselines chose, beyond just being other methods which work for molecular optimization?\n  3. What is success rate? (turns out this is the only metric related to actually succeeding at optimization)\n  4. Why should I care about novelty and diversity? Some readers might just care about finding $n$ molecules with high scores. I think there is a strong case for caring about these but it should at least be mentioned.\n  5. How was the number of oracle calls chosen?\n- Generally weak baselines w.r.t. sample efficiency: the authors claim that needing a small number of oracle calls is a key advantage of their method, and indeed, their method outperforms other methods with fewer oracle calls. Yet, as far as I can tell, none of the methods they compared with would actually be expected to succeed with a small number of oracle calls. GCPN, MolDQN, RationaleRL are all RL methods which only produce good molecules after a long training period. GA+D and MARS are exploratory methods which, as far as I understand, call the oracle indiscriminately on a lot of random points to explore molecular space. I am unfamiliar with LigGPT. To my knowledge, the most sensible baselines to compare with for sample efficiency would be Bayesian optimization. I suggest perhaps references [1-4] (in that order). If object that GPs don't scale and therefore are not applicable here, I would counter by saying that they will easily scale up to 10k points on a reasonable computer, and approximate GPs with 10k inducing points should be manageable with the 50k oracle calls shown in Figure 3.\n- In table 1-2, the number of oracle calls is listed as if it is a fundamental property of each algorithm, which seems strange because it isn't: the number of oracle calls is directly decided by the number of iterations, the number of calls per iteration, etc, which can be chosen by the user. Simply put, all of the methods listed can be run with any number of oracle calls (although with fewer oracle calls performance will almost certainly be worse). It isn't clear whether the number of oracle calls was simply chosen so that DST had the lowest number, whether the number of oracle calls is really just the number of oracle calls necessary to meet some sort of external criterion (e.g. a certain number of successful molecules produced), or something else. In general, I think it is more fair to show comparisons like in Figure 3, which plots performance as a function of the number of oracle calls.\n- Strange choice of metrics in Tables 1-2. I see many disadvantages to the metrics proposed in Tables 1-2. Generally, they are based around evaluating the entire set of molecules produced, rather than just evaluating the top points, which is more normal in optimization (since one generally cares about the optimum reached, not how one gets there). I particularly disagree with the use of \"success rate\" to measure optimization performance. Firstly, the cutoff for \"success\" seems fairly arbitrary. It is not clear how performance changes when this cutoff is changed. For this reason I think people usually report the top N scores or the average of the top N scores, rather than a measure of how many molecules exceed a threshold. Secondly, by measuring success as a \"rate\" the scores of sample inefficient methods seem to be unfairly penalized: even if they produce more molecules of better quality, the presence of many poor molecules (e.g. from the beginning of training in RL) will lower the score. Overall, I think that Figure 3 provides a much more valuable comparison of methods than Tables 1-2\n\n== Suggestions for Improvement ==\n\nOverall I think that the paper would benefit greatly from the following improvements:\n\n1. Some sort of ablation study or comparison with modified GAs to determine whether the added complexity of the DST is actually necessary to achieve their optimization performance.\n2. The addition of stronger baselines (I suggest specifically Bayesian optimization). This will help evaluate the performance of the proposed method.\n3. Re-organization of the experimental section, including more explanations and reporting better metrics.\n\n## References\n\n[1] Korovina, Ksenia, et al. \"Chembo: Bayesian optimization of small organic molecules with synthesizable recommendations.\" International Conference on Artificial Intelligence and Statistics. PMLR, 2020.\n\n[2] Moss, Henry, et al. \"BOSS: Bayesian Optimization over String Spaces.\" Advances in Neural Information Processing Systems 33 (2020).\n\n[3] Moss, Henry B., and Ryan-Rhys Griffiths. \"Gaussian process molecule property prediction with flowmo.\" arXiv preprint arXiv:2010.01118 (2020).\n\n[4] Gómez-Bombarelli, Rafael, et al. \"Automatic chemical design using a data-driven continuous representation of molecules.\" ACS central science 4.2 (2018): 268-276.", "summary_of_the_paper": "Although the term \"genetic algorithm\" is not used in the paper, this paper effectively proposes a fancy genetic algorithm (GA) for molecule optimization. The contribution is chiefly in the mutation step of the GA, which roughly uses a pre-trained graph neural network to take the derivative with respect to the adjacency and node identity matrices of a graph to produce a probability distribution over mutated graphs. They show results on several different molecular optimization tasks which suggests that their algorithm outperforms several baselines.", "summary_of_the_review": "I feel quite conflicted on this paper: on the one hand, the DST is to my knowledge novel and fairly interesting. On the other hand, the lack of \"good\" baselines makes it hard for me to judge the performance of this method, and I am concerned that a lot of the complexity of the DST is not actually necessary. Tentatively I think I will recommend rejection due to these concerns, but I would be happy to raise my score if they are addressed.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "byeclmTp_hr", "reviewer": "Reviewer_5uXz", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Pros:\n- The idea proposed by authors seems novel and interesting.\n- The results reported by authors are really promising.\n- The proposed method allows to use a continuous gradient-based optimization of oracle function. Thus there is no need to train a latent space model before optimization.\n\nCons:\n- In my opinion the section that describes the experimental settings  is not sufficiently developed. E.g. I cannot see how many molecules were generated and how many optimization iterations were performed for every molecule (the authors wrote that DST reached the highest scores within T = 50 iterations, but they do not specify whether they use this amount of steps in every experiments). \n\nQuestions during rebuttal period:\n- How many adam steps were conducted in every iteration? ", "summary_of_the_paper": "The paper proposed differentiable scaffolding tree (DST), a novel method for molecular generation, that allows for a continuous gradient-based optimization of selected molecular properties. The results of de novo molecular optimization are really good. DST outperformed other methods in terms of optimization success rate, at the same time the proposed molecules were diverse and novel. The method also required less oracle calls (and majority of them before starting the optimization), compared to the other methods.", "summary_of_the_review": "The method seems novel and interesting and results are promising. The only drawback is the weak section describing the experimental settings, which, however, can be easily corrected.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "XhxvgB_qLFt", "reviewer": "Reviewer_xvty", "summary": "", "questions": "", "limitations": "", "rating": 10, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "Difference to related works: Most of the existing deep generative models do not scale well on large datasets due to non-differentiable nature of molecule structures and they often heavily rely on brute-force enumeration. The proposed DST significantly increases scalability by employing implicit differentiation for gradient updates, i.e., fully differentiable generation of new molecules that preserves computation time and storage space without introducing any additional latency overhead or stochastic variation being induced by random initialization of weights in DST architecture. \n\nLimitations: More quantitative performance comparisons with related works are needed. \nThe authors should expand on the approximation for formulating the problem as an iterative local discrete search specifically as it relates to equation 10. \nWhy did the authors allow for single atoms to be utilized as \"substructures\"? This aligns against the intuition of chemists who usually consider atoms as sidechains for modification as opposed to substructures for clustering or analyzing data by series.", "summary_of_the_paper": "The molecular optimization is an essential chemical science and engineering task that has important applications such as drug discovery. Deep generative models and combinatorial optimization methods achieved initial success but still struggle with directly modeling discrete chemical structures and often heavily rely on brute-force enumeration. The challenge comes from the discrete and non-differentiable nature of molecule structures. To tackle this, the authors propose a differentiable scaffolding tree (DST) that utilizes a learned knowledge network to encode prior information. \n\nImportant contributions include:\n\n(1) Definition of a local derivative of a chemical graph, we propose the differentiable scaffolding tree. The theory allows for gradient-based optimization of a discrete graph structure.\n\n(2) The authors generalize a chemistry-aware graph network to encode prior knowledge. The network is capable of learning the hierarchy between functional groups and atoms, as well as the topological relations among chemical structure fragments.\n\n(3) A unified graph-based optimization framework is proposed that can tackle both single-task and multi-task objective functions.", "summary_of_the_review": "This paper created a molecular graph locally differentiable by using a differential scaffolding tree (DST). To the best of our knowledge, it does appear as a novel approach to make the molecular optimization problem differiable at the structure level rather than utilizing latent spaces or employing RL/evolutionary algorithms. The authors developed a generic molecular optimization approach based on DST, which was validated through comparison when previous methods and benchmark datasets. The paper appears well written and the results are promising.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "PLcnmyczNm-", "reviewer": "Reviewer_JXBo", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "The paper provides a few interesting ideas although it also contains some weaknesses. \n\nStrengths:\n1. **Differentiable scaffolding tree**: The DST is designed as a learnable graph input of molecules. By back propagating the parameters in A, N and w, it samples out new molecular structure along the direction of optimizing purposed properties. DST implies a novel direction of resolving the problem of lead optimization.\n2. **Scaffold tree instead of molecular graph**: The proposed method smartly uses scaffold tree instead of molecular graph to alleviate invalid molecules.\n3. **Determinantal point process**: The paper leverage DDP to keep expected properties as well as the diversity of generated molecules.\n\nWeaknesses:\n1. **restricted backbone**: The backbone of the oracle model (as GNN used in the paper) must be a differentiable model. This limits the usage of the proposed method. For example, if the oracle response comes from the real wet experiments, the proposed method is unable to get gradients from the wet experiments. Moreover, it may also suffer from the activity cliff of molecular properties.\n2. **Local optimum**: SInce the proposed method generates new molecules step by step, it is possible that the selected pathway of generating new molecule is not an optimal solution. When it falls into a local optimum, the properties of the generated molecules will not be improved more.\n3. **Costly assembling process from scaffold tree**: it is also a problem of JTVAE by practice. Assembling molecular graphs from scaffold tree is very time consuming.\n\nSuggestions:\n1. **Better assembling way of scaffold tree**: the authors could leverage some ideas from [1] to avoid costly assembling process.\n2. **Add constraints to learnable N, A, and w**: Some lead optimization tasks may require un-changeable portions of the molecules. This could be easily achieved by adding constraints to learnable N, A, and w. \n\n[1] Jin, Wengong, Regina Barzilay, and Tommi Jaakkola. \"Hierarchical generation of molecular graphs using structural motifs.\" International Conference on Machine Learning. PMLR, 2020.", "summary_of_the_paper": "The paper designs a novel model to optimize molecular leads with a differentiable scaffolding tree. The idea of learnable parameters in indicator matrix, adjacency matrix, and node weight vector is novel and interesting although the method still contains some weaknesses. The paper still provides a good direction of molecular generation.", "summary_of_the_review": "Overall, it is a good paper with some interesting designs like DST. The proposed method shows a new direction to do lead optimization. I prefer to accept this paper although it contains some cons that hard to avoid.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper was a tough call. The key contribution of the paper is a genuinely useful technique for generating chemical compounds satisfying desired properties. However, there are some key issues with paper.

Reviewer *BjiD* found out that baselines are weak. Most importantly, he run thorough experiments with GraphGA, outperforming by a significant margin the baselines. With minor tweaks (e.g. enabling generating larger molecules). GraphGA achieves comparable though slightly weaker results to DST. Importantly, as pointed out by Reviewer BjiD, there is an important flaw in the experiments that some methods have a cap on the number of atoms they can add. For example, on the logP optimization task, it is possible to optimize the score by just adding carbon atoms. I would like to thank very much Reviewer for going beyond and running these experiments.

All reviewers emphasized the novelty as a key contribution. In internal discussion, I raised concern about novelty and framing of the work. One could argue that any autoregressive model (i.e. adding atoms and bonds at each step) forms a DST. One could also argue that training LSTM to produce the distribution of interest, like in [1], is also a DST because the fine-tuned LSTM encodes the distribution of many molecules and is differentiable with respect to the distribution it encodes.

Despite these flaws, it is a solid contribution, which is likely to be useful for the community. Thank you for your submission, and it is my pleasure to recommend acceptance.

For the camera-ready please: (a) include a well-tuned GraphGA (implementing different tradeoffs of diversity and fitness), (b) include LSTM as implemented in Guacamole as baseline, (c) discuss much more clearly novelty of the work. Additionally please ensure that other baselines are not hampered by limit on number of atoms they can add.

References:

[1] Generating Focussed Molecule Libraries for Drug Discovery with Recurrent Neural Networks, Segler et al, [https://arxiv.org/abs/1701.01329](https://arxiv.org/abs/1701.01329)
Decision: Accept (Poster)