From medicines to materials, small organic molecules are indispensable for human well-being. To plan their syntheses, chemists employ a problem solving technique called retrosynthesis. In retrosynthesis, target molecules are recursively transformed into increasingly simpler precursor compounds until a set of readily available starting materials is obtained. Computer-aided retrosynthesis would be a highly valuable tool, however, past approaches were slow and provided results of unsatisfactory quality. Here, we employ Monte Carlo Tree Search (MCTS) to efficiently discover retrosynthetic routes. MCTS was combined with an expansion policy network that guides the search, and an "in-scope" filter network to pre-select the most promising retrosynthetic steps. These deep neural networks were trained on 12 million reactions, which represents essentially all reactions ever published in organic chemistry. Our system solves almost twice as many molecules and is 30 times faster in comparison to the traditional search method based on extracted rules and hand-coded heuristics. Finally after a 60 year history of computer-aided synthesis planning, chemists can no longer distinguish between routes generated by a computer system and real routes taken from the scientific literature. We anticipate that our method will accelerate drug and materials discovery by assisting chemists to plan better syntheses faster, and by enabling fully automated robot synthesis. Retrosynthetic analysis is the canonical technique to plan the synthesis of organic small molecules, for example drugs, agro-and fine chemicals, and is part of every chemist's curriculum. [1, 2] In retrosynthesis, a search tree is built by "working-backwards", analysing molecules recursively and transforming them into simpler precursors until one obtains a set of known or commercially available building block molecules. [3] The transformations are reversed chemical reactions: A chemist can take the plan, and execute it in the lab in the forward direction to synthesize the target compound. Retrosynthetic analysis is a formidable intellectual exercise that demands broad and deep chemical knowledge. Its central element requires a combination of creativity and a pattern recognition process, in which functional groups (patterns of atoms and bonds) of the target molecule are matched to transformations. [4] Transformations are derived from successfully conducted series of similar reactions with analogous starting materials, and often named after their discoverers. [5] At each retrosynthetic step, out of thousands of known transformations in modern chemistry, chemists intuitively rank the most promising transformations highly, and do not even actively think about the unreasonable ones. 1 Unfortunately, when a transform is applied to a new molecule, there is no guarantee that it will generalize and the corresponding reaction will actually proceed in the expected way. [7] A molecule failing to react in the way predicted by the transform is called 'out of scope'. This can be due to steric or electronic effects, an incomplete understanding of the reaction mechanism, or conflicting reactivity in 1 In cognitive science, this can is explained via "System 1" and "System 2" in dual process theory. [6] the molecular context, which may hamper a reaction and lead to chemo-, regio-or stereoselectivity issues. Predicting which molecules are 'within scope' can be challenging even for the best human chemists, who perform synthesis planning using knowledge and mechanistic reasoning acquired during long years of study and the many intuitive models the chemical community developed over the last 250 years, in combination with extensive literature research for every transform step, which may take hours. [7, 3] Computer-assisted synthesis planning (CASP) could help chemists to find better routes faster, and is a missing component in virtual de novo design and robot systems performing molecular designsynthesis-test cycles. [8] [9] [10] [11] To perform CASP, the knowledge that humans "simply learn" has to be transferred into an executable program. [12, 4, [13] [14] [15] In retrosynthetic analysis, the sequential transformation of the target molecule gives rise to a structure called search tree. Figure 1 translates the traditional chemist's representation to a search tree representation. [13] The nodes in the tree represent the synthetic position. Positions are sets of molecules; they contain all precursors needed to make the molecules of the preceding positions all the way down to the tree's root node (which contains the target molecule). Branches in the search tree then correspond to complete routes, which allows to compare and rank different routes. Search in retrosynthesis poses a significant challenge: The branching factor (number of applicable transformations per step) can be quite large and range from 80 for modest systems [13] up to more than 40000 (see below). Even though the search depth rarely exceeds more than 15-20 steps, the search tree explodes combinatorially if solutions are explored exhaustively. In the past, heuristic best first search (BFS) has been employed for synthesis planning. [13] Heuristic functions determine position values, which can guide the search into promising directions. Unfortunately, unlike in chess, it is difficult to define strong heuristics in chemistry for two reasons: First, chemists tend to disagree about good positions, that is how easy it is to make a molecule. [26, 27] Second, the value of a position depends highly on the availability of suitable precursors. [15, 13] Even complex molecules can be easily made in a few steps if there are readily available building blocks. Therefore, one cannot reliably estimate the value of a synthetic position without completely "playing the molecules until the end of the game". Recently, Monte Carlo Tree Search (MCTS) has emerged as a highly general search technique, which is particularly suited for sequential decision problems with large branching factors without strong heuristics, such as games or automated theorem proving. [28] [29] [30] These properties suggest that MCTS could be effective in synthesis planning as well. MCTS builds a search tree by iterating four distinct phases, see Figure 2a . During selection (1), the most important position to analyse next is selected according to its expected value and the confidence in that value. Expansion (2) performs 2 one step of retroanalysis on the chosen position by adding possible precursor positions. One of these precursor positions is then selected for evaluation by a rollout (3). A rollout is a Monte Carlo simulation, in which random search steps are performed without branching until a solution has been found by transforming the position into building blocks or a maximum depth is reached. During update (4), the position values in the tree are updated to reflect the outcome of the rollout. As the tree grows and many iterations are run, the position values become more accurate, and eventually converge to the optimal solution. For the game of Go, Silver et al. elegantly demonstrated with the AlphaGo system that MCTS can be augmented with neural network policies p(t|s), which predict the probability of taking move (applying transform) t in position s, and are trained to predict the winning move based on grandmaster games. Such policies were used to inform selection and rollout instead of taking random moves. [31] [32] [33] In this work, inspired by the recent progress in Go culminating in AlphaGo, [31] we combine three different neural networks (3N) together with MCTS to perform chemical synthesis planning (3N-MCTS). The first neural network (expansion policy), trained on how molecules were made in the past, transfers its knowledge to novel molecules. It guides the search into promising directions by proposing a restricted number of possible, automatically extracted transformations. A second neural network then predicts if the proposed reactions are actually feasible/in scope. Finally, to estimate the position value, transformations are sampled from a third neural network (policy) during rollout. The neural networks were trained on essentially all reactions published in the history of organic chemistry. We extracted transformation rules from 12.4 million single step reactions from the Reaxys[34] database using our previously reported protocol, which took approx. 3 hours on a laptop. [24] The database contains the bulk of organic and organometallic reactions published since the inception of the field in 1771 until today. [34] Two sets of rules were extracted: The rollout set comprises rules which contain the atoms and bonds changed in the course of the reaction (the reaction center), and the first degree neighbouring atoms. Here, all rules which occurred at least 50 times in reactions published before 2015 were kept. For the expansion rules, a more general rule definition was employed. Here, only the reaction center is extracted. Rules occurring at least 3 times were kept. The two sets encompass 17,134 and 301,671 rules, and cover 52% and 79% of all post-2014 chemical reactions, respectively. 2 Rule extraction associates each reaction, and thus each product with a transformation rule. As previously reported, this allows to train neural networks as policies to predict the best transformations given the product, or in other words, the best reactions to make the product. [24] Importantly, such neural networks also learn about the context in which the reactions can occur (functional group tolerance). [24] As the expansion policy, we employed a deep highway network with exponential linear unit (ELU) nonlinearities. [35, 36] To assess its generalisation ability, we performed a time-split strategy. For training, only reactions published before 2015 were used, while for validation and testing, data from 2015 and later was selected. Time-split validation has been shown to give a more realistic assessment of predictive performance than other splitting strategies, such as random splitting. [37] Table 1 shows the metrics for the expansion policy. The neural network predicts the correct solution out of 301,671 transformations with an accuracy of 31%, which is reasonable. It has to be noted that there are almost always many feasible ways to make a molecule. Top10 and top50 accuracies of 63.3% and 72.5% indicate that the correct transformations are generally ranked highly. Beyond the top50 predicted results, the accuracy only increases marginally. This observation allows to drastically reduce the branching factor. During search tree expansion, we restrict the possible 2 It has to be noted that even when taking all pre-2015 rules without count restriction into account, only 82% of the reactions published after 2014 are covered. The missing 18% are novel reaction types. This highlights the success of chemists inventing novel methodologies, but also implies that eventually a retrosynthesis system should be also able to discover novel reactions on its own. In phase 2) this node may be expanded by processing the molecules of the position A with the Expansion procedure (b), which leads to new positions B,C which are added to the tree. Then, the most promising new position is chosen, and a rollout phase 3) is performed by randomly sampling transforms from the rollout policy until all molecules are solved or a certain depth is exceeded. In the update phase 4), the position values are updated in the current branch to reflect the result of the rollout. b) Expansion procedure: In the step, the molecule m t to retroanalyse is fed into the policy network to predict the best transforms. Then, only the top k transforms are applied to the molecule, which yields complete reactions (Rxn). For each of these reactions, reaction prediction is performed using the transition classifier. Improbable reactions are then filtered out, which leads to the list of admissible actions and corresponding precursor positions B,C. (TF = Transformation, Rxn=Reaction, ECFP = Extended Connectivity Fingerprint) transformations to a maximum of 50. Additionally, we sum up the probabilities of the predicted actions starting from the highest ranked transformation. When the cumulative probability reaches 0.995, we stop further expansion, even if less than 50 actions have been expanded. This allows the system to focus on highly likely transformations in cases where only few good options exists to make a molecule, for example acyl chlorides or Grignard reagents. 3 We observed that the reactions in this reduced top50 are almost exclusively reasonable, and often variations of the correct prediction. For example, a Heck reaction can usually be conducted with bromide, iodide or triflate as the leaving group. The rollout policy network, which is a neural network with one hidden layer, is trained in the same way as the expansion policy. It uses a set of 17,134 rules, which leads to a lower coverage, yet it needs just 10 ms to make a prediction, in contrast to 90 ms for the expansion policy (see Table 1 ). 3 Additionally, this entails that the NP-complete subgraph isomorphism problem, which determines if the corresponding rule can be applied to a molecule and yields the next molecule(s), needs only to be solved for a subset of the rules. The rationale for using two different rulesets is to have a powerful, but slow policy to select the best candidate transformations for expansion, and a fast roll policy to estimate the position values. [31]  After the search space has been narrowed down by the expansion policy to the most promising transformations, we need to predict whether the corresponding reactions will actually work for a particular molecule. We trained a deep neural network as a binary classifier to predict whether the reactions corresponding to the transformations selected by the policy network are actually feasible. [38, 39] The classifier has to be trained on successful and failed reactions. Unfortunately, failed reactions are only seldom reported and not contained in reaction databases. Therefore, starting from reported reactions, negative reactions, for example with incorrect regio-and chemoselectivity were generated using an established protocol [40] [41] [42] together with negative sampling. Using this data augmentation strategy, we generated 100 million negative reactions from reactions prior to 2015 for training and 10 million after 2014 for testing. As positive cases, reported reactions from these periods are used. The classifier only receives the product and the reaction fingerprint as its input -reaction conditions or additional reagents or catalysts are not considered. The classifier thus only serves as rapid filter to predict if a reaction will generally work or not. On the test set, the classifier achieves an area under the Receiver Operation Characteristic (ROC) curve of 0.99, and the area under the precision-recall curve of 0.94, which indicates good performance. [43] The expansion policy network and the in-scope filter network are then combined to a pipeline (Figure  2 b) ). When a position A is to be analyzed, each molecule of the position is fed into the policy network. Then, the transformations with the highest scores are applied to the molecule, which yields the possible precursors and thus full reactions. These reactions are submitted to the in-scope filter, where only transformations and precursors corresponding to positively classified reactions are kept. They represent the "legal moves" available in position A. The expansion procedure and the rollout policy are then incorporated in the respective phases of a Monte Carlo Tree Search algorithm to form our 3N-MCTS ansatz. The four phases are iterated to build the search tree. In the first 3N-MCTS phase, starting at the root node (the target molecule) of the search tree, the algorithm sequentially selects the most promising next position within the tree until a leaf node is reached (Fig. 2a) . The algorithm balances the selection of high-value positions and unexplored positions. If a leaf nodes is visited for the first time, it is directly evaluated by a rollout. If it is visited for the second time, it is expanded by processing via the expansion policy. 2. Expansion In the expansion, the possible transforms determining the follow-up positions of the current position are selected by applying the expansion procedure. The predicted follow-up positions are added to the tree as children of the leaf node, and the most promising position is selected for rollout. 3. Rollout In the rollout phase, first, the status of the position is checked. If the position is already solved, the algorithm directly receives a reward greater than 1. [44] Non-terminal states are subjected to a rollout, where actions are sampled from the rollout network recursively, until the state has been deconstructed into building blocks or a maximal depth is reached. 4. Update If a solution has been found during rollout, a reward of 1 is received. Partial rewards are given if some but not all molecules in the state are solved. If no solution was found, a reward of -1 is received. In the last phase, the tree is updated to incorporate the achieved reward by updating the position values. These four phases of 3N-MCTS are iterated until a time budget or maximal iteration count is exceeded. Finally, to obtain the synthesis plan, we repeatedly select the retrosynthetic step with the highest value until a solved position is reached, or a maximum depth has been exceeded, in which case the problem is unsolved. To evaluate the performance of 3N-MCTS, we compare our algorithm to the state of the art search method, which is Best-First Search (BFS) with the hand-coded SMILES 3 2 heuristic cost function, as described by Grzybowski and coworkers (heuristic BFS). [13] This function assigns the lowest cost to reactions that split up the molecule in equally sized parts. Additionally, we perform BFS with the cost calculated by the policy network (neural BFS). All algorithms use the same set of automatically extracted transformations. The evaluation is performed in a time-split way: The models were trained only on data published before 2015. As test data, only molecules first reported after 2014 are considered, which were not contained in the training data. The different algorithms were provided with a set of target molecules, the algorithms had to find a synthesis route to known building blocks. The building blocks have to be selected before the search is started, and could be molecules on stock in the lab, known in the literature, or commercially available chemicals. Herein, we use a set of 423,731 molecules, containing 84,253 building blocks from three major chemical suppliers (SigmaAldrich, AlfaAesar and Acros), 5 and 339,478 molecules from the Reaxys database, which have been used as reactants at least 5 times before 2015. A researcher may choose different sets of building blocks for each search, e.g. first try to find the solution to a problem with molecules on stock in the lab, and afterwards consider additional molecules from chemical suppliers. Quantitative evaluation Surprisingly, in the past, neither hand-coded nor automatically extracted retrosynthetic systems have ever been validated in a statistical way. To obtain a set of target molecules that contain different scaffolds, first all molecules reported after 2014 were clustered using the extended connectivity fingerprint (ECFP6)-based Butina algorithm [46] . Then, 497 target molecules were randomly selected from amongst the different cluster cores. This test set was employed to quantitatively assess the performance characteristics of the different search algorithms by finding synthesis routes to known building blocks (see Figure 4) Figure 4: Influence of the time per query (max. steps = 100000) on performance. MCTS already solves more than 80% of the test set with a time limit of 5 seconds per target molecule, compared to 40% with neural BFS and 0% for heuristic BFS. MCTS solved 92% of the test set with a limit of 60 s per molecule, whereas neural BFS solved 71%, and heuristic BFS solved 4%. Even at much longer runtimes of 20 minutes per molecule, heuristic and neural BFS are not able to compete with MCTS. Provided with infinite runtime, however, the algorithms will converge to the same performance. To determine which components of MCTS are responsible for its superior performance, we compared MCTS against several related search algorithms (see Table 2 ) at a runtime limit of 3×300 s (three restarts). MCTS in conjunction with the expansion policy network solved the highest number of retrosynthetic targets. On average, MCTS required the least amount of time per molecule to find a solution (entry 1). Plain Monte Carlo search (MC) randomly selects transforms using the expansion policy network, without building a tree. MC (entry 2) solved 89.54% UCT is an MCTS variant which uses the expansion policy network only to narrow down the possible transforms, but not to guide the search via the predicted probability of the transform. This way, 87.12% of the test set is solved. BFS using a cost function based on the expansion policy network only solved 84.24%, highlighting the importance of rollouts. The traditional approach, BFS with a hand-designed heuristic cost function, solves only 45.6% of the test set, and needs 433.4 s on average to find a solution. These results suggest that all components of 3N-MCTS (building a tree, reducing the branching factor via the expansion policy, guiding the search with the expansion policy, and using rollouts) contribute to its superior performance. The central criticism of retrosynthesis systems has been that the proposed routes often contain what chemists immediately recognize as chemically unreasonable steps. Therefore, to assess the quality of the solutions we conducted two AB tests, in which 45 graduate level organic chemists from two 7 world-leading organic chemistry institutes in China and Germany had to choose one of two routes leading to the same molecule based on personal preference and synthetic plausibility. The tests were double blind, meaning that neither participants nor conductors were aware of the origin of the routes. The test molecules were selected randomly from a set of drug-like compounds first published after 2014. Figure 5: a) Chemists did not prefer literature routes over routes found by MCTS b) Chemists significantly prefer routes by 3N-MCTS over routes by heuristic BFS without policy network and in-scope filter. As ratio above 0.5 indicates that more than 50% of participants preferred the MCTS solution. In the first test, the participants had the choice between a route reported by expert chemists in the literature, and a route generated by our 3N-MCTS algorithm for the same target molecule. Routes to nine different target molecules were offered. Routes towards the same molecule were required to have the same number of steps. Here, one might expect that the participants can clearly identify the routes suggested by the machine as inferior. Surprisingly, this is not the case. We found that the experts did not significantly prefer the literature route (43.0%) over our programs' route (57.0%). Figure 5 a) shows the preferences ratios on the individual routes. Here, the preference is generally balanced, with a slight trend towards MCTS. In some cases, the participants have clear preferences (see Fig. 6a and Extended Methods for a examples where the chemists did not prefer MCTS). In the second test, the participants had to report their preferences for either routes found by 3N-MCTS or routes generated by a baseline system, which uses heuristic BFS, the same transformation rules as 3N-MCTS. However, it lacks the policy network to preselect promising transformations and the in-scope filter to exclude unlikely steps. Here, the participants significantly preferred the routes generated by the MCTS algorithm (68.2%) over the baseline system (31.8%). We attribute the preference towards the 3N-MCTS generated routes to lower frequencies of unreasonable steps (see Figure 6 b)). 8 In this work, we have shown that Monte Carlo Tree Search combined with deep neural network policies and an in-scope filter can be used to effectively perform chemical synthesis planning. One advantage of our approach is that it can provide full retrosynthetic pathways within seconds. In contrast to earlier approaches, our purely data-driven ansatz can be initially set up within a few days without the need for tedious and biased expert encoding or curation, and is equally applicable to both in-house and discipline-scale datasets. We have demonstrated that our 3N-MCTS approach has the best performance characteristics in terms of speed and number of solved problems compared with established search methods. In the past, retrosynthetic systems have been criticized for producing more noise than signal. Here, we have shown by double blind AB testing that organic chemists did not generally prefer literature routes over routes found by 3N-MCTS. We observed that heuristic best-first search without neural network guiding did lead to many unreasonable steps being proposed in the routes, while the 3N-MCTS approach proposed more reasonable routes. This is supported by the double-blind AB experiments, where the participating organic chemists showed clear preference towards 3N-MCTS over the traditional approach.
Impressively, the results attained by Segler et al. are shown to perform at a level comparable to human-driven organic molecule synthesis planning REF .
Learning to Plan Chemical Syntheses
Nature 555 (2018), 604-610
Computer Science Physics