FSscore: A Personalized Machine Learning-based Synthetic Feasibility Score

The application allows a user to perform all steps necessary to get a personalized model as intended. These steps involve: pairing molecules, labeling pairs, fine-tuning the pre-trained model, and ultimately scoring molecules. The application is built using Streamlit ¹ but is deployed locally. This means that the speed is limited by the computer on which the user is running the app. All models were trained on one GPU (NVIDIA GeForce RTX 3090) but the application will run also on a computer without GPUs at the cost of a much increased running time. Furthermore, running it locally also requires the user to install the package, which we recommend doing using a conda environment and pip. All the necessary steps and usability are outlined below.
1. Installation
We recommend running the app with the packages stored in a conda environment¹¹1Conda installation guide: https://conda.io/projects/conda/en/latest/user-guide/install/index.html. First, clone the GitHub repository from the link above (alternatively download the zip file):

Next, change the directory to the repository and install all required packages, like so:

2. Start the application
Starting the application only requires one line of code (from the root directory of the repository), which will lead to a window to pop up or return an active link:

From the landing page, select Label molecules in the upper drop-down menu on the left. This will lead to selecting the next steps. If a dataset with labels is available, go to step 4., otherwise continue with the next step. To directly score molecules and a fine-tuned model is available one can jump to step 6.
3. Pair molecules
In this step, a list of molecules provided as SMILES first are paired, and then ranked based on uncertainty of the pre-trained model. This second step will take some time if not performed on a machine with GPU-support. Select the option on the left and provide the path to the csv file (Figure S1). This file should have smiles as a column header for the SMILES structures to be paired.

4. Label molecule pairs
If a dataset with paired molecules is available, place it in streamlit_app/data/unlabeled. The dataset should appear in the drop-down menu on the right. Make sure to provide the paired molecules with column headers smiles_i and smiles_j.
The labeling window (see Figure S2) prompts the user to select the molecule that is harder to synthesize, before going to the next pair. Once the user is finished labeling, the button Send labels can be pressed. The labeling process can also be continued later and the labeled data can be downloaded. Note: The more labeled pairs, the better the performance of the fine-tuned model.

5. Fine-tune
Once data is labeled, one can proceed to fine-tuning. This is by default done using the best pre-trained model from this work. The user can also place their own checkpoint at the same location (folder models). How to pre-train a new model is described on our GitHub repository. Before initializing fine-tuning, one can tune several hyperparameters. If nothing is adjusted, the (recommended) default parameters as used in this work are selected (Figure S3). Fine-tuning in general is very fast, but depending on your system this can take a few minutes. Proceeding to score molecules can be achieved by selecting the respective action on the upper drop-down menu. This process can cause some issues in which case we would recommend returning to the home page before trying again.

6. Score molecules
Lastly, one can deploy the previously fine-tuned model and score molecules. For this, the user may place a csv file in streamlit_app/data/scoring, which contains the column header smiles. Then, the user can select the correct dataset in the drop-down menu on the left as well as the desired fine-tuned model (Figure S4).

Our approach to learning a continuous score to assess the synthesizability is inspired by Choung \latinet al. ², which framed a similar task as a ranking problem using binary preferences. Specifically, every data point consists of two molecules for each of which we predict a scalar in separate forward passes. The minimization of the binary cross entropy between the true preference and the learned score difference $\delta_{ij}:=\hat{f}(m_{i})-\hat{f}(m_{j})$ (scaled using a sigmoid function) constitutes our training objective.

The function $f:\mathcal{M}\rightarrow\mathbb{R}$ learns to parametrize molecules as an expressive latent representation given a set of molecules $m_{1},...,m_{n}\in\mathcal{M}$. We represent molecules as graphs $G=~{}(\mathcal{V},\mathcal{E})$ with atoms as nodes $x_{1},...,x_{n}\in\mathcal{V}$ and bonds as edges $e_{1},...,e_{m}\in\mathcal{E}$ from which we can compute the line graph $L(G)$ offline iteratively to the desired depth. The transformation process to the line graph is defined such that the edges $\{e_{1},...,e_{m}\}$ of graph $G$ are the nodes of the line graph and these nodes are connected if the corresponding edges of $G$ share a node. The graph neural network (GNN) embedding the molecular graph consists of the graph attention network (GATv2) ³ operating on the graph $G$ and the Line Evolution (LineEvo) ⁴ layer operating on the line graph $L(G)$ as message passing schemes. Both GATv2 and LineEvo use an attention mechanism to update the node representations $\{h_{1},...,h_{n}\}$ as follows:

where $a\in\mathbb{R}^{2d^{\prime}}$ and $W\in\mathbb{R}^{d^{\prime}\times d}$ are learned, $\sigma$ denotes an activation function and $||$ denotes concatenation. The GNN operates on a hidden size of 128 and the input dimensions depend on the initial featurization of the edges and nodes (see Tab. S1. GATv2 layers use 8 heads and LeakyReLU as activation function after updating the hidden feature vector while LineEvo layers use ELU. In GATv2, these local attention scores $e_{ij}$ are then averaged across all neighbors $\mathcal{N}$ (see Eq. S2) to obtain a normalized attention coefficient $\alpha_{ij}$ to compute the updated node representation $h_{i}^{\prime}$ as a weighted average (see Eq. S3):

with PReLU ⁵ as nonlinearity $\sigma$. In LineEvo layers, $e(h_{i},h_{j})$ is simply transformed using ELU to obtain the new node representation $h_{i}^{\prime}$.

These transformation layers are stacked so that two GATv2 (G) layers are followed by one LineEvo (L) layer (GGLGGL). Each of these layers is followed by a readout function, which consists of a global max pooling layer and global weighted-add pooling as suggested by Ren \latinet al. ⁴ obtaining a molecular representation through concatenation of the two readouts. The intermediate molecular representations of all layers are summed, and the final score $s_{i}=\hat{f}(m_{i})$ is inferred with a multilayer perceptron (MLP) (3 hidden layers of size 256, ReLU as activation function). The model described above was compared to five other implementations: GATv2 layers only (GGG) and four fingerprint implementations, namely Morgan (boolean), Morgan counts, Morgan chiral (boolean), and Morgan chiral counts, all with radius 4 and 2048 bits. ⁶ The fingerprints are embedded by one linear layer (hidden size of 256) with ReLU as activation function followed by the aforementioned MLP. More detailed information can be found in Section S2.3.

To focus our scoring model on a desired chemical space, we apply a fine-tuning approach inspired by Choung \latinet al. ². This approach was optimized to require few data points in order to limit the time required from expert chemists labeling the training examples in case of labels based on human feedback. Datasets to be used for fine-tuning can be of various origins such as a chemical scope suggested by experimentalists or specific chemical spaces encompassing e.g. natural products. Training our model requires pairs of molecules with a binary preference as label. In order to do so, we first cluster a dataset using the k-means algorithm based on the Tanimoto distances of Morgan counts fingerprint. The number of clusters k was determined for every dataset individually where the mean Silhouette coefficient over a range of k was the smallest. ⁷. The range queried was set from 5 to 29 for datasets larger than 500 pairs and to 3 to 9 for smaller datasets. Secondly, molecules are paired in such a way that they come from different clusters, appear only once in the paired dataset and have opposite labels if available. We also investigated the influence of having overlapping pairs (molecules can appear multiple times) for the SYBA (CP and MC) and the meanComplexity datasets. This approach results in more fine-tuning data but does not satisfy our desire to keep the dataset size small since human labeling could not make use of this advantage. Subsequently, the uncertainty of the prediction of the pre-trained model on those pairs is determined based on the variance of $\delta_{ij}$ obtained using the Monte Carlo dropout method ⁸ with a dropout rate of 0.2 on 100 predictions. The dataset is sorted with descending variance so that the pairs with high uncertainty are labeled first or used first for fine-tuning when selecting a subset. If no label is available the top n pairs (n depends on the dataset size) are submitted to be evaluated by our expert chemist based on their preference with regard to synthesizability.

To center the obtained score around zero a regularization factor of 1e-4 was applied to the predicted score like to obtain the regularization loss $\mathcal{L}_{reg}(\hat{f};\lambda):=\lambda||\hat{f}||^{2}$ as suggested by Choung \latinet al. ². $\mathcal{L}_{reg}$ is added to the cross-entropy loss term as described in Section S2.1 to obtain the final loss $\mathcal{L}=\mathcal{L}_{CE}+\mathcal{L}_{reg}$.

Pre-training: The training set was randomly split in 25 equally sized subsets and the model was trained with each subset for 10 epochs totaling 250 epochs with a batch size of 128. This sequential learning approach was found to work well in increasing speed of training. The model was trained with the Adam optimizer ⁹ and an initial learning rate of 3e-4.

Fine-tuning: We chose a transfer learning type of fine-tuning where all the weights of the pre-trained model can get adapted. All the hyperparameter except learning rate and batch size are kept the same. In order to avoid the forgetting the previously learned we track the performance metrics on a random subset of 5,000 pairs from the hold-out test set of the pre-training data. Of the four tested initial learning rates (1e-5, 1e-4, 3e-4, 1e-3) we chose 1e-4 for the graph-based version and 3e-4 for the fingerprint-based models. These learning rates are a good trade-off between degradation of the previously learned and learning on the new dataset as can be seen in the learning curves in Appendix S6.1. We observed quick conversion with fewer than 20 epochs and in order to balance the aforementioned trade-off a custom early stopping method was applied. Training is stopped once either the validation loss (or training loss if the full dataset is used for production) has increased for 3 epochs (patience of 3) or the accuracy on the 5k subset of the hold-out test set decreased for 3 epochs with a delta threshold of 0.02. We trained for a maximum of 20 epochs. For evaluating the efficiency, we further varied the number of data points used for training of which we took 5 random pairs for validation. The results in Table S2 show improvement in performance the bigger the dataset size. The batch size was set to 4.

To pre-train our model on a large collection of reactions, we combine the USPTO_full ¹¹ patent dataset with a complementary dataset ¹². Reaction data implicitly contains information on the synthetic feasibility through the relation of reactant to product, with the product being synthetically more difficult. The USPTO_full dataset was downloaded according to https://github.com/coleygroup/Graph2SMILES/blob/main/scripts/download_raw_data.py ¹³ using all USPTO_full subsets and "src" as product and "tgt" as reactants. The dataset was split so that every data point is a pair of one reactant and its respective product and the SMILES strings were canonicalized. The datasets do not contain reactions with multiple products. Data points with empty strings, identical reactants and products (isomeric SMILES), either reactant or product containing less than four heavy atoms or containing element types not in {H, B, C, N, O, F, Si, P, S, Cl, Se, Br, I} were removed. The datasets were further deduplicated retaining one instance of replicates. Furthermore, data points were removed ensuring that there are no cycles in the reaction network. This was achieved by removing back edges with the Depth-First Search (DFS) algorithm as implemented by Sun \latinet al. ¹⁴. The filtered dataset consisted of 5,340,704 data points. The train test split was performed so that no molecules (note not just reactant:product pairs) are overlapping resulting in a significant data loss to prevent data leakage. The training set consists of 3,349,455 pairs and the hold out test set of 711,550 pairs (17.5%).

To qualitatively evaluate the pre-trained model the performance on the MOSES ¹⁵ and COCONUT ¹⁶ dataset were assessed. We downloaded the MOSES test set from  Polykovskiy \latinet al. ¹⁵ as is totaling 176,074 SMILES. COCONUT was extracted from  Sorokina \latinet al. ¹⁶, all SMILES returning None with RDKit ¹⁰ were removed and a random subset of 176,000 molecules was selected to match the MOSES fraction.

The data for the fine-tuning case studies come from various sources. The specific processing steps for every dataset is outlined below.

Chirality: The obtain a test set for the chirality case study we filtered the deduplicated molecules from the pre-training training set keeping only those with @ or @@ tokens in the SMILES string, which define the chirality of a tetrahedral stereocenter. Of those, we sample 1,000 SMILES randomly and obtain their partner by getting the non-isomeric SMILES using RDKit rdk ¹⁰.

SYBA sets: The CP and MC test sets were extracted from Voršilák \latinet al. ¹⁷. No further processing was required.

meanComplexity: The meanComplexity dataset was downloaded from Sheridan \latinet al. ¹⁸ and cleaned by removing SMILES whose conversion to an RDKit molecule object failed. This processing step removed 44 SMILES yielding a fine-tuning test set of 1,731 molecules.

PROTACs: For the PROTAC case study we extracted all PROTAC SMILES from the PROTAC-DB Weng \latinet al. ¹⁹ and deduplicated the dataset keeping the first instance resulting in 3,270 SMILES. After pairing and ranking by confidence as described in Appendix S2.2, 100 pairs were labeled by a medicinal chemist with expertise in PROTACs. In order to compare the scores of the individual components to the full PROTAC as done in Figure LABEL:fig:protac_barplot the individual fragments (the two ligands and the linker) have to be extracted. For this, we extracted the collections containing information on anchor and warhead (the two ligands) from PROTAC-DB and tried to match an anchor and warhead to every PROTAC by matching protein target ID and find the biggest substructure match. The linker was obtained by removing the ligands’ atoms. For the 3270 PROTACs in the database we could find matching ligands for 1920 of them. This approach of extraction is necessary because PROTAC-DB does not cross-reference the PROTACs to the fragments and our extraction cannot guarantee concordance with the reported warhead and anchor in the respective publications but still is a good approximation.