Molecule Edit Graph Attention Network: Modeling Chemical Reactions as Sequences of Graph Edits

The central challenge in automated synthesis planning is to be able to generate and predict outcomes of a diverse set of chemical reactions. In particular, in many cases, the most likely synthesis pathway cannot be applied due to additional constraints, which requires proposing alternative chemical reactions. With this in mind, we present Molecule Edit Graph Attention Network (MEGAN), an end-to-end encoder-decoder neural model. MEGAN is inspired by models that express a chemical reaction as a sequence of graph edits, akin to the arrow pushing formalism. We extend this model to retrosynthesis prediction (predicting substrates given the product of a chemical reaction) and scale it up to large datasets. We argue that representing the reaction as a sequence of edits enables MEGAN to efficiently explore the space of plausible chemical reactions, maintaining the flexibility of modeling the reaction in an end-to-end fashion, and achieving state-of-the-art accuracy in standard benchmarks. Code and trained models are made available online at https://github.com/molecule-one/megan.


Introduction
Synthesis planning answers the question of how to make a given molecule. Due to the substantial size and complexity of the reaction space, synthesis planning is a demanding task even for skilled chemists and remains an important roadblock in the drug discovery process 1 . Computer-aided synthesis planning (CASP) methods can speed-up the drug discovery process by assisting chemists in designing syntheses [2][3][4][5] . These methods are also increasingly employed in de novo drug design, in which estimating the synthetic accessibility of a large number of compounds has surfaced as a key challenge to the field 6 .
Designing a synthesis plan involves predicting reaction outcomes for a given set of possible substrates (forward synthesis prediction), or proposing reactions that can simplify a given target molecule (retrosynthesis prediction) 7 . In both cases, it is key for automated synthesis planning tools to model reliably outcomes of a diverse set of chemical reactions. This is especially important when there are additional constraints imposed on the desired synthesis plan such as avoiding certain starting materials or using green chemistry 8 .
Many approaches to predicting reaction outcomes employ a static library of reaction templates.
A reaction template encodes a graph transformation rule that can be used to generate the reaction based on the product or subtrates [9][10][11][12][13][14] . Due to computational limitations, they typically require representing reactions using a restricted number of templates, which necessarily limits the coverage of the chemical reaction space accessible by such methods.
The small coverage and other limitations of template-based methods motivated the development of more flexible deep learning models, which aim to learn transformations and their applicability rules directly from the data 15,16 . Despite the added flexibility and the resulting better coverage of the chemical space, they tend to propose a limited set of plausible chemical reactions for a given input 17 .
Models such as proposed in Schwaller et al. 15 , Karpov et al. 16 belong to the class of sequence to sequence models. They generate a chemical reaction by sequentially outputting individual symbols in the SMILES notation 18 . We hypothesize this makes it challenging to generate more diverse predictions from the model, as there is no natural decomposition of the predictive distribution into plausible but different reactions. Notably, many previous papers have pointed out limited diversity of sequence to sequence models in other domains [19][20][21] .
An arguably more natural idea is to express a reaction as a sequence of graph edits such as bond addition or removal, which is inspired by how chemists describe reactions using the arrow pushing mechanism 22,23 . Bradshaw et al. 22 introduce a model that predicts the reaction mechanism represented as a sequence of bond removal and additions that directly correspond to electron paths.
However, their approach is limited to a certain subset of the chemical reaction space (reactions with linear chain topology) and forward synthesis. In a similar spirit, Do et al. 23 models a reaction as a set of operation on atom pairs. However, due to the lack of support for atom addition, their method also cannot be readily applied to retrosynthesis prediction. We hypothesize that expressing reaction in this way enables to easily search through different alternative reactions because the proposed reaction is obvious usually after the first few generation steps. This allows the generation algorithm to naturally backtrack and explore other possible reactions. We note that related methods have a rich history in parallel applications of generative deep learning such as natural language processing 24 .
In this work, we propose the Molecule Edit Graph Attention Network (MEGAN), an encoderdecoder end-to-end model that generates a reaction as a sequence of graph edits. Designing a novel architecture and training procedure enables us to apply this approach to both retro and forward synthesis, and scale it up to large datasets. More specifically, our two main contributions are as follows: • We extend the idea to model the reaction as a sequence edits to retrosynthesis prediction and scale it to large datasets. This required introducing a novel encoder-decoder architecture and an efficient training procedure that avoids the need to use reinforcement learning.
• We achieve competitive performance on retrosynthesis and competitive performance on forward synthesis as well as state-of-the-art top-k accuracy for large K values on all tested datasets. This serves as evidence that MEGAN achieves excellent coverage of the reaction space.
We open-sourced our code and trained models at https://github.com/molecule-one/megan.

Related work
Computer-aided reaction prediction has a rich history 25 . Early approaches to generating reactions relied on manually crafted rules 26 , which were difficult to apply to novel chemistry. Another family of methods uses physical chemistry calculations 27 . These methods are typically too computationally intensive to be used in synthesis planning software and can be complemented using statistical methods that learn from data.
Statistical approaches to predicting reactions and designing synthesis paths can be broadly categorized into template-based and template-free approaches. Template-based approaches use reaction rules or templates. These templates can be automatically mined from a database of known reactions or defined manually [9][10][11][12][13][14] .
Manually encoding reaction templates necessitates developing complex rules governing their applicability 14,28 . This motivated the recent rise in popularity of machine-learning models used to select or rank most relevant templates 3,29 . Such models can be applied to assess the reactivity of atoms to which reaction rules should be applied 30,31 . An interesting alternative to this approach is to prioritize transformations that have already been applied to similar molecules 32 . A suitable set of templates combined with a strong ranking model can achieve state-of-the-art performance on standard retrosynthesis benchmarks. 13 The small coverage and other limitations of template-based methods motivated the development of template-free models. A particularly successful class of template-free methods is based on sequence to sequence models that sequentially (symbol by symbol) predict the target SMILES 18 string from the input SMILES string 18 . Such approaches can be used out-of-the-box for both forward synthesis 15,33 and retrosynthesis 16,34 . However, they have been shown to make some trivial mistakes 28 , and produce reactions with a limited diversity 17 . These models also act as black-boxes; they do not provide the reasoning behind their predictions and are not able to map atoms between the substrates and the product.
We hypothesize that the limited diversity of sequence to sequence models can be attributed to the fact that they generate reaction as a sequence of SMILES symbols from left to right. Many papers have also pointed out shortcomings in the diversity of the generated output by sequence to sequence models [19][20][21] . In natural language processing, there is a rich parallel history of exploring alternative generation ordering and adding through various means structure to the output. Gu et al. 24 show that using an alternative order of output symbols than left to right enables generating more diverse prediction. Mehri and Sigal 35 propose an alternative decoding algorithm to improve the diversity of sequence to sequence models that starts from generating the middle word. These modeling choices can be viewed as incorporating more suitable inductive biases for the given domain. From this perspective, we study an alternative representation of the output in generative models applied to reaction outcome prediction.
Our work focuses on a recently introduced class of template-free methods, which defines reaction generation as predicting target graph by sequentially modifying the input graph. Such an approach, while remaining template-free, can provide greater interpretability of predictions, which are modeled as direct transformations on molecules. Models that predict graph edits have been successfully applied to forward synthesis 23 and can be fine-tuned to provide especially interpretable results on certain subsets of the chemical reaction space 22 . In this work, we propose a model inspired by these prior works and extend the approach to retrosynthesis prediction.
Concurrently to our work, the graph-edit reaction generation approach has been employed to retrosynthesis in G2GS and GRAPHRETRO models [36][37][38] . The G2GS model presented in Shi et al. 36 , Yan et al. 38 consists of two separate modules for predicting the reaction center and generating final substrates from the disconnected synthons. GRAPHRETRO 37 employs a similar framework but completes the synthons with substructures which are selected from a predefined set found on the training data. In contrast to these concurrent works, we propose an end-to-end model for generating reaction through a sequence of edits and apply it both to retro and forward synthesis prediction.
Our work on reaction outcome prediction is motivated by their wide utility. An important use of reaction outcome prediction models arises in automated synthesis planning; see for example Segler et al. 3 for an excellent discussion. It is also worth noting that reaction outcome prediction has recently risen in importance in the context of de novo drug discovery. Gao and Coley 39 have drawn attention to the fact that many recent de novo drug discovery methods tend to generate hard to synthesize molecules. In this context, automated synthesis planning can be used to filter out these hard to synthesize compounds, and compares favorable to other simpler methods 40,41 . Other works use reaction outcome prediction models to define the space of accessible molecules to de novo models 42 .

Molecule Edit Graph Attention Network
Molecule Edit Graph Attention Network (MEGAN) is an encoder-decoder architecture based on the graph convolutional network architecture that predicts a sequence of actions on atoms and bonds of a chemical compound.

Input and output representation
We begin by describing the input and the output representation of the MEGAN model.  We add an additional node to the graph named supernode 44 , which is connected to all atoms in the graph with a special SUPERNODE bond type. In graph convolutional neural networks, supernode is particularly useful for passing information between connected components of a graph, especially after it was split by deleting a bond. Adding the supernode makes n equal to the number of atoms in the input molecule + 1.
To construct feature vectors for atoms and bonds , we use the following chemical properties calculated using the RdKit package 45 : ATOMIC NUMBER, FORMAL CHARGE, CHIRAL TAG, NUMBER OF EXPLICIT HYDROGEN ATOMS, IS AROMATIC for atoms and BOND TYPE, BOND STEREOMETRY (None, Z or E) for bonds. All these features are set to zero for the supernode and all bonds with the supernode, except for BOND TYPE which has a special supernode value for bonds between atoms and the supernode. We also add a binary flag to indicate whether a given atom or bond was edited or added in any of the previous steps. Finally, we indicate whether a given node is an atom or the supernode with a binary flag on each node. Final atom and bond vectors are constructed by concatenating one-hot vectors of all of their categorical features. We describe the remaining details in Section Featurization in the Supplement.
Output representation The output of MEGAN is a sequence of graph actions, which are then applied to the input graph to produce the desired output graph. An action can be applied to a single atom or a pair of atoms. In the latter case the action modifies the bond between the two atoms.
To train the model, we compare the outputted sequence of actions with the ground-truth sequence based on a mapping between the product and the substrates (see Section Gradient-based training of MEGAN). We define the following graph actions: • Edit atom properties (EditAtom) • Edit bond between two atoms (EditBond) • Add new atom to the graph (AddAtom) • Add new benzene ring to the graph (AddBenzene)

• Stop generation (Stop)
EditAtom changes properties of atoms, such as the formal charge, chirality, or aromaticity.
EditBond adds, edits, or deletes a bond between two atoms. AddAtom adds a new atom of a specified type as a neighbor of another atom already existing in the graph, with a specified bond type. AddBenzene simplifies generation by allowing the model to append a complete benzene ring to a selected carbon atom. Stop action indicates the end of the generation process and can be outputted only by the supernode. We describe the possible actions in detail in the "Graph edit actions" in Supplement. f

Model architecture
Encoder-decoder architecture of MEGAN The MEGAN model consists of two parts: the encoder, which is invoked only once per reaction generation at its beginning, and the decoder, which is sequentially invoked to generate actions. This general encoder-decoder architecture is shared by the most successful sequence to sequence models in deep learning Vaswani et al. 46 .
Both encoder and decoder are variants of graph convolutional networks (see GCN-att layer) that modify hidden vector representation of nodes. They take the same input and output node feature vector size n a and input edge feature vector size n b . The input graph x t is modified according to the action predicted by the decoder at each step t, until a special Stop action is outputted. We stack n e and n d layers of the encoder and decoder, respectively. The architecture and the inference process are depicted schematically in Figure 1.
Generating graph actions At each generation step t, the one-hot encoded features of the current A t = g emb (A OH t ) to convert atoms and bonds to n a and n b -dimensional spaces, respectively. We treat n a and n b as hyperparameters.
In the first generation step (t = 1), we pass the embedded graph through the encoder and the decoder to arrive at atom features H 1 = dec(enc(H 0 1 )). At each subsequent step t, we only use the decoder. We incorporate previous step representation by taking the element-wise maximum operation over the current and previous atom features. Namely, for t > 1, H t = dec(max(enc(H t−1 ), H t−1 )).
We zero-pad features of H t−1 for any node that was added to the graph at step t. The encoder and decoder also take as input A t .
The logits L a t and L b t for atom actions Act a and bond actions Act b are calculated at each step t from H t and A t as follows: , where g atom , g atom , g bond , and g bond are linear layers; σ r is the ReLU activation function. The GCN-att layer The basic building block of the encoder and decoder is an attention-based GCN layer that we call GCN-att. We enhance the GCN layer from Veličković et al. 43 by adding bond features as input information for computing the attention values. Let H t ∈ R n×h denote input node features for the t-th GCN-att layer and N (i) ⊂ Z ≥0 denote set of indices of neighbors of node at index i (where i ∈ N (i)). We calculate new node features H t+1 ∈ R n×h as follows: where as before σ r denotes ReLU activation function and indicates vector concatenation. The scalar K ∈ N + is the number of attention heads and f att : f : R na → R na/K are standard linear layers. Numbers n a , n b , d and K are hyperparameters of the model. We require that n a is divisible by K. We use the same hyperparameter values for all GCN-att layers in the model.

Gradient-based training of MEGAN
In contrast to Do et al. 23 , who use reinforcement learning to train their model, we back-propagate directly through the maximum likelihood objective to train MEGAN. We use teacher forcing to train the model, i.e. to predict each step during reaction generation we use previous steps from the ground-truth as input to the model. This is nontrivial, as computing the gradient of the likelihood objective requires defining a fixed ordering of actions 47 . To solve this issue, You et al. 47 enumerates atoms using breadth-first search.
We adapt a similar idea to reaction generation.
At first, we define a general ordering of action types for forward-and retrosynthesis. This ordering defines which types of actions should be performed first if there are more than one type of actions that could be applied to an atom in a generation step. There are various possible orderings, described in detail in the Section "Ablation study on action ordering". The highest accuracy is achieved in the ordering called BFS RAND-AT. In this ordering, for retrosynthesis, we give bond deletion the highest priority, as it is usually the step that determines the reaction center, which is the first step of reaction prediction in other methods. [36][37][38]48 Analogously, adding a bond has the highest priority for forward synthesis, as this usually determines the reaction center in forward prediction.
Priorities of other types of actions were determined experimentally; we did not observe a significant change in validation error when modifying these priorities. We list these priorities explicitly in Table 1.
The remaining ties between actions are broken as follows. For retrosynthesis prediction, we prioritize actions that act on the least recently modified atoms. Our motivation is to prioritize more difficult actions (such as deciding where the reaction occurs) at the expense of simpler actions (such as completing the synthons resulting from breaking a bond). In forward prediction, we prioritize the most frequently visited atom, which we found to perform better empirically. The remaining ties are broken at random. We describe the algorithm in full detail in the Supplement (Ground truth action ordering).
It is worth noting that the action ordering algorithm is used only to acquire ground truth samples for training MEGAN with a gradient-based method. The model is expected to learn a strong prior towards outputting edits in such order; however, when evaluating MEGAN, we do not explicitly limit the available action space, giving it a chance to determine the order by itself. We observe that MEGAN often predicts the ground truth reaction using a different order of actions than the ground truth training sample.
We use the following rules to resolve conflicts when the action could be outputted at different atoms or bonds: • Stop action can be predicted only by the supernode • Bond actions can be predicted for indices i and j, where i < j and nodes at i and j are atoms Finally, for simplicity, we do not mask out redundant actions, such as deleting a non-existing bond or editing the atom to the same values of properties, as we expect the model to learn not to use such actions.

Experiments
We evaluate MEGAN on three standard datasets for retro and forward synthesis prediction. We first evaluate it on the standard retrosynthesis prediction benchmark USPTO-50k. Next, we investigate how MEGAN scales to the large-scale retrosynthesis task. Finally, we present results on forward synthesis prediction on the USPTO-MIT dataset.

Retrosynthesis
Data First, we evaluate MEGAN performance in retrosynthesis prediction. The goal is to predict correctly the set of reactants based on the product of a reaction. The accuracy is measured by comparing the SMILES 18 representation of the predicted molecules to the SMILES representation of the ground truth target molecules. Before the comparison, we remove atom mapping information from the SMILES strings and canonicalize them using RdKit 45 . We use top-k accuracy computed on reactions from the test set as the main evaluation metric.
We evaluate on the USPTO-50k dataset of approximately 50000 reactions, which was collected by Lowe 49 and classified into 10 reaction types by Schneider et al. 50 . We use the same processed version of the dataset as Dai et al. 13 , where each reaction consists of a single product molecule and a set of one or more reactants, with corresponding atoms between the reactants and the product mapped. Following other studies, we split the reactions into training/validation/test sets with respective sizes of 80%/10%/10%. We use the dataset and split provided by Dai et al. 13 .  51 . We excluded from the comparison some of the concurrent works which had an information leak problem in evaluation (see: Section Concurrent methods excluded from the comparisons in the Supplement).
Results Table 2 reports results on the USPTO-50k benchmark in a variant with and without reaction type information provided. When the reaction type is given, MEGAN outperforms all baselines in all metrics, with the exception of top-1. When the reaction type is unknown, MEGAN beats prior models for K ≥ 3 (and concurrent models for K ≥ 10) and achieves comparable results for lower values of K.
We hypothesize that the advantage of MEGAN for large K stems largely from the fact that MEGAN generates reaction as a sequence of edits. This might help to efficiently search through different plausible reaction centers, hence covering a more diverse subset of the reaction space. It can also enable MEGAN to achieve high coverage of the reaction space, which is indicated by the top-50 accuracy of 93.2% when reaction type is unknown and 95.3% when the reaction type is provided.

Large scale retrosynthesis prediction task
Data The USPTO-50k benchmark, although of relatively high quality, is a small dataset containing only 10 specific types of reactions. We train MEGAN for retrosynthesis on a large benchmark to test its scalability. We use the original data set collected by Lowe 49  Results Table 3 shows top-k accuracy on USPTO-FULL for MEGAN compared to other methods for retrosynthesis prediction. We see that our model achieves competitive performance on a large scale retrosynthesis data set, slightly outperforming other methods in terms of top-10 accuracy. identifies pairwise atom interactions in the reaction center and ranks enumerated feasible bond configurations between these atoms. WLDN5 31 improves this method by combining the problems of reaction center prediction and candidate ranking into a single task. GTPN 23 predicts actions on the graph of substrates, similarly to MEGAN, however, it limits them to actions between existing atom pairs. We include also the concurrent AT model 51 in the comparison.
Experimental Setting We remove the CHIRAL TAG and BOND STEREO features from the model input, as USPTO-MIT has no stereochemical information. Forward synthesis usually takes fewer modifications to predict than retrosynthesis, so we reduce the maximum number of steps to 8 and use a beam size of 20 for evaluation. We also change the strategy to break ties in action ordering. For forward synthesis, we prioritize actions that act on most recently modified atoms. Apart from these changes, we use the same architecture, hyperparameters, and training procedure as for USPTO-50k.
Similarly to Schwaller et al. 15 , we train for two variants of forward synthesis prediction. For the separated variant, compounds that directly contribute to the product are explicitly marked in the set of substrates with an additional atom feature. For the mixed variant, such information is not provided, so the model has a harder task as it has to determine the reaction center from a larger number of possible reactants.
In forward synthesis, MEGAN can generate multiple products for a given set of input substrates set, whereas the USPTO-MIT set contains single-product reactions. We use a simple heuristic to acquire the main reaction product from the MEGAN prediction by selecting the product that has the longest SMILES string. From our observation, this heuristic works correctly for most of the reactions and does not lower the Top K accuracies.
Results training strategy for forward prediction is a promising topic for future work.
In Table 5 we compare the performance of MEGAN and Molecular Transformer on USPTO-MIT in top K predictions. We use the best non-ensemble models provided by Schwaller et al. 15 and set the beam size to 20 during evaluation. We observe that MEGAN surpasses the accuracy of Molecular Transformer for high K values, which again indicates its ability to explore the reaction space efficiently.

Analysis and ablation studies Error analysis
To better understand model performance, we perform here an error analysis focusing on retrosynthesis prediction on the USPTO-50k dataset.
We first visualize in Figure 2 MEGAN prediction on the USPTO-50k test set for random 5 target molecules on which the predicted substrates are different from the ground truth reaction. We observe that the first four reactions are feasible and can be executed using standard methods. The last reaction may pose difficulty because of the nontrivial regioselectivity that depends on the exact reaction conditions utilized.
To better quantify the phenomenon that the top proposed reaction is often correct despite being different from the ground truth, we sampled random 100 reactions from the USPTO-50k test set for which the predicted substrates by MEGAN differed from the ground truth substrates. We then analyzed the correctness, by labeling them into the following categories: (1) no issues detected (seen as a true chemical reaction by a chemist), (2) incorrect chirality of the substrates, (3) the reaction has a low yield and/or important side products, (4) the reaction can be executed but only in multiple stages (e.g. with use of protecting groups), (5) there is a reactive functional group that participates in the reaction instead, (6) incorrect reaction for another reason. All 200 reactions were shuffled and assigned randomly to two experienced organic chemists.
We show the break-down of MEGAN errors by these categories in Figure 3 and examples of reaction in each category in Figure 3. On the whole, in 79.6% cases, the top prediction by MEGAN was deemed correct by chemists even though it differed from the ground truth prediction.
In comparison, 89.5% of ground truth reactions were considered correct. Incorrect chirality of the substrates was the most common reason for reaction to be labeled as incorrect (7.1%, compared to 0.5% of ground truth reactions).
The main outcome of this error analysis is that for retrosynthesis prediction on the USPTO-50k, MEGAN performance seems to be close to human-level performance in the sense that the proposed reactions start to be hard to distinguish from ground-truth reactions in terms of their correctness.
Having said that, there is still significant room for improvement, for example in terms of the chirality of the proposed substrates. Another common source of errors was ignoring the existence of a reactive functional group in the substrates.

Ablation study on action ordering
Our gradient-based method of training MEGAN is based on a fixed ordering on ground-truth actions. Defining this ordering requires breaking ties between different actions of the same type.
For example, if there are several EditBond actions, we have to decide which should be executed first by the model.
In retrosynthesis prediction tasks, we broke ties based on the last time a given atom was modified.
Actions that modify atoms that were least frequently modified are prioritized. The remaining ties were broken at random. We refer here to this ordering algorithm as BFS RAND-AT.
We explore here two modifications to this ordering. First, instead of prioritizing actions that Figure 2: The top-ranked retrosynthesis prediction by MEGAN for 5 random target molecules from the USPTO-50k test set on which the prediction is different from the ground truth reaction. The first four reactions are feasible and can be executed using standard methods. The last reaction may pose difficulty because of the nontrivial regioselectivity that depends on the exact reaction conditions utilized.   Figure 3: Categorization of the top-ranked retrosynthesis predictions by MEGAN on a subset of USPTO-50k test set on which the top prediction differs from the ground truth reaction. For reference, we also include labels assigned to ground truth reaction from the test set. We find that the top-ranked prediction is often (in 79.6% cases) deemed correct by a chemist even though it differs from the ground truth ("No issues detected"). The most common source of error is incorrectly predicting chirality in the substrates ("Incorrect chirality", 7.1%).  Figure 3).

Target product
Ground truth reactants

Predicted reactants
Only possible in multiple steps Low yield or side products Incorrect chirality A reactive functional group ignored modify the least frequently modified atom, we experiment with prioritizing actions that modify the most frequently modified atom. We refer to this as DFS ordering. Second, we investigate a strategy to break remaining ties based on the canonical SMILES. In this strategy, we prioritize actions that modify atoms that are earlier in the ordering determined by the canonical SMILES. We refer to this as CANO-AT. By combining these different choices, we arrive at four different ordering strategies: BFS RAND-AT, DFS RAND-AT, BFS CANO-AT, and DFS CANO-AT.
We present the performance of all variants on the USPTO-50k dataset in Table 6. We observe that BFS RAND-AT ordering achieves the highest performance across most values of K, which motivates our choice to use it for training MEGAN on retrosynthesis prediction tasks. It is worth emphasizing that the choice of action ordering can strongly impact the performance of MEGAN.
For example, the RANDOM action ordering achieves over 4% lower top-1 accuracy compared to DFS RAND-AT, which we used in the experiments. Performance with respect to reaction popularity Next, we analyze how MEGAN performance depends on reaction popularity. To approximate reaction popularity, we count the number of times that the corresponding template occurred in the training set. We use code provided by Coley et al. 31 to extract templates for each ground-truth reaction from the datasets. Crucially, MEGAN generates reaction as a sequence of graph edits, which is inspired by how chemists represent chemical reactions. This inductive bias might explain the strong empirical performance of the model. In particular, we argued that it enables the model to more efficiently search through the space of plausible reactions. Looking forward, generating reactions as a sequence of edits is promising for building more intuitive human-computer interaction in synthesis planning software.
A natural topic for the future is reducing the reliance of MEGAN on the mapping between products and substrates. We also found that incorrectly predicting chirality of the output was the most common source of errors on the USPTO-50k dataset. Improving these aspects has the potential to further push state-of-the-art in retro and forward synthesis prediction.

Data and Software Availability
We open-sourced our code and trained models at https://github.com/molecule-one/megan. Our implementation is based on freely available software and is released under the MIT license. All datasets are available online, with references given in the main text.

Supplementary Material Featurization
We featurize atoms and bonds with one-hot encoded vectors of features calculated using RdKit 45 .
On each dataset, we select features that allow for correct reconstruction of all products and substrates from the development (training+validation) set. In Tables   For both atoms and bonds, we add a special IS EDITED feature that marks all bonds and atoms that have been modified by actions. This aims to help the decoder to focus on these atoms and bonds, as they are the most probable candidates for the next generation steps.

Graph edit actions
In Table 8 we show all possible graph actions found on the USPTO-50k development set. The actions were found during the generation of the training samples, which we describe in the next section. The actions have different sets of parameters, depending on the action type. BondEdit actions are bond actions, that is they are predicted for a pair of atoms. All other types of actions are atom actions and are predicted for a single atom in the graph.

Ground truth action ordering
An important design choice in training of MEGAN is the algorithm that determines ordering of ground-truth actions that transform the product into the substrates (or vice versa in forward prediction). This ordering allows us to backpropagate through the maximum likelihood objective on the training set.
Action type priorities At first, we define a general ordering of action types for forward-and retrosynthesis in Table 1. This ordering defines which types of actions should be performed first if there are more than one type of actions that could be applied to an atom in a generation step. For retrosynthesis, we give bond deletion the highest priority, as it is usually the step that determines the reaction center, which is the first step of reaction prediction in other methods. [36][37][38]48 Analogously, adding a bond has the highest priority for forward synthesis, as this usually determines the reaction center in forward prediction. Priorities of other types of actions were determined experimentally; we did not observe a significant change in validation error when modifying these priorities.
Atom priorities Given the ordering of action types described in Table 1, we still need to decide on which action to choose in a scenario when two of the ground-truth graph edits are possible for different atoms. For example, we might need to decide on the order of adding leaving groups to two atoms that were disconnected by deleting a bond between them. Ordering of such actions can be determined by two factors: 1. Whether an atom has already been modified (BFS vs DFS) 2. Order of the atom in the canonical SMILES of the target compounds (CANO-AT vs RAND-AT) In the case of EditBond actions that act on pairs of atoms, their priority is determined based on the highest of the priorities of these atoms. We present the whole algorithm for determining the next action in the Algorithm 1. Let Act(a) = all actions that lead to T that could be performed on a Let Act (a) = all actions from Act(a) with the action type that has the highest priority Let ActionT ype = action type of actions from Act (a) if ActionT ype = EditAtom then By definition, there could only be a single EditAtom action performed on an atom. return The single EditAtom action from Act (a) else if ActionT ype = EditBond then Let b = the second atom from actions in Act (a) ranked highest according to the atom ordering return The EditBond action between a and b else if ActionT ype = AddAtom then Let b = the second atom from actions in Act (a) ranked highest according to the atom ordering return The AddAtom action between a and b else if ActionT ype = AddBenzene then return A random AddBenzene action from Act (a) end if Hyperparameter search Table 9 shows the hyperparameter values used for the final model on the USPTO-50k. We found these hyperparameters manually based on the performance on the USPTO-50k validation set. For other datasets, we only slightly modified these hyperparameters (as described in the main text).

Concurrent methods excluded from the comparisons
After a detailed examination, we decided to exclude from the benchmark comparisons some of the concurrent methods of retrosynthesis prediction that suffer from an information leak problem between the training and test sets. This problem is described in details at https://github.com/utasmile/RetroXpert (accessed on May 7th, 2021). The methods with the information leak problem include GraphRetro 37 and RetroXpert. 38