Predicting higher activity catalysts via regression modeling in asymmetric catalysis: rebuilding the Denmark 2019 machine-learning workflow from scratch

Introduction

For those new to machine learning in chemistry, the high-level summary of how to build a regression model to predict molecular properties is: (1) assemble training data and compute descriptors, (2) train machine learning models using your descriptors and some target output variable like enantiomeric excess (ee), and (3) make predictions on unseen molecules or experiments to guide synthesis and testing.

Those excited about this space should be aware of Professor Scott Denmark's group's breakthrough 2019 publication demonstrating that machine-learning models can predict the ee of unseen asymmetric transformations, even when the test data falls outside the range of the training data. They demonstrated this on an enantioselective N,S-acetal formation using imine electrophiles, thiol nucleophiles, and chiral BINOL phosphoric acid catalysts. The experimental dataset contains a library of 43 chiral BINOL phosphoric acid catalysts, five imines, five thiols, and a total of 1,075 experimental ee measurements (Figure 1).

At a high level, Denmark and co-workers featurized catalysts using the grid-based average steric occupancy (ASO) descriptor, which gives a measure of which points on a 3D grid a catalyst's conformer ensemble occupies. ASO is a set of descriptors, or descriptor vector, where each grid point is a unique descriptor across catalysts, so the representation is high-dimensional. Imines and thiols were featurized using the same ASO approach on their respective conformer ensembles, with each substrate class aligned to its own common core. Electronic effects were captured via a substituent-based descriptor in which the variable group of interest was appended to a tetramethylammonium cation and the electrostatic potential surface was computed with DFT at B3LYP/6-31G*, with the minimum and maximum ESP values serving as electronic descriptors. For thiols, two additional electronic descriptors were included: the sulfur natural bond orbital (NBO) charge, also called NPA charge, on sulfur, and molecular orbital energies. These descriptor sets were used to train different types of ML models, but the most productive were ultimately support vector machines and feed-forward neural networks.

Part 1. Computing catalyst, imine, and thiol descriptors from scratch

To do our own modeling on this dataset, we started by calculating DFT-level descriptors. First, in CDXML files, we drew out molecules of interest and labelled atoms of conserved substructures that we wanted consistent atom- and bond-level descriptors for. We uploaded three separate CDXML files to Chemetrian and used them to make three separate libraries: a library of 43 labeled BINOL phosphoric acids, a library of 5 labeled imines, and a library of 5 labeled thiols.² Each library was then submitted to the DFT-level descriptor calculation pipeline, which includes 2D→3D conversion, conformer searching,³ geometry optimizations on GPU using Meta's UMA-s-1.2 neural network potential,⁴ and GPU-accelerated DFT at the r²SCAN-D4/def2-TZVP level of theory⁵ (Figure 2).

This provided three descriptor sets, one for each reaction component. Total descriptor pipeline calculation times were 20 minutes for five imines, 13 minutes for five thiols, and only 4 hours of wall-clock time for 43 catalysts, with some having > 180 atoms, because we parallelized the library across five GPUs. We then downloaded the N,S-acetal formation ee data from the Denmark group's GitHub, and used Claude to efficiently create a CSV of combined catalyst, imine, and thiol descriptors as well as the experimental ee labels for the various combinations of reaction components tested in the paper.

Part 2. Building regression models on the full dataset

With the descriptors in hand, we were able to start building regression models. We started by building the same model in Fig. 6 of the Denmark 2019 paper: a random train/test split repeated 10× to obtain average metrics, with the main success metric being mean absolute error (MAE) of ΔΔG‡ predictions (ΔΔG‡ represents the energy difference between the two diastereomeric transition states that yield different enantiomers of product, and is used instead of ee as the target). Following the paper, we trained a support vector machine (SVM) regressor with a second-order polynomial kernel and got essentially the same result Denmark and co-workers reported: a test MAE of 0.150 ± 0.010 kcal/mol across the 10 splits, against their 0.152 ± 0.005 (Figure 3, left).

This was exciting, considering we used a totally different featurization approach than the Denmark lab. Despite different conformer search, geometry optimization, and DFT engines, our pipeline produced molecular descriptors that support the same level of predictive accuracy as the original work, and these features were accessible through a simple user interface in less than five hours from start to finish. An important note: we did no feature selection. We simply inputted the 335 features per experiment into the models. Feature selection may have improved our results.

From there we tested a few other models. Plain linear regression, ridge, and Lasso were all clearly worse, underfitting relative to SVM. A random forest model came in at 0.136 ± 0.011 kcal/mol and XGBoost at 0.135 ± 0.010 kcal/mol (Figure 3, middle and right), both essentially matching the SVM on the test set, though with more aggressive fitting on the training data. Notably, all models tend to underpredict the ΔΔG values for the top-performing catalysts, which was also observed by the Denmark group.

Part 3. Predicting > 80% ee experiments by training only on < 80% ee data

Now, we will skip to the model described in Fig. 8 of the paper. The authors used < 80% ee data to train a neural network to make predictions on experiments that yielded > 80% ee. This is a notable challenge and requires model extrapolation instead of simply interpolation. This also represents one of the "holy grails" of regression modeling in catalysis, which is the idea of using ML to discover higher-activity catalysts. The authors found that a feed-forward neural network was able to extrapolate to the test dataset and provided 0.33 MAE, but still underpredicted ΔΔG values. We emulated this prediction exercise but opted for classical machine learning models instead of neural networks due to the substantially faster rate that they can be tested and optimized. SVM regression with a second-degree polynomial kernel (test MAE = 0.43) came close to matching the performance of Denmark's neural network (Figure 4A). It also had the same issue of underpredicting ΔΔG at the high values of the dataset.

Denmark and co-workers emphasized that the real test of these models isn't accuracy at the high-ee extreme, where both their model and ours systematically underpredict, but whether the rankings would guide a chemist toward the best catalysts. By this logic, both Denmark's neural network and our SVM were successful (Figure 4B). For our SVM, all eight catalysts that experimentally yielded > 80% ee were predicted above 80% ee, meaning the model successfully extrapolated outside the training distribution and correctly identified high-performing catalysts. More importantly, the top three catalysts observed to experimentally provide > 95% ee are also the top three predicted by the model.

Denmark's neural network is more accurate in absolute terms and avoids some ranking errors that our model makes. But considering that our SVM was built on DFT descriptors computed end-to-end through a user interface, that no coding was required, that no feature selection was performed, and that the model was only tuned with a quick grid search over hyperparameters, the fact that we correctly identified the top three catalysts and all eight > 80% ee molecules represents exciting and actionable results. A chemist running this workflow in practice would synthesize the same three catalysts our model and Denmark's both point to, and would arrive at > 95% ee in the first round of screening. That is the optimal outcome for an ML-driven catalyst design campaign.

Part 4. Feature analysis to gain chemical insights on the system

Machine-learning models can be interrogated to gain insights on the fundamentals driving the chemistry of a system. Specifically, understanding which descriptors are most responsible for the model's predictive accuracy can inform molecule design and yield knowledge that translates to future systems outside the scope of the present study. To gain such insights, we calculated feature permutation importances for the model trained on < 80% ee data. Using this method, features are permuted and the effect on the model is quantified. The greater the change in the model's accuracy, the more important a feature is. Our DFT-calculated descriptors are well-suited for this analysis because they are interpretable to chemists (Figure 5).

Analysis of the top ten features by permutation importance highlights that geometric, steric, and electronic catalyst features are important for performance. The highest importance is the H–I bond length, the distance between the naphthyl C3 (atom H) and the first atom of the 3,3′-substituent (atom I). Top catalysts cluster tightly at ~1.49 Å, which is standard for a Csp²–Csp² bond between two aryl rings. Lowest-performing catalysts sit at 1.7–1.9 Å, the range characteristic of C–Si, sp³ C–C (–CH₂–Ar spacers), or C–halogen bonds in the dataset. The descriptor thus acts as a readout of linker identity, and the associated charge features on atoms H and I (NPA ≈ −0.07 vs −0.15, Mulliken on the ortho-C ≈ +0.1 vs −0.2) reinforce the same distinction, since an aryl–aryl link leaves the ipso carbon mildly electron-poor whereas silyl or alkyl groups push electron density onto it.

Sterically, four features describe how bulk is distributed across the chiral pocket: lower buried volume around phosphorus, higher buried volume on the substituent's ipso carbon, and larger Sterimol B1 along F–G and Sterimol L along G–H.⁶ The best catalysts keep the Brønsted acidic site unhindered at short range but more hindered further away. Electronically, the active site is tuned to be moderately, not maximally, acidic: the Fukui dual descriptor on the acidic proton (atom C) is essentially zero in top performers and the Mulliken charge on phosphorus is slightly less positive for better catalysts than worse catalysts. In summary, permutation importance analysis reveals three design principles: the 3,3′-substituent must be directly aryl-linked, steric bulk should be projected away from the Brønsted acidic site rather than crowded onto it, and the catalyst's acidity should be moderate rather than maximal.

Conclusion

We have outlined a full descriptor calculation and machine-learning pipeline accessible via the Chemetrian platform through a simple user interface. Our goal is to accelerate the pace at which frontier machine-learning approaches for molecular design and optimization can be implemented, and to make them accessible to chemists without machine-learning or software-engineering backgrounds. The case study presented here, from raw catalyst structures through DFT-calculated descriptors, a trained model, and the chemically interpretable insights extracted from it, was executed end-to-end in a few hours. The same workflow would typically require weeks of effort from a computational chemist comfortable with conformer searching, electronic structure packages, descriptor extraction codes, and machine-learning frameworks. By collapsing that timeline and removing the software-engineering barrier, Chemetrian puts predictive modeling on the same timescale as the experimental campaigns it is meant to inform, and brings it within reach of chemists able to act on its predictions.