A user-tunable machine learning framework for step-wise synthesis planning

Abstract

We introduce MHNpath, a machine learning-driven retrosynthetic tool designed for computer-aided synthesis planning. Leveraging modern Hopfield networks and novel comparative metrics, MHNpath efficiently prioritizes reaction templates, improving the scalability and accuracy of retrosynthetic predictions. The tool incorporates a tunable scoring system that allows users to prioritize pathways based on cost, reaction temperature, and toxicity, thereby facilitating the design of greener and cost-effective reaction routes. We demonstrate its effectiveness through case studies involving complex molecules from ChemByDesign, showcasing its ability to predict novel synthetic and enzymatic pathways. Furthermore, we benchmark MHNpath against existing frameworks using the PaRoutes dataset, achieving a solution rate of 85.4% and replicating 69.2% of experimentally validated “gold-standard” pathways. Our case studies show that the tool can generate shorter, cheaper, moderate-temperature routes employing green solvents, as exemplified by molecules such as dronabinol, arformoterol, and lupinine.

---

1 Introduction

The integration of machine learning (ML) into chemical synthesis has transformed the field of computer-aided synthesis planning (CASP), providing chemists with powerful tools to design and execute synthetic routes more efficiently. Traditional retrosynthetic analysis, which relies heavily on expert intuition and experience,^1,2 is increasingly being augmented by data-driven approaches³ to break down complex molecules into simpler precursors using ML techniques to predict viable synthetic pathways.^1,4 These advancements are particularly crucial in addressing the growing complexity of modern synthetic challenges, where exploring vast chemical spaces catalyzes the discovery of novel molecules. Examples of such challenges include the synthesis of complex natural products with multiple stereocenters⁵ and the need for highly selective functional group transformations in pharmaceutical development.⁶

In CASP, clean chemical data is the foundation for computationally designing synthetic routes. These data enables training of a predictive model that finds precursors while an efficient search algorithm navigates the expansive chemical space to propose feasible retrosynthetic pathways. State-of-the-art ML models for CASP can be categorized into template-based and template-free approaches where templates are generalized representations that encapsulate the core chemical transformation patterns inherent to each reaction. Complementing the ML model and search algorithm is an integrated scoring system that evaluates routes based on multiple criteria, such as cost and reaction temperature. Together, these three components establish a cohesive framework that underpins modern computer-aided synthesis planning.

Recent advancements in ML, including the development of deep neural networks and transformer-based models, have significantly enhanced the predictive capabilities of template-free methods in CASP.^7–11 Transformer models have emerged as a powerful tool in retrosynthetic prediction due to their ability to handle sequence-to-sequence tasks effectively.¹² The molecular transformer model,⁸ and its variants,^9–11 for instance, has demonstrated remarkable accuracy in predicting both reactants and reaction conditions by treating retrosynthesis as a language translation problem from products to reactants. These models leverage self-attention mechanisms to capture complex dependencies within chemical reactions, making them highly effective for single-step retrosynthesis.^13,14 Some graph-based methods^15,16 predict the products from reactants in an auto-regressive fashion, while others^17–19 predict the products by anticipating the final states of bonds or electrons. However, the reliance of template-free approaches on extensive labeled reaction data, computational resources, and single-step focus has been a significant limitation in CASP.

On the other hand, template-based approaches have also been fundamental in CASP.²⁰ These approaches directly apply predefined reaction rules to fragment target molecules into simpler structures, making them more interpretable and chemically intuitive compared to template-free approaches. Segler et al.²¹ and Coley et al.²² were among the first to use automatically extracted rules to predict outcomes of organic reactions using neural networks while Zhang et al.²³ extend it by also continually evolving them and Rho et al.²⁴ give a global outlook by effectively abstracting the detailed substructures. Chen et al.²⁵ proposed a chemistry-motivated graph neural network that uses templates to describe the net changes in electron configuration between the reactants and products. Recent advancements, such as the use of modern Hopfield networks for template prediction, have improved the performance of template-based CASP tools.²⁶ Modern Hopfield networks allow for the efficient prioritization of reaction templates by leveraging associative memory mechanisms, significantly improving the speed and reliability of retrosynthetic planning. Nevertheless, template-based methods often struggle with the rigidity of predefined rules, analyzed by Choe et al.,²⁷ which can limit their ability to propose novel synthetic routes or adapt to emerging reaction types.

Computational retrosynthesis has traditionally been treated as a tree search problem, where each step involves searching for chemically feasible precursors to derive the product molecule. Early works have used greedy search,⁴ but more recent efforts have adopted Monte Carlo tree search²⁸ and A*-like algorithms.^29,30 These search methodologies attempt to optimize the exploration of synthetic routes by balancing the exploitation of high-confidence reaction pathways with the exploration of novel or underutilized transformations. However, their effectiveness is contingent on the quality and comprehensiveness of the reaction databases they utilize.

In addition to the different components of CASP, the integration of biocatalysis in computational retrosynthetic planning has gained attention for its potential to enhance sustainability. Biocatalysis leverages enzymatic transformations to achieve highly selective reactions under mild conditions, reducing the environmental footprint of chemical synthesis.³¹ Incorporating enzymatic steps into synthetic pathways is particularly advantageous for pharmaceutical and fine chemical industries, where stringent purity and selectivity requirements are critical.^32,33 Tools like RetroBioCat³¹ have demonstrated the ability to design biocatalytic cascades that complement synthetic routes by introducing enzymatic steps that are often more selective and environmentally friendly. However, even these specialized tools are limited by the availability of comprehensive enzymatic reaction data and their focus on narrow reaction types.

In this work, we introduce MHNpath, an ML-driven retrosynthetic tool designed to help chemists explore greener, cost-effective, and moderate-temperature synthetic pathways. First, we develop a robust template-prioritization model based on modern Hopfield networks,³⁴ enhanced with Xavier initialization, dropout, and L2 regularization to improve stability and generalization. We additionally introduce two new evaluation metrics that enable a more comprehensive comparison of template prioritizers. Second, we implement a user-tunable scoring system that integrates practical constraints and include precursor cost, reaction temperature, and solvent toxicity, allowing chemists to steer the search toward more sustainable and operationally feasible routes. Third, we demonstrate MHNpath's effectiveness through case studies from PaRoutes³⁵ and ChemByDesign,³⁶ showcasing its ability to identify viable, interpretable multi-step pathways. Together with its curated enzymatic and synthetic template datasets and a global greedy tree-search algorithm, MHNpath provides a flexible and data-efficient framework that addresses the limitations of existing CASP tools, particularly their reliance on extensive training and rigid template-based systems. The following sections will detail the methodology behind MHNpath, including the data processing pipeline (Section 2.1), the implementation of the modern Hopfield network-based template prioritizer (Section 2.2), and the tree search methodology and custom scoring employed for retrosynthetic exploration (Sections 2.3 & 2.4), as outlined in Fig. 1. Furthermore, we will present results (Section 3) comparing MHNpath against existing CASP tools and provide case studies demonstrating its practical applications.

Fig. 1 (a) Data processing and model architecture. The model development follows three primary stages. Step 1 (data curation): raw enzymatic (BKMS-react, Rhea) and synthetic (USPTO) databases are cleaned and filtered to extract unique reaction templates. Step 2 (encoding): the target molecule m and the extracted reaction templates t are converted into ECFP and RDKit fingerprints, respectively, and passed through dedicated neural network encoders employing Xavier initialization and batch normalization to produce dense, continuous representations m_h and T_h. Step 3 (template prioritization): a modern Hopfield network projects m_h and T_h into a shared associative memory space, where an attention-like update mechanism produces the new molecule state ξ^new by retrieving the most applicable reaction templates, with the model optimized via cross-entropy loss. (b) Tree search methodology. A weighted global greedy tree search explores the retrosynthetic space iteratively in Step 4 (search and scoring). At each iteration, the algorithm (i) selects the highest-scoring node from a priority queue, (ii) expands it by applying the prioritized templates to generate new precursor nodes (grey circles for enzymatic, yellow for synthetic), and (iii) updates the priority queue using a multi-objective scoring function P′ = f(C′, T′, S′) that promotes low precursor cost, green solvent usage, and moderate reaction temperatures. This cycle continues until buyable precursors (light blue squares, <$100 per g) are found, the maximum search depth is reached, or the allotted time limit is exhausted.

2 Methods

The MHNpath framework is presented in Fig. 1. A modified modern Hopfield network³⁴ based architecture by Seidl et al.²⁶ is adapted. Each prioritizer (also referred to as model) accepts a target molecule as input and predicts a ranked list of the most applicable reaction templates, which are applied to the target molecule to get precursors (nodes). The individual scores for precursor cost, reaction temperature, and solvent toxicity are calculated and summed using user-tunable weights. The precursors are further explored in a recursive manner using a global greedy tree search algorithm where the best-scoring unexplored node from the entire tree is explored first. Fig. 1(b) gives an example of the multi-modal-guided global greedy tree search approach.

2.1 Dataset processing

Our study leveraged two comprehensive training datasets to develop the predictive component of the MHNpath framework. Additionally, specific external datasets were curated for benchmarking and evaluation tasks. In this context, we selected training data that provided a comprehensive representation of both enzymatic and synthetic reactions. The training dataset of enzymatic reactions was compiled from the BKMS³⁷ dataset and the Rhea³⁸ database, initially containing about 68 000 reactions. Following a rigorous preprocessing procedure that involved removing duplicates, invalid SMILES representations, and unbalanced reactions, we obtained 35 289 unique and valid reaction SMILES. From this refined enzymatic subset, we extracted 17 047 reaction templates using the RDKit³⁹ library. This subset was used to train a single enzymatic template prioritizer, utilizing an 80 : 10 : 10 random split for training, validation, and testing, respectively.

We also used a dataset of synthetic organic reactions for training which was sourced from the USPTO⁴⁰ dataset. Starting with 1 808 937 reactions and applying a series of preprocessing steps resulted in 1 693 109 refined reactions, from which 301 242 reaction templates were extracted in the synthetic subset. To manage computational resources and memory overheads due to the large number of data points in the synthetic subset, we employed a scaffold split technique using Datamol⁴¹ to ensure high chemical diversity while dividing the subset further into five equally sized subsets. This partitioning strategy allowed us to train five distinct synthetic template prioritizers (one per divided synthetic subset). Each divided subset underwent an 80 : 10 : 10 random split for model development. During inference, the predictions from these five models are aggregated using an ensemble approach. This strategy significantly enhances robustness with a single synthetic template prioritizer achieving a top-1 accuracy of 36.6% on its own test set and the ensemble method reaching a top-1 accuracy of 42.2% on the combined test set. In this context, top-1 accuracy denotes the model's ability to accurately prioritize the correct ground-truth template as the number one rank in the prediction list.

Once the predictive components of MHNpath were developed, we rigorously assessed the framework's performance in Sections 3.2 and 3.3 using three tasks as described below. These tasks required curating distinct datasets from the literature and provide unbiased performance metrics for systems not present in the training set.

• PaRoutes benchmark: to evaluate pathway reconstruction against patent literature, we utilized the PaRoutes³⁵ dataset. We applied a scaffold split to identify the ten most commonly occurring scaffolds within PaRoutes. From the clustering analysis, we selected 130 diverse target molecules to test our framework's ability to replicate and improve upon patent-derived “gold-standard” routes.

• Novelty assessment (ChemByDesign): to test generalization to unseen data, we selected 5 target molecules from ChemByDesign.³⁶ We specifically filtered for pathways published after 2021 with fewer than six steps. This filtering ensured that the ground-truth reactions were not present in our USPTO training data (which was last updated in September 2016), providing a test of zero-shot pathway prediction.

• Hybrid comparative set: to benchmark against existing hybrid planners, we selected specific complex targets (e.g., dronabinol, arformoterol, 4-ethenyl-2-fluorophenol) used in prior studies by Levin et al.⁴ and RetroBioCat.³¹ These molecules were manually curated to allow for a direct qualitative comparison of pathway length, cost, and enzyme usage.

2.2 Model development

We utilize the modern Hopfield network³⁴ based template prioritization architecture for our study. This architecture, initially introduced by Seidl et al.,²⁶ consists of three main components: a molecule encoder, a reaction template encoder, and one or more stacked or parallel Hopfield layers.

The molecule encoder function, denoted as h^m_w(·), learns a relevant representation for the input molecule m. We utilize a fingerprint-based approach, specifically the extended connectivity fingerprint (ECFP), coupled with a fully connected neural network with weights w. This encoder maps a molecule to a dense representation m_h = h^m_w(m) of dimension d_m.

Similarly, the reaction template encoder function, denoted as h^t_v(·), learns relevant representations of reaction templates. We employ a fully connected neural network with RDKit template fingerprints as input. The function is applied to the set of all templates , and the resulting vectors are concatenated column-wise into a matrix with shape (d_t, K), where K is the number of templates. Both encoders utilize Xavier Initialization⁴² and Batch Normalization⁴³ to improve convergence and prevent vanishing gradients.

The core of the model consists of Hopfield layers, denoted as g(·,·), which associate the molecule with the memory of templates. To perform the retrieval, the encoded molecule m_h and the template matrix T_h are projected into a common associative space to form the state pattern ξ and the stored patterns X, respectively. The Hopfield layer then updates the molecule representation via an attention-like mechanism. The update rule yielding the new state ξ_new is defined as:

ξ_new = Xp = X·softmax (βX^⊤ξ) (1)

where β is a learnable scaling parameter (inverse temperature), p is the vector of associations (probabilities) over the templates, and softmax is applied column-wise.

For the loss function, given a training pair (m, t) and the set of all templates , the model aims to maximize the probability assigned to the correct template t. We employ the negative log-likelihood as the loss function, optimizing parameters using stochastic gradient descent via the AdamW optimizer. To prevent overfitting, we apply dropout and L2 regularization. Additionally, a post-processing fingerprint-based substructure screen is used during inference to filter out chemically non-applicable templates.

We train two sets of models on our datasets: the synthetic template prioritizer (consisting of five models) and the enzymatic template prioritizer. We utilized PyTorch⁴⁴ version 1.9.0 and Python version 3.8 to implement the model. These models take a target molecule as input and rank the most applicable rules within the collected dataset. Given the same inputs, the outputs from the five synthetic template prioritizers are taken individually, and a final ranking of the templates is collated using the highest predicted score. Furthermore, to fine-tune the models for optimal performance, hyperparameter tuning was performed using a one-factor-at-a-time (OFAT) approach. We tuned number of epochs, concatenation threshold, dropout rate, learning rate and some hopfield parameters, more information on the parameters and the values chosen can be found in SI Sections 1–4. Each experiment involved modifying a single hyperparameter from a baseline configuration while keeping all other settings constant. The performance of each model was assessed using an evaluation score, with lower validation loss indicating better performance. Training and validation loss curves and accuracy metrics for the enzymatic and synthetic template prioritizers are provided in SI Fig. 1 and 2.

2.3 Searching algorithm

We employ a weighted global greedy tree search approach to explore synthesis pathways for a target molecule. This method leverages the template prioritizer models to identify the most applicable reaction rules. Specifically, the MHN generates a probability distribution over the predefined reaction templates, and we apply a cutoff by selecting only the top-k (also user-tunable) most confident applicable rules. These filtered rules are then applied to the target molecule to derive precursor molecules. In this search framework, precursor molecules are represented as nodes in a tree, while reactions and their associated conditions serve as edges. However, not all precursors may be readily available or affordable. We define “buyable” as being available for purchase at a cost of under $100 per g. A representative synthetic pathway tree illustrating the search output for 2-phenoxyethanamine is provided in SI Fig. 3. Algorithm 1 outlines the weighted global greedy tree search algorithm. Excluding API latency, the runtime complexity is , where b is the effective branching factor bounded by the user-tunable top-k parameter, d is the maximum search depth, and log V reflects the cost of priority queue operations over V nodes. The search algorithm proceeds as follows:

• Initialization: we begin by initializing the search with the target molecule as the start_node. This node is added to a priority queue with high priority.

• Global greedy search: the main search loop iterates until the priority queue is empty. At each iteration, the node with the highest priority is popped from the queue.

• Goal check: if the node meets the goal criteria (e.g., low cost or maximum depth), the algorithm continues to the next iteration.

• Rule application: for the current node, we retrieve the top-k transformation rules ranked by the MHN's predicted probabilities. We apply each valid rule to generate new chemical structures. The physical and operational properties of these new structures (such as cost, temperature, and solvent score) are then calculated.

• Node insertion: each new structure is encapsulated in a new_node, which is added to the current node's subtrees and inserted into the priority queue based on its multi-objective score.

• Scoring: nodes are strictly prioritized in the queue based on a user-tunable criterion P′ = f(C′, T′, S′) that includes precursor cost, reaction temperature, and solvent toxicity. The highest-scoring nodes are explored first.

• Termination: the search continues iteratively until buyable precursors are found, the tree is fully explored up to a specified maximum depth, or the allotted search time is exhausted.

2.4 Scoring system

We introduce a user-tunable scoring methodology to evaluate synthesis pathways. In this approach, individual scores for three features are obtained and summed based on user-tunable weights, allowing users to prioritize features according to their preferences during the pathway search. A detailed workflow diagram of our inference and scoring workflow is shown in Fig. 2, and representative ablation studies illustrating the effect of the user-tunable scoring weights on route selection are presented in SI Fig. 11–13. The first feature of the score is the cost. We utilize the ASKCOS⁴⁵ buyable dataset, Molport,⁴⁶ Mcule,⁴⁷ and Chemspace⁴⁸ APIs with the ChemPrice⁴⁹ library to determine the cost and availability of molecules. A molecule is considered buyable if it can be purchased for less than $100 per g.^4,50 Users have the flexibility to adjust this threshold based on their objectives. The score is calculated as , as we aim to maximize the score, thus exploring reactions involving cheaper precursors first. The value 500 is chosen as a normalizer because we define anything costing above $500 per g as non-buyable; the user can modify this threshold. These threshold align with the affordability criteria for research-grade building blocks utilized in major chemical databases.⁵⁰

Fig. 2 Inference and scoring workflow. This schematic illustrates the step-by-step pipeline executed during each iteration of the weighted global greedy tree search (Algorithm 1). First, the highest-scoring target node (N) is extracted from the global priority queue. The node is concurrently passed to the synthetic and enzymatic template prioritizers. For the synthetic pathway, an ensemble of 5 MHN models predicts applicable rules. These ranked lists are normalized, merged, and deduplicated (retaining the highest score per unique rule) to isolate the top n_syn templates. Simultaneously, the single enzymatic MHN model yields the top n_enz templates. During the applicability check, this combined set of n_syn + n_enz rules is applied to N to generate chemically valid precursor nodes (N′). Each valid precursor undergoes multi-objective scoring to calculate its specific cost (C′), temperature (T′), and solvent (S′) constraints, culminating in a final priority score (P′). Finally, a priority queue update inserts these newly scored precursors into the global queue, directing the next search iteration toward the most highly optimized pathways.

The second feature of the score is the reaction temperature. We employ a CASP tool developed by Gao et al.⁵¹ to predict the temperature at which a reaction might occur. This prediction is essential for exploring pathways where reactions may not have been previously documented or experimented on. To find the predicted temperature, we take the weighted average of the top ten predictions made by the tool, a strategy suggested by the authors to improve accuracy. The score is calculated as , as we aim to maximize the score, thus exploring reactions involving lower temperatures first. The value 300 is chosen as a normalizer because we define any reaction requiring over 300 °C as non-practical; this cutoff reflects the upper operating limits of standard synthetic laboratory equipment.⁵²

The third feature of the score is the solvent and reagent greenness-toxicity score. We use the same CASP tool developed by Gao et al.⁵¹ to predict the solvent and reagent required for a reaction. We curate a toxicity and greenness dataset, assigning a score of −1 for toxic molecules, 0 for neutral molecules, and +1 for green or natural molecules. We use the ACS Solvent Selection Guide⁵³ to classify 100 commonly used solvents. We also use the SuperNatural 3.0 dataset,⁵⁴ containing 350 000 natural products, and the T3DB dataset⁵⁵ to classify 4000 toxic molecules. This ternary scoring system serves as a heuristic to steer the tree search toward environmentally favorable regions of chemical space. While this approach does not currently account for process-specific metrics such as atom economy, E-factors, or reagent concentration, it provides a computationally efficient means of penalizing more hazardous precursors. This modular design allows the −1/0/+1 baseline to be substituted with more granular, continuous toxicity models as they become available. Furthermore, because these scores are additive, users can adjust the global weight of the greenness feature to balance environmental safety against precursor cost and temperature constraints.

3 Results and discussion

In this section, we evaluate the performance of the MHNpath framework in prioritizing reaction templates and generating feasible retrosynthetic pathways. We benchmark our predictive models against established baselines and analyze the practical utility of the proposed scoring system in guiding pathway discovery.

3.1 Model performance

To evaluate the effectiveness of our synthesis planning framework, we employed a comprehensive set of metrics designed to assess both the accuracy and applicability of the predicted reaction templates. Table 1 presents the performance of our model in comparison to existing state-of-the-art models, using three primary metrics:

Table 1 Performance metrics for template prioritization. The , , and results for the enzymatic dataset are color coded. The table is divided into three sets of metrics. The first set (columns 2–5) shows the accuracy of the presence of the literature rule in the top predictions (T1, T10, T50, T100). The second set (columns 6–9) represents the average number of applicable rules in the top predictions. The third set (columns 10–13) indicates the accuracy of the presence of at least one applicable rule in the top predictions

a The synthetic model is trained and tested on a separate, larger dataset, while the other three models are trained and tested on the same dataset.

---

1. Accuracy of the presence of the literature rule in the top n predictions (T1, T10, T50, T100): this metric evaluates the model's ability to prioritize reaction templates documented in the literature. The “literature rule” refers to the reaction template extracted from the corresponding reaction in our dataset. This metric is crucial for assessing how well the model replicates known synthetic pathways. This is the commonly used metric in recent works.^4,26 While this metric is important for evaluating the replication of known synthetic pathways, assessing the model's ability to propose novel and feasible pathways is equally critical. Therefore, we also consider the following metrics, which provide a more comprehensive view of the model's performance.

2. Average number of applicable rules in the top n predictions: this metric provides insight into the diversity and feasibility of the proposed reaction templates. A rule is considered “applicable” when the RDKit library successfully applies the transformation to the target molecule, resulting in valid precursor structures. This metric is particularly important as it indicates the model's capability to explore broader chemical space and suggest alternative synthetic routes that may not be present in the literature but are chemically plausible. This is critical for our tree search, as a higher density of valid rules in the top predictions (T50, T100) expands the branching factor of the search tree with high-quality and newer candidates.

3. Accuracy of the presence of at least one applicable rule in the top n predictions: this metric complements the other two by measuring the model's practical utility in retrosynthetic analysis. It assesses how often the model suggests at least one viable synthetic step, ensuring the progression of retrosynthetic analysis in chemically plausible directions. High scores here ensure the search algorithm rarely hits “dead ends”.

The introduction of these additional metrics (2 and 3) addresses the limitations of traditional evaluation methods, which often focus solely on predicting literature rules. While the accuracy of predicting literature rules is important, it does not fully capture the model's ability to suggest novel pathways. By incorporating metrics that consider rule applicability and diversity, we provide a more comprehensive evaluation that aligns with the exploratory nature of retrosynthetic planning.

Our enzymatic model (Ours (Enz)) demonstrates good performance across most metrics compared to the baseline DNN (Deep Neural Network)⁴ and MHN²⁶ models. It is important to contextualize the baseline performance, the DNN achieves a top-1 accuracy of only 10%, which underscores the inherent complexity of the retrosynthetic prediction task. Unlike standard classification problems, template prioritization involves selecting the correct chemical transformation from a massive search space of valid reaction templates. In this high-dimensional context, low absolute accuracy scores are standard, and incremental gains represent significant practical improvements.

First, compared to the standard DNN baseline,⁴ our model delivers a large improvement. The DNN achieves a top-1 literature accuracy of only 10% and a top-1 applicability presence of 12.8%. Our model nearly doubles these figures to 18.3% and 20.1%, respectively. Given the high-dimensional search space (selecting 1 out of >17 000 templates), this jump represents a significant improvement in the model's ability to prioritize relevant chemistry. This trend is mirrored in the applicability metrics, indicating that our model not only retrieves ground truth better but also generates a higher number of valid chemical precursors.

Second, while our architecture shares similarities with the MHN model,²⁶ our specific enhancements (Xavier initialization, rigorous dropout, and hyperparameter tuning) yield better convergence and robustness, particularly when predicting precursors for structurally complex or challenging target molecules. For example, in the avg. number of applicable rules metric at T50, our model outperforms the MHN baseline (6.125 vs. 5.647). This indicates that our model retrieves approximately 8.5% more valid chemical options in the top-50 candidates. This increased density of applicable rules provides the global greedy search algorithm with a richer pool of precursors, reducing the likelihood of missing a viable pathway.

Finally, the synthetic ensemble model (Ours (Syn)) demonstrates exceptional performance on the USPTO dataset, achieving a T1 accuracy of 42.2%. Most notably, it exhibits a massive applicability rate: the top-10 predictions yield, on average, 6.315 valid reaction rules. This high applicability ensures that the synthetic branch of our hybrid planner can almost always identify multiple feasible precursors, allowing the scoring system to aggressively optimize for cost and environmental impact without running out of chemical options.

It is important to note that the baseline DNN,⁴ the MHN model,²⁶ and our enzymatic model (Ours (Enz)) were trained and evaluated on the exact same enzymatic dataset splits, ensuring a direct and fair comparison.

3.2 Comparison with pathways from the literature

We assessed our synthesis planning framework by benchmarking it against established datasets and frameworks, including PaRoutes³⁵ and ChemByDesign.³⁶ Precisely, we assess the pathway lengths generated by our framework and report metrics on the number of replicated pathways and the average number of predicted pathways per molecule. We also present some representative examples involving a hybrid approach combining enzymatic and synthetic pathways in Fig. 3(a) and SI Fig. 4 and 5. This comparative analysis aimed to assess our framework's accuracy, efficiency, and versatility in predicting viable synthetic routes for complex target molecules.

Fig. 3 (a) Tree of reaction pathways. The tree shows a representative example for a pathway presented in PaRoutes.³⁵ Our predicted pathway is producible using cheap precursors and less toxic, naturally occurring solvents. (b) Performance metrics for literature comparison. These plots present the number of molecules solved, the average number of pathways predicted, and the distribution of predicted pathway lengths.

PaRoutes³⁵ is a robust benchmarking framework for evaluating multi-step retrosynthesis methods. It comprises two datasets of 10 000 synthetic routes derived from the patent literature alongside a curated list of purchasable molecules and reactions suitable for training retrosynthesis models. For this study, we utilized a scaffold split technique to select pathways associated with the ten most commonly occurring scaffolds in the PaRoutes dataset, which are given in SI Section 5. This approach enabled us to efficiently assess our framework's ability to predict synthetic routes that align with those documented in the patent literature while ensuring chemical diversity.

ChemByDesign,³⁶ on the other hand, is an online platform that organizes experimentally verified reaction pathways by name, year, and author. To ensure an unbiased evaluation, we focused on pathways discovered after 2021 that are under six steps long. These pathways were intentionally excluded from our training dataset to test the predictive capabilities of our framework on novel and unexplored synthetic routes. This strategy provided a unique opportunity to evaluate how effectively our framework generalizes to unseen data.

The results of this comparative analysis are summarized in Fig. 3. Fig. 3(a) presents a representative tree of reaction pathways generated by our framework for a target molecule in the PaRoutes dataset. This evaluation was conducted on an expanded test set comprising 130 molecules from the PaRoutes dataset and 5 novel targets from ChemByDesign. Notably, our predicted pathway utilizes inexpensive precursors and environment-friendly solvents, such as ethanol and methanol, demonstrating its potential for sustainable synthesis planning. The tree also highlights alternative routes that leverage naturally occurring molecules while avoiding toxic solvents like dichloromethane (DCM); further examples of these green alternatives are illustrated in SI Fig. 7–9. These results underscore the framework's ability to prioritize green chemistry principles without compromising synthetic feasibility and cost.

Fig. 3(b) provides quantitative performance metrics for our framework compared to PaRoutes and ChemByDesign. The first bar plot illustrates the number of molecules attempted and successfully solved by each method. Our framework solved 114 out of 135 molecules attempted. A representative example of a successfully replicated PaRoutes pathway, annotated with reaction conditions and precursor costs, is illustrated in SI Fig. 7. Furthermore, our framework replicated 91 known pathways from the literature while suggesting a high number of pathways on average per molecule (4.6 and 2.2). These results demonstrate its ability to explore diverse chemical spaces and propose multiple viable options for synthetic planning.

Finally, the third bar plot compares the lengths of predicted pathways relative to known literature routes. Our framework identified 55 shorter or 39 equivalent-length pathways, with 488 cases yielding routes two steps shorter than those documented in PaRoutes, including two notable molecules N-(3-[5-chloro-2-(difluoromethoxy)phenyl]-1-2-[(4-pyridinylmethyl)amino]ethyl-1H-pyrazol-4-yl)pyrazolo[1,5-a]pyrimidine-3-carboxamide and 4-(1-tert-butyl-4-oxo-5H-pyrazolo[4,3-c]pyridin-3-yl)thiophene-2-carboxamide, and one molecule (lupinine) yielding a route two steps shorter than that documented in ChemByDesign. This capability to optimize pathway length is particularly valuable in industrial settings where shorter synthetic routes can lead to significant cost savings and improved process efficiency.

SI Fig. 9 illustrates an alternative pathway for lupinine⁵⁶ (found in ChemByDesign) synthesis predicted by MHNpath. Our framework identified a streamlined three-step route with moderate reaction conditions (12.58 °C, 11.1 °C, and −23.81 °C) and low precursor costs ($87.36 per g and $0.10 per g), contrasting with the five-step approach by Wang et al.⁵⁶ that spans temperatures from −78 °C to 85 °C. The test set of utilized from PaRoutes has been provided in the SI Data.

3.3 Comparison with other models

We benchmark the performance of MHNpath against two retrosynthetic planning tools: RetroBioCat³¹ and the hybrid enzymatic-synthetic planner developed by Levin et al.,⁴ as illustrated in Fig. 4. Precisely, we assess the pathway lengths generated by our framework and report metrics on the number of replicated pathways and the average number of predicted pathways per compound. We also present some case studies involving a hybrid approach combining enzymatic and synthetic pathways in Fig. 4(i) and SI Fig. 6. These comparisons highlight the advantages of our framework in terms of pathway length, cost-effectiveness, and the ability to replicate and improve upon previously reported pathways.

Fig. 4 (a–g) Published pathways. Overview of existing pathways to produce dronabinol. (h) Previously predicted pathways by Levin et al.⁴ Four-step reaction pathway to produce dronabinol as predicted by Levin et al.⁴ (i) Our predicted pathway. Three-step reaction to produce dronabinol from cheap precursors in ambient temperatures. We also replicated some of the other pathways. (j) Performance metrics for comparison with other models. These plots present the number of molecules solved, the average number of pathways predicted, and the distribution of predicted pathway lengths.^4,31

Although synthetic planning tools such as AiZynthFinder,⁵⁷ ASKCOS,⁴⁵ Retro*,²⁹ and MHNreact²⁶ represent significant benchmarks in the field, this section prioritizes comparisons with hybrid planners. The predictive capabilities of the underlying MLP and MHN-based engines used by these synthetic tools have already been quantitatively evaluated in Section 3.1 (Table 1). We exclude a full tree search comparison with these platforms as they are restricted to purely synthetic routes and lack the integrated multi-objective scoring criteria central to MHNpath, specifically real-time cost, toxicity, and temperature optimization. Furthermore, we do not report search time metrics because our framework relies on live API calls to retrieve dynamic pricing information. Since the total search duration is dominated by network latency rather than algorithmic efficiency, direct runtime comparisons with tools utilizing static building block datasets would be unrepresentative. Consequently, comparing global success rates between conventional planners, which optimize primarily for route brevity and model confidence, and MHNpath would be fundamentally misaligned, as our algorithm routinely bypasses traditionally ‘successful’ shortest paths in favor of longer, but demonstrably greener and more cost-effective alternatives.

RetroBioCat³¹ is a widely used platform for designing biocatalytic cascades. RetroBioCat facilitates the construction of selective and efficient biocatalytic pathways by leveraging an expanding enzyme toolbox and encoded reaction rules. Its strength lies in its ability to identify promising enzyme-specific routes, validated through several literature examples. However, as shown in Fig. 4(i), our framework finds shorter pathways than RetroBioCat for many molecules like 4-ethenyl-2-fluorophenol and 2-phenylpiperidine. Specifically, MHNpath successfully solved all five molecules attempted using RetroBioCat data, with an average of 4.2 pathways per molecule. The molecules attempted are listed in SI Section 18. Moreover, our framework excelled in uncovering multiple shorter pathways and various pathways of equivalent length, offering users a range of options encompassing different enzymes, starting molecules, and temperature conditions.

SI Fig. 6 compares the pathways generated by RetroBioCat and MHNpath. In the top panel, while the RetroBioCat pathway for a fluorinated compound involves multiple enzymatic steps (TPL, TAL, and DC), our approach achieves the same transformation using a single enzymatic step with TPT at 99 °C, thereby significantly reducing both the complexity and the precursor cost to $6.03 per g. In the bottom panel, for an amine compound, the RetroBioCat route requires a three-step enzymatic cascade (CAR, TA, and IRED) that depends on cofactors such as NADPH and ATP. In contrast, our hybrid pathway utilizes only two enzymatic steps (PT at 54 °C and PP at 10.1 °C), with considerably lower precursor costs ($0.1 per g and $1.8 per g, respectively).

The hybrid enzymatic-synthetic planner by Levin et al.⁴ represents another state-of-the-art approach that combines neural networks trained on extensive reaction databases with enzymatic transformations. Levin's model offers hybrid strategies for complex molecules such as THC and R,R-formoterol, showcasing the potential of enzyme-synthetic integration. However, our results demonstrate that MHNpath not only replicates several of Levin's proposed pathways but also identifies novel alternatives that are shorter and more cost-effective.

As shown in Fig. 4(h) and (i), Levin's predicted pathway for dronabinol involves a four-step synthesis. In contrast, our framework proposes a three-step pathway that reduces synthesis costs to $0.12 per g at ambient temperatures. This significant cost reduction is achieved by incorporating an optimized enzymatic step that eliminates the need for high-temperature reactions, underscoring the practical advantages of our approach. Although the final step of the predicted pathways remains unchanged, our framework excels in exploring diverse routes to the penultimate molecule, owing to its global-greedy tree search algorithm and scoring methodology. In addition to replicating published pathways from Levin et al., our framework identifies four shorter additional pathways.

Furthermore, SI Fig. 8 presents a novel four-step synthetic pathway for the synthesis of arformoterol predicted by our model. It starts from inexpensive precursors ($3.11 per g, $1.30 per g, $0.51 per g, and $0.24 per g). In contrast, Levin et al.’s⁴ approach involves a more complex five-step biocatalytic cascade that requires multiple enzymes and cofactors. The detailed annotations of reaction conditions in our pathway, such as temperature, solvent, and reagent specifics, demonstrate the practical applicability of our hybrid strategy.

A key advantage of MHNpath is its ability to balance exploration and exploitation during retrosynthetic planning. This is evident from the higher average number of pathways generated per molecule (6.5 and 4.2), which provides chemists with a broader range of options for optimizing synthesis strategies based on reaction temperature, cost, or solvent toxicity.

4 Conclusions

The development and implementation of the MHNpath tool represent progress over existing retrosynthetic frameworks in computer-aided synthesis planning. By leveraging machine learning, we have created a framework that can efficiently predict retrosynthetic pathways and facilitate the exploration of diverse chemical spaces. Our template prioritizer outperformed existing methods, as demonstrated by its higher accuracies across specific metrics and increased number of applicable rules. Importantly, our framework was able to replicate gold-standard pathways from PaRoutes as well as novel experimental reaction routes reported in ChemByDesign. In some instances, the model not only reproduced known synthesis routes but also identified alternative pathways that were shorter and more cost-effective. For example, our approach discovered a three-step synthetic route for dronabinol, reducing the synthesis cost to $0.12 per g as compared to a previously reported four-step pathway. The framework also demonstrated partial retrosynthetic decomposition of more structurally complex targets, as illustrated in SI Fig. 10, though reaching fully buyable terminal states for such targets remains an open challenge.

MHNpath underscores how ML can be utilized in reducing the manual workload and minimizing the trial-and-error strategies traditionally associated with chemical synthesis. Combined with a tree search-based strategy, the machine learning model allows for the automated and efficient prediction of optimal reaction pathways. The tool's user-tunable criteria allow researchers to prioritize different aspects of the synthesis process, such as cost, reaction temperature, and toxicity of solvents and reagents. This adaptability is particularly beneficial for experimental chemists who can tailor the tool to meet their research needs and constraints.

Recent years have seen a surge in chemistry-focused large language models (LLMs) such as BatGPT-Chem⁵⁸ and ChemDFM,⁵⁹ which bypass rigid template libraries to offer broad generalization across diverse chemical spaces. However, several properties of these models make them poorly suited as the core engine of a constraint-driven CASP tool. First, their non-deterministic, autoregressive generation processes introduce a non-negligible risk of producing chemically invalid intermediates or synthetically implausible steps, errors that propagate and compound across iterations in a deeply branched tree search. Second, state-of-the-art chemistry LLMs typically require on the order of 13 to 15 billion parameters, creating substantial computational bottlenecks when integrated into iterative search algorithms that may query the model hundreds of times per target molecule. Third, and most critically, it remains exceptionally difficult to algorithmically enforce hard physical constraints, such as strict precursor cost ceilings or solvent toxicity indices, through text-based prompting alone.

By contrast, the MHN-based template prioritization approach used in MHNpath explicitly maps target molecules and the full template library into a shared associative memory space, producing deterministic and structurally valid outputs at every retrosynthetic step by construction. Although MHNpath operates on a predefined template library, it is not strictly limited to historical precedent. The extraction of 301 242 localized transformation rules from approximately 1.7 million USPTO reactions allows the model to apply generic reaction templates to novel substrates, effectively generating novel multi-step pathways without sacrificing intermediate chemical validity. The multi-objective scoring function P′ = f(C′, T′, S′) enforces physical constraints algorithmically and with full transparency, granting synthetic chemists a degree of precise operational control that generative LLMs currently struggle to guarantee.

Beyond architectural differences, MHNpath's utility is designed to function as a decision-support system rather than a black-box generator. While the multi-objective optimization identifies pathways that are computationally “greener” or more “cost-effective,” these outputs serve as prioritized hypotheses for expert evaluation. By quantifying the trade-offs between temperature, cost, and solvent toxicity, the framework provides synthetic chemists with a structured “shortlist” of routes that align with specific laboratory constraints, significantly narrowing the vast search space before high-cost experimental validation is attempted.

Despite its capabilities, MHNpath has some limitations that outline future work directions. One potential issue is the reliance on predicted pathways without sufficient experimental validation, meaning that while the tool can suggest plausible synthetic routes, it cannot guarantee their success in practice and should be used as a starting point for further experimental investigation rather than a definitive solution. Additionally, the tool cannot effectively address enantiomer selectivity issues, as it does not inherently account for the stereoselective outcomes of reactions involving chiral molecules. Furthermore, MHNpath is limited in its ability to provide detailed mechanistic insights into the predicted reactions, a drawback for users who require a deeper understanding of the underlying processes to optimize and troubleshoot reactions. Its overall accuracy and reliability are also heavily dependent on the quality and comprehensiveness of the underlying datasets; incomplete or biased data can lead to inaccurate predictions and missed opportunities for novel synthetic routes. To address these challenges, future work will focus on integrating enantioselective predictions by incorporating enantiomer-specific data and developing algorithms capable of predicting stereoselective outcomes, providing more mechanistic information through the integration of mechanistic databases, and continuously expanding and diversifying the datasets via collaborations and feasibility tests with experimental chemists.

In summary, MHNpath demonstrates the potential of ML-driven tools in advancing CASP. Our results highlight this promise: for example, the enzymatic model achieved a top-1 accuracy of 18.3%, a marked improvement over the 10% of the DNN baseline and comparable to the 18.1% of MHN. On the other hand, the synthetic ensemble model reached impressive 42.2% top-1 accuracy. Moreover, our framework delivered an average of approximately 6 applicable rules in the top 10 predictions, highlighting its ability to explore diverse chemical spaces. In comparisons with pathways in the literature, MHNpath successfully solved 114 out of 135 molecules attempted and even identified a three-step synthetic route for dronabinol that reduces costs to $0.12 per g at ambient temperatures, outperforming a competing four-step pathway. However, addressing its limitations and mitigating potential misuse remains essential to ensure the long-term success and reliability of MHNpath for assisting in synthesis of complex organic molecules.

Code availability

The code used for data processing, model training, inference and the instructions to run our framework are available at https://github.com/MSRG/mhnpath.

Conflicts of interest

There are no conflicts to declare.

Data availability

The synthetic dataset utilized in this study was obtained from the USPTO⁴⁰ repository, while the enzymatic dataset was sourced from the BKMS³⁷ and RHEA³⁸ databases. All datasets are publicly accessible and open-source. The processed datasets generated and analyzed during the study and the trained model weights are available on Figshare (https://doi.org/10.6084/m9.figshare.28673540).

Supplementary information (SI) is available. See DOI: https://doi.org/10.1039/d5dd00562k.

Acknowledgements

This research has received funding from the research project: “Quantum Software Consortium: Exploring Distributed Quantum Solutions for Canada” (QSC). QSC is financed under the National Sciences and Engineering Research Council of Canada (NSERC) Alliance Consortia Quantum Grants #ALLRP587590-23. The authors also acknowledge the funding received from the University of Toronto’s Data Sciences Institute (DSI) via its Catalyst Grant as well as the Canada Research Coordinating Committee's (CRCC) New Frontiers in Research Fund (NFRF) for their continued support. S. P. and V. K. P. thank DSI for providing financial support during the summer of 2025 via the Summer Undergraduate Data Science Research Opportunities Program and the Postdoctoral Fellowship, respectively. S. P. would also like to thank MolPort, Mcule and ChemSpace for providing API keys and access to a cost related database of chemical compounds. V. K. P. is grateful for the computational resource support provided by the Digital Research Alliance of Canada.

References

1. S. Szymkuć, et al., Computer-assisted synthetic planning: The end of the beginning, Angew. Chem., Int. Ed., 2016, 55, 5904–5937.
2. D. A. Pensak and E. J. Corey, Lhasa—logic and heuristics applied to synthetic analysis, in Computer-Assisted Organic Synthesis, vol. 61 of ACS Symposium Series, American Chemical Society, 1977, pp. 1–32,  DOI:10.1021/bk-1977-0061.ch001.
3. C. Avila, et al., Chemistry in a graph: modern insights into commercial organic synthesis planning, Digital Discovery, 2024, 3, 1682–1694,  10.1039/D4DD00120F.
4. I. Levin, M. Liu, C. A. Voigt and C. W. Coley, Merging enzymatic and synthetic chemistry with computational synthesis planning, Nat. Commun., 2022, 13, 7747.
5. K. C. Nicolaou, T. Montagnon, Molecules that Changed the World, Wiley-VCH, Weinheim, 2008.
6. D. G. Brown and J. Boström, Analysis of past and present synthetic methodologies on medicinal chemistry: Where have all the new reactions gone?, J. Med. Chem., 2016, 59, 4443–4458.
7. P. E. Hart, N. J. Nilsson and B. Raphael, A formal basis for the heuristic determination of minimum cost paths, IEEE Trans. Syst. Sci. Cybern., 1968, 4, 100–107.
8. P. Schwaller, et al., Molecular transformer: A model for uncertainty-calibrated chemical reaction prediction, Chem. Sci., 2019, 10, 369–377.
9. D. Kreutter, P. Schwaller and J.-L. Reymond, Predicting enzymatic reactions with a molecular transformer, Chem. Sci., 2021, 12, 8648–8659,  10.1039/D1SC02362D.
10. I. V. Tetko, P. Karpov, R. V. Deursen and G. Godin, State-of-the-art augmented nlp transformer models for direct and single-step retrosynthesis, Nat. Commun., 2020, 11, 5575,  DOI:10.1038/s41467-020-19266-y.
11. R. Irwin, S. Dimitriadis, J. He and E. J. Bjerrum, Chemformer: a pre-trained transformer for computational chemistry, Mach. Learn. Sci. Technol., 2022, 3, 015022,  DOI:10.1088/2632-2153/ac3ffb.
12. E. Granqvist, R. Mercado and S. Genheden, Retrosynformer: planning multi-step chemical synthesis routes via a decision transformer, Digital Discovery, 2026, 5, 348–362,  10.1039/D5DD00153F.
13. X. Wang, et al., Retroprime: A diverse, plausible and transformer-based method for single-step retrosynthesis predictions, Chem. Eng. J., 2021, 420, 129845 ), https://www.sciencedirect.com/science/article/pii/S1385894721014303.
14. M. Andronov, N. Andronova, M. Wand, J. Schmidhuber and D.-A. Clevert Fast and scalable retrosynthetic planning with a transformer neural network and speculative beam search, arXiv, 2025, preprint, arXiv:2508.01459,  DOI:10.48550/arXiv.2508.01459.
15. J. Bradshaw, M. J. Kusner, B. Paige, M. H. S. Segler and J. M. Hernández-Lobato, A generative model for electron paths, in International Conference on Learning Representations, 2019, https://openreview.net/forum?id=r1x4BnCqKX.
16. M. Sacha, et al., Molecule edit graph attention network: Modeling chemical reactions as sequences of graph edits, J. Chem. Inf. Model., 2021, 61, 3273–3284,  DOI:10.1021/acs.jcim.1c00537.
17. C. Coley, et al., A graph-convolutional neural network model for the prediction of chemical reactivity, Chem. Sci., 2019, 10, 370–377,  10.1039/C8SC04228D.
18. W. W. Qian et al. Integrating Deep Neural Networks and Symbolic Inference for Organic Reactivity Prediction ( 2020).
19. H. Bi et al., Non-autoregressive electron redistribution modeling for reaction prediction, in Proceedings of the 38th International Conference on Machine Learning, vol. 139 of Proceedings of Machine Learning Research, ed. M. Meila and T. Zhang, PMLR, 2021, pp. 904–913, https://proceedings.mlr.press/v139/bi21a.html.
20. E. Corey, A. Long and S. Rubenstein, Computer-assisted analysis in organic synthesis, Science, 1985, 228, 408–418.
21. M. H. S. Segler and M. P. Waller, Neural-symbolic machine learning for retrosynthesis and reaction prediction, Chem.–Eur. J., 2017, 23, 5966–5971,  DOI:10.1002/chem.201605499.
22. C. W. Coley, R. Barzilay, T. S. Jaakkola, W. H. Green and K. F. Jensen, Prediction of organic reaction outcomes using machine learning, ACS Cent. Sci., 2017, 3, 434–443,  DOI:10.1021/acscentsci.7b00064.
23. X. Zhang, H. Lin, M. Zhang, Y. Zhou and J. Ma, A data-driven group retrosynthesis planning model inspired by neurosymbolic programming, Nat. Commun., 2025, 16, 192.
24. J. Roh, et al., Higher-level strategies for computer-aided retrosynthesis, ChemRxiv, 2025, preprint,  DOI:10.26434/chemrxiv-2025-21zvt-v3.
25. S. Chen and Y. Jung, A generalized-template-based graph neural network for accurate organic reactivity prediction, Nat. Mach. Intell., 2022, 4, 772–780,  DOI:10.1038/s42256-022-00526-z.
26. P. Seidl, et al., Improving few-and zero-shot reaction template prediction using modern hopfield networks, J. Chem. Inf. Model., 2022, 62, 2111–2120.
27. J. Choe, H. Kim, Y. T. Chok, M. Gim and J. Kang, Retrosynthetic crosstalk between single-step reaction and multi-step planning, J. Cheminf., 2025, 17, 130.
28. X. Wang, et al., Towards efficient discovery of green synthetic pathways with monte carlo tree search and reinforcement learning, Chem. Sci., 2020, 11, 10959–10972,  10.1039/D0SC04184J.
29. B. Chen, C. Li, H. Dai and L. Song, Retro*: Learning retrosynthetic planning with neural guided a* search, in The 37th International Conference on Machine Learning (ICML 2020), 2020.
30. K. Yu et al., Double-ended synthesis planning with goal-constrained bidirectional search, arXiv, 2024, preprint, arXiv:2407.06334,  DOI:10.48550/arXiv.2407.06334.
31. W. Finnigan, L. J. Hepworth, S. L. Flitsch and N. J. Turner, Retrobiocat as a computer-aided synthesis planning tool for biocatalytic reactions and cascades, Nat. Catal., 2021, 4, 98–104.
32. R. A. Sheldon, The e factor: fifteen years on, Green Chem., 2007, 9, 1273–1283,  10.1039/B713736M.
33. I. T. Horváth and P. T. Anastas, Innovations and green chemistry, Chem. Rev., 2007, 107, 2169–2173,  DOI:10.1021/cr078380v.
34. H. Ramsauer et al., Hopfield networks is all you need, arXiv, 2020, preprint, arXiv:2008.02217,  DOI:10.48550/arXiv.2008.02217.
35. S. Genheden and E. Bjerrum, Paroutes: towards a framework for benchmarking retrosynthesis route predictions, Digital Discovery, 2022, 1, 527–539,  10.1039/D2DD00015F.
36. C. Draghici and J. T. Njardarson, Chemistry by design: A web-based educational flashcard for exploring synthetic organic chemistry, J. Chem. Educ., 2012, 89, 1080–1082,  DOI:10.1021/ed2006423.
37. M. Lang, M. Stelzer and D. Schomburg, Bkm-react, an integrated biochemical reaction database, BMC Biochem., 2011, 12, 1–9.
38. P. Bansal, et al., Rhea, the reaction knowledgebase in 2022, Nucleic Acids Res., 2021, 50, D693–D700,  DOI:10.1093/nar/gkab1016 , https://academic.oup.com/nar/article-pdf/50/D1/D693/42058388/gkab1016.pdf.
39. G. Landrum., Rdkit: Open-source cheminformatics software, 2016, https://github.com/rdkit/rdkit/releases/tag/Release_2016_09_4.
40. D. Lowe, Chemical reactions from US patents (1976-Sep2016), 2017,  DOI:10.6084/m9.figshare.5104873, https://figshare.com/articles/dataset/Chemical_reactions_from_US_patents_1976-Sep2016_/5104873.
41. H. Mary et al. datamol-io/datamol: 0.12.3, 2024,  DOI:10.5281/zenodo.10535844.
42. X. Glorot and Y. Bengio, Understanding the difficulty of training deep feedforward neural networks, in Proceedings of the Thirteenth International Conference on Artificial Intelligence and Statistics, vol. 9 of Proceedings of Machine Learning Research, ed. Y. W. Teh and M. Titterington, PMLR, Chia Laguna Resort, Sardinia, Italy, 2010, pp. 249–256, https://proceedings.mlr.press/v9/glorot10a.html.
43. S. Ioffe and C. Szegedy, Batch normalization: Accelerating deep network training by reducing internal covariate shift, in International conference on machine learning, PMLR, 2015, pp. 448–456.
44. A. Paszke et al., PyTorch: an imperative style, high-performance deep learning library, Curran Associates Inc., Red Hook, NY, USA, 2019.
45. C. W. Coley, et al., A robotic platform for flow synthesis of organic compounds informed by ai planning, Science, 2019, 365, eaax1566.
46. Molport, Easy compound ordering service - molport, https://www.molport.com/shop/index.
47. R. Kiss, M. Sandor and F. A. Szalai, http://mcule.com: a public web service for drug discovery, J. Cheminf., 2012, 4, P17,  DOI:10.1186/1758-2946-4-S1-P17.
48. ChemSpace, ChemSpace, https://chem-space.com, accessed: 2024-02-21.
49. M. C. Sorkun, B. Saliou and S. Er, Chemprice, a python package for automated chemical price search, ChemRxiv, 2024, preprint,  DOI:10.26434/chemrxiv-2024-1bxgg.
50. T. Sterling and J. J. Irwin, Zinc 15 – ligand discovery for everyone, J. Chem. Inf. Model., 2015, 55, 2324–2337,  DOI:10.1021/acs.jcim.5b00559.
51. H. Gao, et al., Using machine learning to predict suitable conditions for organic reactions, ACS Cent. Sci., 2018, 4, 1465–1476.
52. T. Razzaq and C. Kappe, Continuous flow organic synthesis under high-temperature/pressure conditions, Chem.–Asian J., 2010, 5, 1274–1289,  DOI:10.1002/asia.201000010.
53. ACS Green Chemistry Institute® Pharmaceutical Roundtable. Solvent selection guide: Version 2.0, 2011, Retrieved 12th May 2024 from: https://www.acs.org/content/acs/en/greenchemistry/research-innovation/tools-for-green-chemistry.html.
54. K. Gallo, et al., SuperNatural 3.0—a database of natural products and natural product-based derivatives, Nucleic Acids Res., 2022, 51, D654–D659,  DOI:10.1093/nar/gkac1008 , https://academic.oup.com/nar/article-pdf/51/D1/D654/48440479/gkac1008.pdf.
55. D. Wishart, et al., T3DB: the toxic exposome database, Nucleic Acids Res., 2014, 43, D928–D934,  DOI:10.1093/nar/gku1004 . https://academic.oup.com/nar/article-pdf/43/D1/D928/7311219/gku1004.pdf.
56. J. Wang, et al., Enantioselective synthesis of the 1, 3-dienyl-5-alkyl-6-oxy motif: Method development and total synthesis, Angew. Chem., Int. Ed., 2024, 63, e202400478.
57. S. Genheden, et al., Aizynthfinder: a fast, robust and flexible open-source software for retrosynthetic planning, J. Cheminf., 2020, 12, 70,  DOI:10.1186/s13321-020-00472-1.
58. Y. Yang, et al., Batgpt-chem: A foundation large model for chemical engineering, Research, 2025, 8, 0827,  DOI:10.34133/research.0827.
59. Z. Zhao, et al., Developing chemdfm as a large language foundation model for chemistry, Cell Rep. Phys. Sci., 2025, 6, 102523 ). https://www.sciencedirect.com/science/article/pii/S2666386425001225.

---