Mechanism of covalent binding of ibrutinib to Bruton's tyrosine kinase revealed by QM/MM calculations

Ibrutinib is the first covalent inhibitor of Bruton's tyrosine kinase (BTK) to be used in the treatment of B-cell cancers. Understanding the mechanism of covalent inhibition will aid in the design of safer and more selective covalent inhibitors that target BTK. The mechanism of covalent inhibition in BTK has been uncertain because there is no appropriate residue nearby that can act as a base to deprotonate the cysteine thiol prior to covalent bond formation. We investigate several mechanisms of covalent modification of C481 in BTK by ibrutinib using combined quantum mechanics/molecular mechanics (QM/MM) molecular dynamics reaction simulations. The lowest energy pathway involves direct proton transfer from C481 to the acrylamide warhead in ibrutinib, followed by covalent bond formation to form an enol intermediate. There is a subsequent rate-limiting keto–enol tautomerisation step (ΔG‡ = 10.5 kcal mol−1) to reach the inactivated BTK/ibrutinib complex. Our results represent the first mechanistic study of BTK inactivation by ibrutinib to consider multiple mechanistic pathways. These findings should aid in the design of covalent drugs that target BTK and other similar targets.


S1.1 Method validation
We performed 1D potential energy scans of a model system comprising methyl thiolate and an acrylamide warhead attached to a piperidine linker as an ibrutinib mimic using Gaussian16. 1 Geometry optimisations were performed along the covalent C-S bond formation reaction coordinate, that was varied between 1.6 and 3.0 Å in steps of 0.01 Å. The results suggested that B3LYP does not predict a stable minimum for the carbanion intermediate ( Figure S1). This is in agreement with previous studies that show significant failings for B3LYP and that Range-Separated DFT Functionals are Necessary to Model Thio-Michael Additions. 2,3 In addition, DFTB2 and DFTB3//DFTB2 incorrectly predict a carbanion intermediate at S-C distances of 2.45 Å ( Figure S2). DFTB2 geometries are used because there is no DFTB3 implementation in Gaussian16. The energy jump observed in the potential energy surface scans for AM1, PM3 and PM6 at large S-C distances is the result of isomerisation of the olefin group in acrylamide. Figure S1. Potential energy profiles of the nucleophilic attack step in thiol-Michael addition between methyl thiolate and acrylamide at the DFT/aug-cc-pVTZ level. All methods tested show a stable enolate intermediate, apart from B3LYP/aug-cc-pVTZ in accordance with previous studies. 2,3 Figure S2. Potential energy profiles of the nucleophilic attack step in thiol-Michael addition between methyl thiolate and acrylamide computed using semi-empirical quantum chemistry methods. The potential energy profile at the MP2 level of theory is shown for comparison. The energy jump observed in the potential energy surface scans for AM1, PM3 and PM6 at large S-C distances is the result of isomerisation of the olefin group in acrylamide.
To test the accuracy of the geometries predicted for C-S bond formation predicted by PM6 and DFTB3//DFTB2 along the reaction path, we performed single point energy calculations on the structures along the DFTB3//DFTB2 and PM6 reaction pathways with the wB97X-D density functional in combination with the aug-cc-pVTZ basis set (denoted wB97X-D//DFTB2 and wB97X-D//PM6), Figure S3. The potential energy scan for S-C formation at the wB97X-D//DFTB2 level is very similar to the potential energy scan at the wB97X-D/augcc-pVTZ level, indicating that the geometries predicted by DFTB2 are reliable and are suitable for modelling thiol addition reactions. This was not the case at the wB97X-D//PM6 level, where the PM6 geometries led to approximately 10-15 kcal mol -1 higher energies across the surface. Next, we carried out a transition state (TS) search for the reaction shown in Figure   S6(B) with a variety of density functional and semi-empirical QM methods. We confirmed the presence of a TS by checking for a single imaginary frequency corresponding to the reaction coordinate and performed intrinsic reaction coordinate (IRC) calculations to investigate the pathways predicted by different methods for thiol addition by plotting the value of the reaction coordinate for S-C bond formation and proton transfer from methyl thiol to the carbonyl oxygen of the truncated ibrutinib model ( Figure S4). All of the density functional methods and DFTB2 predict a similar reaction pathway along the nucleophilic attack and proton transfer reaction coordinates. AM1 was in reasonably close agreement, whereas PM3, PM6 and PM7 were in poor agreement with the density functional methods and predicted significantly different reaction pathways in which the C-S bond formation occurs prior to proton transfer.
The corresponding energy plots for each IRC pathway ( Figure S5) suggest that DFTB2 underestimates the barrier to S-C formation and proton transfer. These results confirm the suitability of using the DFTB3 method for modelling thiol addition reactivity of this type. Figure S4. Plot to show the nucleophilic attack and proton transfer reaction coordinates along the reaction path predicted by IRC calculations with several different QM methods for mechanism 3 (the lowest energy pathway), where the thiol proton transfers to the carbonyl oxygen of acrylamide to form an enol. The location of transition states for each method are shown as circles on each pathway. Figure S5. Energy plot of the IRC calculations suggest that DFTB2 underestimates the barrier to S-C formation and proton transfer in a model BTK/ibrutinib complex ( Figure S6(B)).
The IRC pathway at the DFTB3 level ( Figure S4) was in better agreement with higher level methods than that for S-C bond formation alone ( Figure S2). Examination of the structures and important bond lengths of the TSs for DFTB3, MP2 and wB97X-D ( Figure S6 and Table S1) shows a high degree of structural similarity between the TSs. The S-C and O-H distances predicted by DFTB3 are closer to the values predicted by MP2 than the values predicted by wB97X-D (Table S1). In addition, the IRC for solvent-assisted tautomerisation confirms that DFTB3 gives energetics of this reaction step similar to wB97X-D, showing that DFTB3 provides a good description of the reaction ( Figure S7).  Table S1. Tabulated bond lengths in the optimised TS corresponding to the model BTK/ibrutinib system ( Figure S6(B)). . Potential energy profiles from IRC calculation of the rate-limiting solvent-assisted tautomerisation step on a model BTK/ibrutinib system ( Figure S6(B)) in the gas phase at the DFTB2 and level and wB97X-D/6-31G(d) level. The DFTB2 energies are in good agreement with wB97X-D/6-31G(d) for this solvent-assisted tautomerisation reaction.

S1.2 Model setup
Snapshots for the QM/MM reaction simulations were selected from 500 ns of classical molecular dynamics simulations using the AMBER ff14SB force field. Structural coordinates of the covalently bound BTK/ibrutinib complex were taken from the crystal structure with PDB ID 5P9J. 4 Missing residues were included by using the Modeller program. 5 Hydrogens atoms were added using the protein preparation wizard in Maestro, 6 where PROPKA3.1 7,8 was used to assign protonation states of titratable amino acid side chains and optimise the hydrogen bonding network. For the covalently bound C481-ibrutinib complex, RESP charges were generated using the REDServer 9 (RESP ESP charge Derive Server) and any missing FF parameters were generated using Antechamber, distributed the AMBER 2018 package. 10 The The free energy surfaces were generated using the weighted histogram analysis method (WHAM), using code distributed by the Grossfield lab. 13,14 Scheme S1. Schematic of the four mechanistic pathways explored for the covalent inhibition of C481 in BTK by ibrutinib.

S1.4 Umbrella sampling convergence
To check convergence of the umbrella sampling pathways, histograms of the number of counts of each sampled reaction coordinate value in each reaction coordinate window were produced. 15,16 An example of the histograms along the proton transfer and S-C reaction coordinates from mechanism 3 (the lowest energy pathway) are shown below. These show good overlap between neighbouring umbrella sampling windows and uniform heights of each histogram. Figure S9. Umbrella sampling histograms along the proton transfer and nucleophilic attack (S-C formation) reaction coordinates for mechanism 3 after 25 ps of sampling in combination with a biasing potential of 200 kcal mol -1 Å 2 in each umbrella sampling window.
We also checked that the barrier heights had converged by performing WHAM analysis on the complete dataset and 80% of the dataset). The free energies were considered converged if the change in the barrier height (DG ‡ ) was less than 0.1 kcal mol -1 . The free energy change in each reaction coordinate window when comparing 100% and 80% of the sampling data across the entire free energy surface for the solvent-assisted tautomerisation step is less than 0.4 kcal mol -1 , and the free energy change in the barrier height is 0.02 kcal mol -1 ( Figure S10).
These tests show that the QM/MM umbrella sampling results are well converged. Figure S10. Heat map to show the difference in free energy in each reaction coordinate window for the solvent-assisted tautomerisation free energy surface (mechanisms 3 and 4) when 100% and 80% of sampling data are compared at the DFTB3/MM level of theory. S1.5 Alternative mechanistic pathways S1.5.1 Mechanism 1 Figure S11. Approximate transition state, taken from the final frame of the highest energy reaction coordinate window (RC1 = 2.6 Å, RC2 = 0.6 Å) along the reaction pathway of Mechanism 1. (B) Representative snapshot of the covalently bound BTK/ibrutinib keto product (EP) from mechanism 1. The geometry around the newly formed S-C bond differs from the crystal geometry (PDB 5P9J, 4 shown in transparent green), and results in a high energy product state. Figure S12. Free energy surface produced at the DFTB3 level, with 30 ps of sampling per window along the minimum energy path for the direct addition of ibrutinib to C481 in BTK. The high reaction barrier of 47.7 kcal mol -1 is consistent with the reaction being thermally forbidden.

S5.1.2 Mechanism 2
For mechanism 2, we had trouble generating a productive reaction pathway and were only able to successfully generate an approximate 2D free energy pathway along the diagonal between each reaction coordinate. Reaction coordinate (RC) 1 was defined by proton transfer from the a-carbon of the inhibitor to the oxygen of a water molecule and S-C bond formation.
RC2 was defined by an additional proton transfer from the same water molecule to the C481 sulfur atom. Exploration of the remainder of the free energy surface resulted in additional proton transfers, even when additional restraints were added. As a result, the free energy barrier for mechanism 2 is approximate, but the numerous failed attempts to obtain a productive pathway, and the high reaction barrier suggest this is not a feasible reaction mechanism for covalent modification of C481 by ibrutinib. Figure S13. (A) Approximate free energy profile for mechanism 2, where solvent assisted thiol addition occurs to result in the keto product. The approximate free energy barrier is 39.0 kcal mol -1 and 2 ps of sampling was carried out in reach umbrella sampling window. (B) Approximate transition state for mechanism 2. The structure is the final frame from the highest energy reaction coordinate window along the reaction path.

S5.1.3 Mechanism 3
In order to confirm that the formation of the C481-S -/N484-NH2 interaction is dependent on the orientation of the amide side chain of N484, rather than the QM method used, free energy surfaces were produced at the DFTB3 and wB97X-D/6-31G(d) levels of theory ( Figure S14). The starting points for both simulations were a single snapshot in which N484 was oriented away from C481. The reaction was modelled in the forwards direction (from ES to EI2). Comparison of the free energy surfaces produced at the DFTB3 and wB97X-D/6-31G(d) levels of theory reveals that no stable EI1 intermediate is predicted by either method as a result of the absence of a C481-S -/N484-NH2 interaction. The barrier to enol formation is 23.9 and 29.2 kcal mol -1 for DFTB and wB97X-D/6-31G(d) respectively, indicating that DFTB underestimates the barrier for S-C formation by approximately 5 kcal mol -1 . The reaction barrier is much higher when the reaction is modelled as a concerted process due to unfavorable geometry of the TS, and the absence of the stabilising interaction from N484. Figure S17. Free energy barrier heights of the four mechanistic pathways considered for the reaction between ibrutinib and C481 in BTK. Mechanism 3 is the lowest energy pathway with a rate-limiting step that corresponds to solvent-assisted keto-enol tautomerisation (DG ‡ = 10.5 kcal mol -1 ).