Mechanistic aspects of thiol additions to Michael acceptors: Insights from computations

Computational studies have delivered valuable mechanistic insights into thiol Michael additions, which are important CS bond‐forming reactions used in biological and materials chemistry. The field has delivered a wealth of understanding about the ways in which substituents, catalysts, and the local environment influence the addition pathway. Several mechanistic scenarios are now recognized, differing with respect to the energies and timing of the bond‐forming processes. While technical challenges still exist, the field has advanced to such an extent that full‐scale simulations of the additions of Michael acceptors to protein thiol groups are now possible.


| UNLIKELY REACTION MECHANISMS
Before discussing the preferred reaction mechanisms of thiol Michael additions, we will discuss works that have explored (and largely ruled out) alternative reaction mechanisms. Unlikely mechanisms include concerted, fourcentered 1,2-additions across the acceptor C C bond. A concerted 1,2-addition of the neutral thiol has been explored by several groups. [5][6][7][8] Cronin and coworkers 7 computed transition states (TSs) for 1,2-additions (5, Figure 2) of MeSH and cysteine to 22 acceptors, while Schüürmann and coworkers 8 computed TSs for 1,2-additions to both the neutral (5) and O-protonated (6) forms of 35 acceptors. The computed barrier heights correctly predicted the experimentally measured trends in reactivities (e.g., rate constants for glutathione additions). However, the authors recognized that the TSs were not necessarily representative of the actual mechanism in the solution.
For reactions occurring under basic conditions, the O-protonated TS 6 can readily be ruled out on simple chemical grounds. The neutral TS 5, on the other hand, is ruled out on energy grounds. Cronin's calculations predicted large barriers of 50-59 kcal/mol for this type of TS with B3LYP. A recent high-accuracy CCSD(T) calculation by Rowley and coworkers 9 gave an even higher barrier of 65 kcal/mol.

| General mechanistic scenarios
It is generally accepted that a base-catalyzed thiol Michael addition involves three key processes (Figure 1b). 10 First, the thiol is deprotonated to generate the reactive thiolate anion. 2 Second, nucleophilic addition to the Michael acceptor gives an enolate intermediate (4). Third, protonation gives the neutral adduct (3). Within this general framework, computational studies have uncovered a diversity of energy landscapes and mechanistic subtleties. Two important scenarios are illustrated in Figure 3. The blue curve represents the most common case: a reaction in which the nucleophilic addition step is the rate-determining step. The pink line represents an alternative scenario in which the protonation of the enolate is the rate-determining step. In this section, we will discuss the various mechanistic processes, one step at a time.

| Identity of the base
The identity of the base that deprotonates the thiol represents an interesting question. The base catalyst could perform this function directly. Alternatively, it could act via an indirect mechanism termed "nucleophilic" initiation. 2 The F I G U R E 1 Additions of thiols to Michael acceptors F I G U R E 2 Transition states for 1,2-additions of thiols across the C C bond of Michael acceptors mechanistic differences between base catalysis and nucleophilic initiation of a thiol addition are summarized in Figure 4. In the nucleophile-initiated mechanism (orange portion of the figure), the catalyst first adds to the C C bond of the acceptor to generate a zwitterion (7). This zwitterion, rather than the base, is the species that deprotonates the thiol.
Two commonly used catalysts/initiators for thiol Michael additions are phosphines and amines. Computational studies have supported the finding that phosphines, being weakly basic, favor the nucleophile-initiated pathway. 11 For amines, however, experiments have shown and computations have supported a variety of mechanistic possibilities. [11][12][13] For example, while NMe 3 tends to act as a base catalyst, N-methyl imidazole tends to act as a nucleophilic catalyst. Qi and Wang 11 reported DFT studies to examine these preferences, focusing on additions to vinyl sulfones. Their computations with M06-2X//B3LYP in PCM dichloromethane correctly predicted the preferred activation mode for each catalyst. Importantly, they also revealed that the energy differences between the two pathways could be as small as a few kcal/mol or less. A similar finding was reported by Bowman and coworkers. 14 The above-described results suggest that both base-catalyzed and nucleophile-initiated pathways may operate, depending on conditions. Such a conclusion has been further supported by several integrated experimental and theoretical studies. These include a microkinetic modeling study performed by Reyniers and co-workers 13 (which also gave information about diffusional limitations on the reaction kinetics) and a subsequent kinetic modeling study by Bowman and co-workers 15 (which yielded further insights into the contributions of enthalpy and entropy to different mechanistic steps).
Northrop and coworkers 12 explored whether NEt 3 acted as a base or a nucleophile in the addition of MeSH to N-methylmaleimide. Computations with M06-2X//B3LYP in PCM chloroform indicated that nucleophilic addition of F I G U R E 3 Two alternative energy landscapes for base-catalyzed thiol additions F I G U R E 4 Base catalysis versus nucleophilic initiation in thiol Michael additions NEt 3 to the C C bond was surprisingly facile. Its barrier was only 2 kcal/mol higher than the barrier for the addition of MeS À from the ion pair [Et 3 N + ][MeS À ]. Nucleophilic initiation was, however, not a viable pathway for the overall reaction, because the subsequent deprotonation of the thiol by the zwitterion had a prohibitively high barrier.

| The nucleophile
In principle, either the free thiolate anion RS À or the ion pair [BH + ][RS À ] could serve as a nucleophile. In the works of Qi 11 and Northrop 12 mentioned above, the nucleophile was modeled as the ion pair. This approach is appropriate for reactions conducted in nonpolar solvents, where ion-pairing is expected to be important. Other computational studies have used the free thiolate ion as the model nucleophile. Such an approach is perhaps more appropriate for reactions taking place in an aqueous solution where ions are separated. It does not fully model the solvation of the ions, however. The molecular interactions of thiolate ions (and of thiols) with water were studied in detail by Awoonor-Williams and Rowley, 16 who performed hybrid quantum mechanical/molecular mechanical (QM/MM) molecular dynamics (MD) simulations on MeS À and MeSH in explicit water. Their simulations used a special mixed quantum/classical approach for modeling the solvent, in which the 12 water molecules nearest the solute were treated with QM while the others were treated classically. For MeSH, hydrogen bonding interactions with the solvent were found to be weak. For MeS À , stronger hydrogen bonding was observed, with a first coordination shell comprising six water molecules interacting with the sulfur. The solvation of the thiolate was found to have a strong impact on the energy landscape for its nucleophilic addition to a Michael acceptor (vide infra).

| The transition state for nucleophilic attack
Rowley and coworkers reported that numerous DFT methods were unable to locate TSs for the addition of a thiolate to a Michael acceptor in the gas phase. 9,17 For example, in the reaction of MeS À with the s-trans conformer of methyl vinyl ketone (9, Figure 5), the functionals ωB97X-D, PBE0, and M06-2X predicted tiny barriers or no barrier at all.
However, solvation has a major impact on the shape of the energy landscape. Figure 6 illustrates this point by displaying the energy profiles for the addition of MeS À to methyl vinyl ketone (now in its s-cis conformation) with and without aqueous solvent, as modeled with the CPCM implicit solvent model. 18 As would be expected for a reaction involving a small ion and a neutral molecule, the barrier in solution is substantially larger than in the gas phase. Accordingly, TSs for thiolate attack tend to be easier to locate in solution. Most of the functionals mentioned in the previous paragraph have been successfully used to calculate transition states for thiolate attack when coupled with implicit solvent models. Even B3LYP, whose deficiencies for modeling thiol Michael additions are more significant and are discussed in more detail below, has proven useful for locating thiolate addition TSs in solution.
In a sophisticated study of solvent effects, Rowley and coworkers 9 used QM/MM MD simulations to compute the reaction of MeS À with methyl vinyl ketone in explicit water. Their methodology included the wellbenchmarked ωB97X-D3BJ functional to model the QM layer. In the gas phase, the TS lay 9 kcal/mol below the separated reactants. The conversion of the non-covalently bound reactant complex into the covalently bound enolate had a small barrier (7 kcal/mol), which was entirely ascribed to the loss of entropy. In water, by contrast, there was no significant energy minimum associated with the non-covalently bound reactant complex. The nucleophilic addition TS lay 13 kcal/mol higher in energy than the separated reactants. The loss of one of the water molecules from the solvent shell of the thiolate was, in part, responsible for the higher barrier observed for the reaction in solution compared to the gas phase.
F I G U R E 5 Conformers of methyl vinyl ketone Wong et al. 19 used explicitly solvated models to study systems where solvent molecules play an integral role in the TS. They used a combination of MN15 density functional and GFN2-xTB semiempirical calculations (using their QMTSDock program) to analyze the impact of water molecules on transition states of catalyst-free, water-assisted thiol additions. It was found that the hydrogen bonding networks formed by explicit water molecules were vital in positioning the reactants in favorable positions to shuttle the proton from the thiol to the acceptor. These water molecules substantially lowered the reaction barrier heights by 10-17 kcal/mol, giving ΔG ‡ values in line with fast experimentally observed reaction rates.
Many studies have successfully employed low-cost implicit solvent models to compute thiol addition barriers, with useful accuracy. For example, Liu, Brummond, and coworkers studied thiol additions to 12 α-methylene-γ-lactams (10, Figure 7). 20 They compared the functionals M06-2X and B3LYP, the solvent models SMD and PCM, and the effect of including a quasi-harmonic correction of low-frequency vibrational modes. The best predictive accuracy (R 2 = 0.95) was obtained with M06-2X in conjunction with the SMD solvent model and quasiharmonic corrections.
When examining TS geometries, substantial variation has been observed in the forming C S bond lengths computed with different theoretical methods. To illustrate, Figure 8a shows the TSs for the addition of MeS À to methyl vinyl ketone as calculated with five different QM methods in CPCM implicit water. The forming bond length varies from 2.3 Å to 2.6 Å. The three DFT methods (B3LYP, ωB97X-D, and M06-2X) predict similar bond lengths within 0.1 Å of each other. The two other methods, HF (short C S bond) and MP2 (long C S bond), are not routinely used nowadays but have historical interest as they were used in one of the earliest theoretical studies of thiol Michael additions. 21 Of more contemporary interest is the way in which the forming bond length varies between different Michael acceptors. Substantial variation is observed. Figure 8b illustrates this point by showing the TSs for reactions of MeS À with two medicinally relevant acceptors, an acrylamide and a β-aryl-α-cyanoacrylamide. The acrylamide has a short TS bond length (2.25 Å) while the cyanoacrylamide has a long bond (2.81 Å). Interestingly, this trend goes contrary to the Hammond postulate. The acrylamide reacts more exothermically than the cyanoacrylamide, yet it has a more product-like TS with a shorter bond length.
The inclusion of a solvent model can have a significant impact on the geometry of the TS. For example, Liu and Brummond reported 20 that the forming bond length in the TS for MeS À adding to an α-methylene-γ-lactam (10, F I G U R E 6 Effect of solvent on the energy profile for the addition of MeS À to methyl vinyl ketone. The diagrams show CBS-QB3 free energies with and without B3LYP solvation energy corrections Figure 7) was 0.5 Å shorter in PCM implicit water than in the gas phase. Rowley noted a small shortening (0.1-0.2 Å) of the forming TS bond length on going from the gas phase to an explicit water model. 9 Sensitivity to the environment would also be expected to influence the TS geometries for additions to protein thiols, which take place within electronically inhomogeneous binding sites.
Another geometrical feature relevant to protein binding is the preferred orientation of the thiolate ion with respect to the C C bond. As illustrated by the comparison between Figure 8a,c, the syn conformation is typically preferred relative to anti. 18,20,22 The syn conformer avoids repulsive electrostatic interactions between the lone pairs of sulfur and oxygen while also benefitting from a stabilizing interaction between the thiolate protons and the enolate π-cloud. Further conformational flexibility also exists within the Michael acceptor itself: for example, s-cis or s-trans conformations of the α,β-unsaturated carbonyl compound, as well as side-group rotamers. When a given Michael acceptor binds to a protein thiol, the geometrical constraints within the binding site may mean that the TS is unable to adopt the conformation preferred in the solution.

| The enolate
The enolate intermediate has attracted a great deal of attention in theoretical studies. Several popular DFT methods display unexpected shortcomings in modeling the enolate. Its existence (or not) depends on the method used. Early studies focused on the enolate formed by the addition of MeS À to methyl vinyl ketone 9 ( Figure 5) in its s-trans conformation. For this enolate, Rowley showed that several popular functionals including B3LYP and PBE failed to predict an energy minimum. 17 Instead, a non-covalently bound charge-transfer complex was formed, having an electronic structure partway between [RS À …MA] and [RS•…MA• À ]. In the sense that CCSD(T) did predict a stable enolate, the B3LYP and PBE functionals appear qualitatively deficient. Rowley et al. attributed the deficiency to DFT delocalization error 9,173 and showed that it depended on a functional's treatment of exact exchange. Functionals with a large component of exact exchange globally, including M06-2X and PBE0, predicted a stable enolate, as did the range-separated functional ωB97X-D, which has a large component of exact exchange at long range. These results provide useful guidance about suitable choices of methodology for studying thiol Michael additions.
Notwithstanding the deficiencies observed for methyl vinyl ketone, a trend toward "correct" behavior has been observed for more electrophilic acceptors such as the cyanoacrylamide 11 ( Figure 9). Even B3LYP predicted the formation of a stable enolate from this acceptor and MeS À . 9 It is worth pointing out, however, that the results on enolates described in this paragraph and the one above are related to gas-phase computations. To some degree, the correct qualitative behavior can be recovered by the inclusion of a solvent, for example by using an implicit model such as CPCM or SMD. Implicit solvent is routinely used nowadays when modeling thiol Michael additions.

| Protonation of the enolate
The generic energy landscapes in Figure 3 depicted the enolate intermediate being protonated by the conjugate acid of the catalyst, BH + . However, BH + is not the only species capable of delivering a proton to the enolate. 13,15,23 Computational studies have devoted significant attention to exploring the roles of different proton donors, different protonation sites, and different timings of proton transfer relative to C S bond formation.
In an early study, Brinck, Berglund, and coworkers modeled a thiol addition that took place in the active site of a lipase. 24 Based on docking and MD simulations, they built a model system containing the substrate (acrolein), thiol (MeSH), and four key active-site residues. A histidine residue served to deprotonate the thiol and protonate the enolate. DFT calculations with B3LYP in an implicit solvent (dielectric constant 4.0) suggested that the enolate intermediate was strongly stabilized by hydrogen bonding groups in an oxyanion hole. Due to this strong stabilization, the subsequent proton transfer had such a high barrier that it was the rate-determining step of the overall addition.
There are two potential sites for protonation-carbon and oxygen. Protonation at carbon directly gives the final adduct. Protonation at oxygen gives an intermediate enol, which tautomerizes to the final keto form with assistance from the base and/or solvent. Engels and coworkers 5 explored the site-selectivity of protonation for three simple acceptors. Their computational model consisted of four components: a thiol (MeSH), acceptor (acrolein, methyl vinyl ketone, or methyl acrylate), base (NH 3 ), and proton donor (NH 4 + ), and the computations used B3LYP//BLYP with a COSMO implicit water model. For acrolein and methyl vinyl ketone, protonation took place at oxygen. For methyl acrylate, protonation at oxygen was very unfavorable; direct protonation at carbon occurred instead, in concert with the C S bond formation. Of note, the barriers for nucleophilic addition and protonation were similar in magnitude. This led to the proposal that if a proton donor or solvent molecules were not available in the required location-as might be the case inside the binding site of a protein-then the reaction may stop at the enolate or enol rather than going on to form the keto adduct. Northrop and coworkers 12 examined the importance of different proton donors. They computed the enolate derived from the addition of [Et 3 N + ][MeS À ] to N-methylmaleimide, and compared the barriers for proton donation to it by the Et 3 NH + counterion versus a molecule of MeSH. In implicit chloroform, the barrier for protonation by the thiol was only 2 kcal/mol larger than by Et 3 NH + . Given this small difference, Northrop et al. discussed how the relative importance of the two proton donors in the overall transformation would depend on their concentrations. Despite the thiol having the higher barrier, it could still be an important proton donor, especially at the beginning of a reaction when its concentration is much higher than that of the base.
Keillor and coworkers provided interesting insights into how the protonation energy landscape depends on the thiol. 25 In a combined experimental and theoretical study exploring reactions of maleimides (ωB97X-D in SMD water), aryl thiols and alkyl thiols were found to have different energy landscapes. Protonation of the enolate was high in energy for aryl thiols, while for alkyl thiols protonation and the nucleophilic attack had similar barriers. These studies modeled the proton donor as an explicit water molecule.

| The adduct
As for the enolate, computations of the neutral thiol adducts have also revealed unexpected shortcomings in certain DFT methods. For example, Krenske and Houk 18 computed the ΔH values for additions of MeSH to simple α,β-unsaturated ketones with a range of popular methods. The B3LYP, PBE0, and MPW1PW91 functionals gave large F I G U R E 9 A cyanoacrylamide errors (2-10 kcal/mol) compared to high-accuracy CBS-QB3 benchmark values. Better accuracy (<1 kcal/mol) was obtained with B2PLYP-D and SCS-MP2, as well as with M06-2X (although in this latter case only if a triple-zeta basis set was used). 4 Rowley and coworkers have also studied different methods' performance in detail. 9,17 For MeSH additions to simple Michael acceptors, they found that B3LYP gave low accuracy compared with CCSD(T) (errors 6-7 kcal/ mol), but PBE0, M06-2X, and ωB97X-D showed close agreement (errors mostly <1 kcal/mol).
Despite these recognized errors, DFT calculations have a demonstrated capacity for reproducing experimentally measured reaction energies. [27][28][29] For example, Wang and coworkers 27 reported that the M06-2X functional in conjunction with the 6-31G(d) basis set and SMD implicit solvent model of water was able to reproduce the experimental ΔG values for thiol additions to five chalcones (12, Figure 10) in phosphate-buffered saline with an accuracy of <1 kcal/mol when the experimental thiol, glutathione, was modeled as MeSH. Similarly, Krenske and Houk 18 found that M06-2X, when used in conjunction with the 6-311+G(d,p) basis set and CPCM water model, was able to reproduce the experimental ΔG values for thiol additions to five cyano-substituted Michael acceptors (e.g., 13) with an accuracy of <1 kcal/mol. 5 Thiol Michael additions-even simple ones-represent highly challenging reactions for theory. Methods that can accurately treat electron delocalization and nonbonded interactions are essential for obtaining high-quality results. The works cited in this review suggest that methods that show good performance include PBE0, PBE0-D3BJ, M06-2X, and ωB97X-D, coupled with appropriately large basis sets and an appropriate solvent model.

| Aminomethyl-assisted thiol addition
Several Michael acceptor-based drugs contain an aminomethyl substituent at the β-position, which activates the drug toward thiol addition. A longstanding question has been whether the aminomethyl group deprotonates the thiol in concert with the C S bond formation, as depicted in 14 ( Figure 11). Birkholz and coworkers 30 explored this question in the context of a series of β-aminomethyl-substituted cyanoacrylamides. Using ωB97X-D in PCM water, they compared the experimental barriers for GSH with the computed barriers for MeSH. They concluded that the trends in the experimental barriers were more consistent with a mechanism in which the aminomethyl group deprotonated the thiol ahead of the C S bond formation, not in concert with it. The protonated aminomethyl group accelerated the thiolate addition through inductive electron-withdrawing effects.

| Importance (or not) of the enol
Several groups have explored whether or not the enol intermediate is mechanistically important. Usually, the keto tautomer is significantly more stable than the enol. Tautomerization is expected to be facile in protic media. As discussed F I G U R E 1 0 Substituted Michael acceptors F I G U R E 1 1 Putative aminomethyl-assisted thiol addition above, however, Engels and coworkers 5 suggested that if there are no appropriately positioned solvent molecules available to catalyze the tautomerization, as might occur in the binding site of a protein, then the reaction might stop at the enol stage. One example of a molecule where the enol and keto tautomers of a thiol adduct were genuinely found to be similar in energy is the anti-inflammatory molecule bardoxolone methyl 28 ; the enol (15, Figure 12) was computed to be only 0.3 kcal/mol less stable than the keto tautomer (16; M06-2X in CPCM water).

| Product decomplexation
In an intriguing contribution to the field, Alabi, Clancy, and coworkers 31 examined Michael additions of fluorous thiols containing a long-chain hydroxylated substituent. Based on a combination of experimental and computational data (B97-D3), they suggested that the rate-determining step was "product decomplexation": that is, the dissociation of the ternary initiator-adduct-thiolate complex 17 ( Figure 13) to make the thiolate ion (purple) available for the next round of nucleophilic addition. The decomplexation was affected by π-stacking and hydrogen bonding interactions. Although this feature may not occur in reactions of simpler thiols, it nonetheless adds a new dimension to the mechanistic diversity of thiol Michael addition chemistry.

| Regiochemical switch to favor attack at the α-carbon
Certain substituents on a Michael acceptor can alter the regioselectivity such that addition to the α-carbon becomes favored. In 2005, Lewandowska and Chatfield 32 examined this possibility by performing experiments and computations on the additions of thiols to two β-heteroaryl-substituted acrylates. Calculations with B3LYP correctly captured the preference for α-addition in one of the two cases considered (18 vs. 19, Figure 14) but did not predict the correct selectivity for the second (a pyrimidine derivative). As was common at the time, the calculations only included solvent effects via single-point energy calculations (PCM); this may have impacted on the agreement between experiment and theory.

| "Silent" binding
In drug discovery, spectroscopic assays with simple thiols are used to determine whether a Michael acceptor is likely to bind to thiols in vivo. An interesting insight arose when exploring the naturally occurring Michael acceptor α-santonin F I G U R E 1 2 Enol and keto tautomers of thiol adducts of bardoxolone methyl F I G U R E 1 3 Ternary initiator-adduct-thiolate complex in a phosphine-mediated thiol addition (20, Figure 15). Based on an NMR spectroscopic assay using cysteamine as the thiol, α-santonin had been designated as unreactive toward thiols. 33 However, computations with ωB97X-D//M06-2X showed 34 that the barrier for thiolate attack was only 14-15 kcal/mol in implicit chloroform, which is an unremarkable barrier height for this type of process and suggested a facile addition. It was proposed that thiol addition to α-santonin actually takes place quite readily. The reason for the apparent lack of thiol binding in the spectroscopic assay was ascribed to thermodynamics, as the computed ΔG was positive. This computational finding suggested that acceptors which appear unreactive in chemical assays might still engage thiols transiently in a spectroscopically silent way. When the acceptor interacts with a protein rather than with the simple thiol, a range of additional interactions could then potentially alter the thermodynamics of binding, enabling "chemically unreactive" acceptors (i.e., those showing no evidence of binding in spectroscopic assays) to bind to proteins in a thermodynamically favorable way.

| Role of the base in controlling reversibility
In chemical synthesis, temporary thiol adduct formation provides a way to protect a C C bond. 35 Thiol additions to bicyclic and tricyclic natural products such as 22 (Figure 16) were found experimentally to have ΔG values close to zero. This allowed the position of the addition/elimination equilibrium to be tipped in a certain direction simply by manipulating the strength and concentration of the base. For example, the elimination of PhSH from 21 to give 22 was successfully achieved when stoichiometric DBU was used as the base, but not with Et 3 N. The ΔG for the elimination was +2 kcal/mol for Et 3 N but À3 kcal/mol for DBU, a difference of 5 kcal/mol. Computations with M06-2X//B3LYP successfully explained this result by showing that the extra driving force in the case of DBU arose from the reaction of DBU with the thiol. In the ion-pairing solvent dichloromethane, the computed ΔG for the acid/base equilibrium of B, PhSH, and [BH + ][PhS À ] was 4 kcal/mol more favorable for DBU than for Et 3 N, in close agreement with the experimentally observed difference of 5 kcal/mol in the ΔG values for the two eliminations.

| Asymmetric hydrogen-bonding catalysis
Asymmetric thiol Michael additions have been the subject of several computational studies. [36][37][38][39] A DFT study of an asymmetric thiol-enone addition, reported by Grayson and Houk, 37 revised the understanding of the way in which cinchona alkaloids activate substrates ( Figure 17). M06-2X calculations showed that after the thiol transferred its proton to the tertiary amine of the catalyst, the resulting ammonium cation bound to the enone oxygen atom and the hydroxyl group of the catalyst directed the thiolate attack onto a certain face of the C C bond (23). This result marked a departure from the previously assumed activation mode of the cinchona alkaloid. Previously, it had been proposed that the hydroxyl group bound to the Michael acceptor and the ammonium ion directed the attack by the thiolate (24). Subsequent computational studies have built upon this mechanistic model for thiol Michael additions catalyzed by other cinchona alkaloid catalysts. 38,39

| REACTIONS OF PROTEIN THIOLS
Given the methodological pitfalls encountered in modeling small-molecule thiol additions described above, it would be expected that the reactions of protein thiols with Michael acceptors would present a formidable task for computational modeling. Such processes combine all the above-mentioned technical difficulties as well as introducing additional challenges of their own. Nonetheless, computational investigations 40-56 of protein thiol-binding reactions have generated valuable insights into the ways in which the binding site environment affects the timing and energies of the bondforming steps.
An early study was conducted by Grazioso and coworkers, 40 who modeled the reaction of an acrylate ester with the cysteine protease falcipain-2 (FP-2). They performed ONIOM calculations on a small model of the active site and identified a mechanism in which a histidine residue acted as the catalytic base ( Figure 18a). Interestingly, the enolate underwent protonation at oxygen giving the enol in a chemically unreactive conformation. In order to generate the final (keto) product, C C bond rotation was required, so that the α-carbon could be positioned closer to the histidinium ion enabling proton transfer.
Subsequent works moved beyond small model systems to study entire proteins. For example, Lodola, Mor, and coworkers 41 reported QM/MM umbrella sampling calculations with SCC-DFTB/AMBER99SB to examine the binding of an acrylamide to the epidermal growth factor receptor (EGFR). The mechanism that emerged from their work is summarized in Figure 18b. An aspartate residue in the binding site was proposed to serve as the catalytic base. It deprotonated the cysteine residue (Cys797) and then delivered a proton to the α-carbon in concert with C S bond formation, thereby bypassing an enolate or enol intermediate. Interestingly, the computations suggested that the ratedetermining step was not any of the covalent bond-forming processes but instead was the desolvation of the intermediate thiolate ion.
Mulholland and coworkers studied the binding of the covalent anti-cancer drug ibrutinib to its biological target, Bruton's tyrosine kinase (BTK; Figure 18c). 42 BTK differs from the enzyme discussed in the above example (EGFR) in that it has no obvious residue in the binding site which could serve as the catalytic base. Furthermore, the pK a of the nucleophilic cysteine (Cys481) had been previously estimated 57 to be 10.4, which is consistent with only a very small fraction existing as the thiolate. By means of QM/MM MD umbrella sampling simulations with DFTB3/ F I G U R E 1 7 Asymmetric thiol addition AMBER, a range of possible reaction pathways were explored. The lowest-energy pathway involved direct proton transfer from Cys481 to the acrylamide oxygen. This process generated a thiolate ion which was stabilized by a nearby asparagine residue (in the same position as the Asp of EGFR mentioned above) as well as two water molecules. Subsequent C S bond formation gave an enol. Water-assisted tautomerization to the keto form was the ratedetermining step.
Developing this field further, Awoonor-Williams and Rowley 43 studied the covalent binding of a cyanoacrylamide to the same enzyme (Figure 18d). In contrast to ibrutinib, the cyanoacrylamide was designed to bind reversibly rather than irreversibly to Cys481. A combination of computational approaches including alchemical free-energy perturbation, constant-pH MD simulations, and QM/MM MD umbrella sampling was used to calculate the different steps of the reaction, including the initial non-covalent binding step, thereby mapping out the full free energy profile. A typical threestep mechanism for the covalent bond formation was found, involving proton transfer from Cys to the bulk solvent followed by nucleophilic addition of the thiolate and protonation of the enolate.
The influence of protein conformation on the reaction pathway was the topic of a study by Gleeson and coworkers. 44 They used QM/MM ONIOM calculations to study the covalent binding of an acrylate to TAK1 kinase. Four different pathways were identified, in which the thiolate, enolate, and enol intermediate assumed different levels of importance by being either present or not present as distinct intermediates. The catalytic base was an aspartate residue located relatively distant from the acceptor. Water molecules within the binding site mediated the various proton transfer steps. Importantly, Gleeson's work revealed that different conformational states of the protein favored different variants of the mechanism because each conformation positioned the catalytic base and other stabilizing groups at different distances away from the thiol and acceptor.

| CONCLUSION
Herein, we have provided an overview of mechanistic insights obtained from computational studies of thiol additions. 6 Even the simplest thiol Michael additions (e.g., MeSH plus methyl vinyl ketone) still pose difficult challenges for today's computational methods. Nonetheless, computational chemists have gained useful insights into the ways in which the energies and timing of the bond-forming and bond-breaking events vary depending on the reactants and the chemical environment, even for large systems such as protein thiols.

F I G U R E 1 8 Additions of Michael acceptors to protein thiols
Further technical developments to improve the treatments of delocalization error and nonbonded interactions should be of significant benefit in improving the ability of theory to model thiol addition kinetics and energetics. Another challenge is to bridge the gap between small (model) systems and large (experimentally studied) systems. Only a few studies have gone beyond the popular model nucleophile MeSH to tackle larger thiols, such as glutathione, used in experimental structure-reactivity studies. Explorations of the mechanistic roles of solvents, until now primarily studied in an approximate way through implicit or microsolvated models, represent a further valuable avenue for future work, to understand the significant impacts of the solvent on addition kinetics and reversibility.
Computer-based predictive tools such as machine learning 58 may provide a new way to uncover additional mechanistic details. For example, when the experimentally observed behavior of an acceptor does not match the computer prediction, detailed QM/MM MD simulations may uncover the mechanistic origins of the difference. Computational tools are already advancing the rational design of thiol acceptors in drug discovery. 59

RESEARCH RESOURCES
Computational resources for the generation of structures in Figure 8 were provided by the Australian National Computational Infrastructure under the National Computational Merit Allocation Scheme.

FUNDING INFORMATION
We thank the Australian Research Council for funding (DP180103047 to Elizabeth H. Krenske).

CONFLICT OF INTEREST
The authors have declared no conflicts of interest for this article.

DATA AVAILABILITY STATEMENT
The data that support the findings of this study are openly available in UQ eSpace, reference number UQ746e22a.

RELATED WIREs ARTICLES
Interrogating chemical mechanisms in natural products biosynthesis using quantum chemical calculations From quantum-derived principles underlying cysteine reactivity to combating the COVID-19 pandemic ENDNOTES 1 We note also that there is a parallel body of literature on radical-based mechanisms for thiol additions to C C bonds, which are not covered here. 2 An interesting computational study by considered an addition of the neutral thiol rather than a thiolate. 60 Even in the presence of a hydrogen bond acceptor capable of interacting with the SH group, this addition was uphill in energy and had no discernable TS on the potential energy surface. 3 Inclusion of an empirical dispersion correction (D3BJ) did not improve the two functionals' ability to locate the covalent enolate. 4 In an earlier study on conjugate additions, P apai and coworkers had similarly noted relatively poor accuracy for B3LYP and good accuracy for M06-2X and SCS-MP2. 26