A broadly applicable quantitative relative reactivity model for nucleophilic aromatic substitution (SNAr) using simple descriptors

We report a multivariate linear regression model able to make accurate predictions for the relative rate and regioselectivity of nucleophilic aromatic substitution (SNAr) reactions based on the electrophile structure. This model uses a diverse training/test set from experimentally-determined relative SNAr rates between benzyl alcohol and 74 unique electrophiles, including heterocycles with multiple substitution patterns. There is a robust linear relationship between the experimental SNAr free energies of activation and three molecular descriptors that can be obtained computationally: the electron affinity (EA) of the electrophile; the average molecular electrostatic potential (ESP) at the carbon undergoing substitution; and the sum of average ESP values for the ortho and para atoms relative to the reactive center. Despite using only simple descriptors calculated from ground state wavefunctions, this model demonstrates excellent correlation with previously measured SNAr reaction rates, and is able to accurately predict site selectivity for multihalogenated substrates: 91% prediction accuracy across 82 individual examples. The excellent agreement between predicted and experimental outcomes makes this easy-to-implement reactivity model a potentially powerful tool for synthetic planning.

. General competition experiment design. Two electrophiles 1 and 2 in excess but equimolar amounts were reacted with benzyl alcohol as the limiting reagent, using potassium tert-butoxide as the base. Outcomes were quantified by changes to electrophile concentration, giving a ratio of pseudo first order rate constants k1/k2.
Three typical sets of LC method parameters are shown in Tables S1 to S3 (the 3 min and 10 min ones are on the UPLC and the 6 min one is on the LCMS). The LC inlet methods were adjusted as needed for different pairings of substrates, including the mobile phase gradient, type of columns, runtime and flowrate. To ensure reliable quantitative analysis for every single competition reaction, the inlet method was set up to satisfy the following criteria: the resolution between adjacent peaks should be no less than 1.5, the tailing factor of the peaks of interest should be within 0.8-2.0, and the peak height of the peaks of interest should be above the LOQ (limit of quantification) and below the risk of overloading (usually below 1800 mAu). The UV data was collected at all wavelength range (190 -800 nm), then the spectra at either 220 nm or 254 nm were extracted to be processed for peak area integration. A representative set of LC chromatograms for a competition experiment is shown in Fig. S2.     Three of the SNAr reactions were chosen as the touchstone reactions and their absolute reaction rates were measured by kinetic analysis. The ∆G ‡ SNAr (obtained from the pseudo first order rate constant) determined for one of those touchstone reactions (2-chloropyridine and benzyl alcohol in DMSO) was used to calibrate all the relative rates determined by competition reactions, giving the absolute ∆G ‡ SNAr values for the entire array of substrates. The other two touchstone reactions were used to validate the accuracy of the rate constants determined from this competition experimentation approach, by comparing the reaction rates determined from kinetic analysis to those obtained from competition experiments.
All of the following procedures were performed inside an MBraun glovebox under an N2 atmosphere. The reactions were conducted with the same concentrations of each species as for the competition reactions (0.1 M of substrate, 0.02 M of benzyl alcohol, 0.03 M of tBuOK and 0.003 M of naphthalene or 0.02 M 1,3,5-trimethoxybenzene as internal standard in 10 mL of DMSO). All the components except for tBuOK were added into a 20 mL glass vial charged with a stir bar and mixed well by stirring, then a small amount (200 µL) of this solution was transferred into an UPLC vial followed by dilution with 800 µL of acetonitrile as the time zero (t0) point. A stopwatch was used to keep track of the reaction time immediately after the tBuOK stock solution was added. After a certain time period, a portion of the reaction mixture (500 µL) was transferred into a 2 mL UPLC vial containing 500 µL of 0.01 M hydrochloric acid in acetonitrile; the excess amount of acid was added to quench the reaction by consuming the base in the reaction solution. Reaction mixture aliquots were withdrawn at 10 time points to construct the reaction progress curves.
These reaction mixture aliquots were analyzed by UPLC and the peak areas of benzyl alcohol and the internal standard were determined. A [BnOH] versus time plot was constructed, and an exponential fit applied. Using a pseudo first-order rate law assumption, the rate constant (k) was calculated using Eq (S1). This k value was substituted into the Eyring equation to obtain the Gibbs free energy of this reaction (∆G ‡ SNAr) according to Eq (S2). The results of the three reaction rate measurement experiments are summarized in Fig. S3.
According to the pseudo first-order rate law assumption, the rate constant (k) can be determined by Eq (S1): Then the Gibbs free energy (∆G ‡ SNAr) is determined by the Eyring equation Eq (S2) (transmission coefficient κ assumed to be 1): Calculating ∆G ‡ SNAr for the entire array of Substrates All of the competition reactions were conducted under pseudo first-order conditions, with the two electrophiles in excess but equal amounts ([substrate]:[nucleophile] = 5:1). The ratio of the reaction rates can be expressed as Eq (S3): The concentrations of the two substrates are assumed to stay constant throughout the reaction; then, Eq (S1) can be simplified to Eq (S4): is equal to the percent change of the peak area of the two substrates from the reaction mixture at tstart to that at tend. By substituting the Eyring equation (Eq (S5)): into Eq (S4), the relative activation energy (∆∆G ‡ SNAr) of the two SNAr reactions from the competition can be calculated by Eq (S6): The ∆∆G ‡ SNAr values determined from the competition experiments were related to the absolute ∆G ‡ SNAr determined from kinetic analysis of the touchstone reaction between 2chloropyridine and benzyl alcohol in DMSO; therefore, all other ∆G ‡ SNAr values are given relative to 2-chloropyridine.
The SNAr rate ratio of 94 competition reactions is shown in Table S4, and the experimental ∆G ‡ SNAr for the 74 substrates determined from the 94 competition reactions in DMSO is summarized in Table S5. Substrates highlighted are those with ∆G ‡ SNAr determined by multiple competition experiments with different substrate pairings. The given ∆G ‡ SNAr values for these substrates are averages of those determined by at least two different competition experiments, and the relative standard deviation (RSD) as well as the ∆G ‡ SNAr from each individual measurement for these examples are also given.    Table S4).

Preparative Scale Synthesis of SNAr products
Five representative products from SNAr reactions were isolated, purified, and characterized using the following general procedure.
In the glovebox, a 50 mL round bottom flask containing a stir bar was charged with 3 mmol of the substrate (1.0 equiv.), 3.9 mmol of benzyl alcohol (1.3 mmol) and 15 mL DMSO. 4.5 mmol tBuOK was dissolved into 5 mL DMSO and the solution was added dropwise into the round bottom flask. The reaction mixture was stirred at room temperature overnight. On the next day, the reaction mixture was transferred outside the glovebox for workup and isolation.
The reaction mixture was poured into 30 mL water, then extracted with ethyl acetate (3×20 mL). The combined ethyl acetate phase was washed using saturated NaCl solution (30 mL) then dried over Mg2SO4, filtered and evaporated under vacuum. Then the crude product was purified by automated flash chromatography (Biotage Selekt) on silica gel using hexanes/ethyl acetate to obtain the purified product.

S4
The SNAr product S4 was prepared by the general procedure using 4-chloro-6,7dimethoxy-quinoline. During the flash chromatography purification, the targeted product could not be eluted with hexanes/ethyl acetate gradient flow, and the flash column was flushed with 100% methanol to collect the targeted product. 212.3 mg of bright yellow crystal was obtained (24% yield).

S5
The SNAr product S5 was prepared by the general procedure using 2-chloro-3methoxypyridine. 168.5 mg of colorless liquid was obtained (26% yield).

Computational Determination of Molecular Descriptors
All molecular geometry optimizations were performed using ORCA version 4.0.1.2 6 and all electrostatic potential (ESP) calculations were performed using Multiwfn version 3.7 7,8 .
Initial substrate structures were either downloaded from the ChemSpider database 9 or generated using Avogadro version 4.1. 10 Geometry optimizations were performed for the uncharged species (N-electron system) using a B3LYP/def2-TVZPD approach. The energy of the anionic species (N+1-electron system) was calculated at the optimized geometry of the uncharged species.
Using the conceptual link between LUMO energy and electron affinity (EA) related to Koopmans's theorem 11 , we calculated EA for each substrate according to: The corrections to the EA values for substrates in the main text Fig. 5 that proceed through overlap of the LUMO+1 orbitals were conducted by adding the orbital energy difference between the LUMO+1 and LUMO, as output by the DFT calculations, to the EA values.
A .gbw file that contains a binary summary of the calculation was generated automatically during the execution of ORCA, then transformed to a .wfn file by the utility program ORCA_2aim. The .wfn file contains electron density information and can be used for molecular surface analysis by the Multiwfn program. The utility program ORCA_2mkl was used to generate a .fch file, containing information on the full range of calculated MOs from the .gbw file.
The .wfn file was imported into Multiwfn for ESP calculations. The molecular ESP was calculated at a 0.004 au isosurface of electron density. By selecting Output surface properties of each atom, the surface area occupied by each atom, as well as the maximal, minimal and average ESP values at that surface were displayed. The average ESP values at the reactive center (ESP1) and sum of the ESP values at the ortho and para atoms in relation to the reactive center (ESP2) were used as the electronic descriptors in construction of the SNAr predictive model.
The electrostatic potential maps and molecular orbitals were plotted and rendered in VMD, a molecular visualization program. 12 To do so, the electron density cube file and the ESP cube file were generated in Multiwfn from the .wfn file. The two cube files can be read by VMD and the ESP map can be plotted, by projecting the ESP values at each grid point to a defined electron density isosurface. To obtain molecular orbital pictures, orbital cube files were generated in Multiwfn using .fch file as input, then plotted in VMD at a defined electron density isosurface.

Construction of the Multivariate Linear Regression Model
The selection of the molecular descriptors used to correlate ΔG ‡ SNAr with substrate structures was guided by the mechanistic features of nucleophilic aromatic substitution, 13,14 DFT calculations on the transition states, and iterative refinement of the included descriptors based on our experimental observations. As summarized in Fig. S18, there are three descriptors that lead to accurate predictions of SNAr reactivity: Electron affinity (EA) of the substrate (a measurable molecular property to approximate the LUMO energy), average electrostatic potential (ESP1) at the reactive carbon, sum of the average electrostatic potential at the ortho and para atoms in relation to the reactive carbon (ESP2). Linear regression analysis with the normalized descriptors has been performed to obtain the contribution of the individual descriptor to the reactivity of the substrate, and the results are also summarized in Fig. S19.  To validate the inclusion of these descriptors and evaluate simpler potential models, we have compared the three-descriptor model to a series of alternatives (Table S7).
First, we evaluated three univariate linear regression models constructed by each individual descriptor. The linear regression with ESP1 gives fairly good performance with an R 2 of 0.77 and an MAE of 3.0 kJ/mol, as shown in entry 2. In contrast, the univariate models with the other two descriptors only give a weak linear relationship (EA, entry 4) and almost no linear relationship (ESP2, entry 3).

Ar-F (9)
Ar-Cl (51) Ar-Br (14) ESP1 is the descriptor that contributes the most to predicting SNAr reactivity among the three, as indicated by the performance of the univariate models. Then we evaluated two bivariate linear regression models as shown in entries 5 and 6. Compared to the univariate model in entry 2, the inclusion of one more descriptor, either EA or ESP2 improves the prediction accuracy, but still in a slightly lower level compared to the prediction accuracy of the multivariate model with all the three descriptors included.
The inclusion of ESP2 (sum of the ESP at the ortho and para atoms) in the model was explored based on our DFT transition state calculations and the experimental observations. From our iterative analysis of the linear regression model, we discovered that the two ESP values -ESP at the ortho atoms and the ESP at the para atomhave approximately equal contributions when linear regression fitting is done with ESPortho and ESPpara as separate descriptors. This alternate, four-descriptor model is shown in entry 7. An R 2 of 0.92 and a MAE of 1.9 kJ/mol were obtained from the four-descriptor linear regression model, which are almost identical to those obtained from the original three-variable regression model (entry 1). The two models give the same level of prediction accuracy, and thus the ESPortho and ESPpara are treated as one summed value for simplification.
The model with steric effect as an additional descriptor (steric A value 15 ) has also been evaluated. The statistic measures in the entry 8 show that the prediction accuracy is almost identical compared to the initial model shown in Fig. S19. Linear regression analysis with normalized descriptors shows that the contribution of the steric A value is only 2%, which is much lower than any other descriptor. As a result, the steric descriptor is not included because of its insignificant impact to the reactivity of the substrates in our training dataset.
Overall, the alternatives presented in Table S7 demonstrate that all of the three descriptors are necessary and sufficient to predict ΔG ‡ SNAr for this diverse set of substrates.

Cross Validation and Out-of-Sample Prediction
To further evaluate the linear regression model, we performed cross-validation by doing five random 60/40 training/test data splits (Figs S20-S24). Excellent linear correlation was achieved between the observed and predicted ΔG ‡ SNAr, as indicated by the range of R 2 from 0.89 to 0.93 for training set, Q 2 from 0.86 to 0.93, and MAE from 1.63 kJ/mol to 2.30 kJ/mol for test set. The good agreement between the observed and predicted ΔG ‡ SNAr obtained from this random split cross-validation has indicated that our multivariate linear regression model is appropriately fitted, with no overfitting issues.
We also calculated the 95% prediction intervals of the predicted ΔG ‡ SNAr for the test set for one of the 60/40 random split model (Fig. S20). The narrow prediction intervals (±5.1 kJ/mol to ±5.5 kJ/mol) indicate that the model makes reliable prediction of the SNAr rates.         The model performance was further evaluated by an out-of-sample prediction. To test if the model can give reliable predictions for molecules with a variety of structural features, we split the data set into a training set containing only the chloropyridine substrates, and a test set containing all other substrates. The model has achieved an excellent performance with a R 2 of 0.92 for the training set, and a MAE of 1.81 kJ/mol and a Q 2 of 0.93 for the test set (Fig. S25).   We use the electron affinity of a substrate as an approximation to its LUMO energy. Electron affinity is a measurable physical property, it is defined as the amount of energy released after an electron is attached to a neutral molecule. Electron affinity is calculated as the energy difference between a neutral molecule and its anionic radical by DFT. The orbital energy calculation requires less computational effort since only the DFT calculation on the neutral molecule is needed. While LUMO is a virtual orbital, and the calculated virtual orbital energy can introduce uncertainties since it is sensitive to the selection of the basis set.
We compared the LUMO energies to the electron affinity calculated at the same level of theory and basis set (B3LYP/def2-TVZPD), obtaining a linear correlation as shown in Fig.  S26. The strong correlation (R 2 =0.94) indicates that the LUMO energies may be used as an alternative to the electron affinity to reduce computational cost. Then we constructed a multivariate linear model using LUMO energies and ESP values as shown in Fig. S27, and the model performance (R 2 =0.92 and MAE=1.91 kJ/mol) is equivalent to the model using electron affinity (Fig. S19).  LUMO energy at GFN1-xTB level by Entos Envision To assess whether a rapid method at a lower level of theory could be suitable for generating LUMO energies that give accurate predictions, we evaluated Entos Envision 16 , an interactive web-based platform for molecular simulation and visualization developed by Entos, Inc. It requires the molecular structure as input to run computational calculations using GFN1-xTB, 17 a semi-empirical tight-binding method, and outputs a series of properties including molecular orbital energies. Calculation at this theory level only takes a few seconds for the size of a typical substrate in SNAr reaction, much faster than hybrid functional DFT calculations, which may run for hours. We compared the LUMO energies computed in Entos to electron affinity at B3LYP/dev2-TVZPD level, obtaining a linear correlation as shown in Fig. S28. This good correlation, with an R 2 = 0.88, indicates that these semi-empirical calculations give reasonable LUMO energy values. We used the LUMO energies obtained from Entos in our multivariate regression model, and compared its performance with the model built using DFTcalculated electron affinity (Fig. S19). The two models are almost identical in their performance, in terms of R 2 and MAE values (Fig. S29).     Electron Affinity (kJ/mol)

Transition State Calculations of SNAr by DFT
All transition state calculations were performed in ORCA 4.0.1.2. Geometry optimization and frequency calculations for starting materials, intermediates, transition states and products were performed at B3LYP/def2-SVPD level with D3BJ dispersion for all atoms and implicit solvation model CPCM when solvent effect is included. Additional single point calculations were performed at B3LYP/def2-TVZPD level with D3BJ dispersion and CPCM when needed.
Transition state geometries were sought through relaxed scans along a simplified reaction coordinate (the breaking or forming bond length), at a B3LYP/def2-SVPD/CPCM (when needed) level with D3BJ dispersion. The geometry at the maximum energy point of the reaction coordinate scan was used as input geometry for continuing the transition state search via saddle-point optimization.
The optimization to a saddle point was performed using a B3LYP/def2-SVPD/CPCM (when needed) approach with D3BJ dispersion. The transition state structure was confirmed by vibrational analysis. The vibrational mode corresponding to the single imaginary frequency was visualized to ensure that it lied along the reaction coordinate bond stretch.
The free energies of the starting materials, intermediates and transition states were calculated as the sum of the electronic energies (single point energy at B3LYP/def2-TVZPD level) and the Gibbs free energy correction (obtained from frequency calculations at B3LYP/def2-SVPD level). Gibbs free energies of activation (ΔG ‡ SNAr) for the reactions were calculated as the difference in free energies between the transition state and the starting materials.
These labor-intensive transition state calculations were pursued for 6 di-halogenated heterocycles at both reactive sites except for 2,4-dichloro-1,3,5-triazine, for which the two sites are identical. In implicit solvation model CPCM(DMSO), concerted pathways were found at all reactive centers, and no sigma-complex intermediates could be located; the starting materials transformed into the products via only one intrinsic reaction coordinate (the bond length between the reactive carbon at the substrate and the negatively charged oxygen at deprotonated benzyl alcohol). Transition state structures and energy diagrams are shown in Figs. S30 -S39. All molecular images were generated in Vesta. 18 In gas phase, stepwise mechanisms were found for the carbon-fluorine sites at 2-chloro-3-fluoro-pyridine and 2-chloro-4-fluoro-pyridine (Figs. S38 -S39). The stable sigma-complex was located from geometry optimization and confirmed by frequency calculations, that no negative frequencies were found. The first step, addition of the nucleophile, is an exothermic reaction that involves a very early transition state. In these cases, the transition state could not be found as no maximum point was accessible from the potential energy surface scan. This computational challenge has been reported by Jacobsen and co-workers' in their SNAr mechanistic study, where transition state searches failed in similar exothermic steps. 13 Transition states for the second step (the departure of the fluorine) were located and the structures were optimized and confirmed by vibrational analysis for both of the carbon-fluorine sites.
HOMO diagrams for the transition state were calculated at both the major and minor reactive sites for the following 3 substrates: 2,4-dichloropyridine (Fig. 6, main text), 2-chloro-4-bromopyridine (Fig. S32) and 2,4-dichloropyrimidine (Fig. S34). One motivation to conduct transition state calculations is to look for the mechanistic factors that determine the intramolecular selectivity, and our orbital analysis shows that the frontier orbital interaction is certainly among those important factors. Transition state frontier orbital analysis shows that different sites of a substrate use different antibonding orbitals to interact with the nucleophile, where those orbitals are accessible through thermal fluctuations. In these cases, the major site undergoes SNAr corresponding to the LUMO orbital, while the minor site corresponds to LUMO+1.

Fig. S30.
Calculated reaction coordinates for SNAr of 2,4-dichloropyridine and the anion of benzyl alcohol in DMSO. The blue pathway is the SNAr coordinate at C4, the major reactive site; and the red pathway is the SNAr coordinate at C2, the minor reactive site.   Frontier orbitals of 2-chloro-4-bromopyridine and the anion of benzyl alcohol involved in SNAr, and HOMO orbitals of the transition states at C2 and C4. The major reactive site C4 undergoes SNAr with the LUMO orbital, as shown as the blue reaction pathway in Fig. S31. The minor reactive site C2 undergoes SNAr with the LUMO+1 orbital, as shown as the red reaction pathway in Fig. S31.

Fig. S33.
Calculated reaction coordinates for SNAr of 2,4-dichloro-pyrimidine and the anion of benzyl alcohol in DMSO. The blue pathway is the SNAr coordinate at C4, the major reactive site; and the red pathway is the SNAr coordinate at C2, the minor reactive site.  Experimental details for assessing outliers from external case studies (Fig. 7C) To assess the model's applicability on new SNAr reactions beyond our training dataset, we conducted a series of external validation on both SNAr rate correlation (main text, Fig. 7) and site-selectivity prediction (main text, Fig. 8-12). We also performed experimental checks on the one significant outlier that we observed from the SNAr rate correlation test (Substrate 6, Fig. 7C).
The substrate 6 from Fig. 7C (1-bromo-4-fluoro-2-nitrobenzene) has two potential SNAr sites (Ar-F and Ar-Br) and is expected to form multiple products from reaction with piperidine.
The substrate was purchased from a commercial supplier and used as is. The reaction was conducted using the same reaction conditions from literature 22 (which is also shown in Fig.  S40A) for two days, then the reaction solution was analyzed by LCMS. The molecular ion peaks from the LCMS analysis gives experimental evidence that the reaction of substrate 6C and piperidine leads to two SNAr products. The product ratio of 1.5:1 was obtained from the peak area ratio from the UV spectrum at 254 nm, with the Ar-Br site slightly favoured.
Excellent linear correlation (R 2 = 0.90) was obtained for the dataset in Fig. 7C with all the substrates included, and exclusion of the outlier 6C gives rise to a stronger correlation, with an R 2 of 0.94. Fig. S40. A). Reaction condition of substrate 6 in Fig.7C reacting with piperidine. B). UV chromatogram at 254 nm extracted from LCMS analysis of the reaction solution, and the m/z of the two SNAr products.