Profiling MAP kinase cysteines for targeted covalent inhibitor design

QM/MM Simulations of Afatinib-EGFR Addition: The Role of β-Dimethylaminomethyl Substitution

Acrylamides are the most commonly used warheads of targeted covalent inhibitors (TCIs) directed at cysteines; however, the reaction mechanisms of acrylamides in proteins remain controversial, particularly for those involving protonated or unreactive cysteines. Using the combined semiempirical quantum mechanics (QM)/molecular mechanics (MM) free energy simulations, we investigated the reaction between afatinib, the first TCI drug for cancer treatment, and Cys797 in the EGFR kinase. Afatinib contains a β-dimethylaminomethyl (β-DMAM) substitution which has been shown to enhance the intrinsic reactivity and potency against EGFR for related inhibitors. Two hypothesized reaction mechanisms were tested. Our data suggest that Cys797 becomes deprotonated in the presence of afatinib, and the reaction proceeds via a classical Michael addition mechanism, with Asp800 stabilizing the ion-pair reactant state β-DMAM+/C797- and the transition state of the nucleophilic attack. Our work elucidates an important structure-activity relationship of acrylamides in proteins.

Evaluation of a Covalent Library of Diverse Warheads (CovLib) Binding to JNK3, USP7, or p53

Purpose

Over the last few years, covalent fragment-based drug discovery has gained significant importance. Thus, striving for more warhead diversity, we conceived a library consisting of 20 covalently reacting compounds. Our covalent fragment library (CovLib) contains four different warhead classes, including five α-cyanoacacrylamides/acrylates (CA), three epoxides (EO), four vinyl sulfones (VS), and eight electron-deficient heteroarenes with a leaving group (SNAr/SN).

Methods

After predicting the theoretical solubility of the fragments by LogP and LogS during the selection process, we determined their experimental solubility using a turbidimetric solubility assay. The reactivities of the different compounds were measured in a high-throughput 5,5'-dithiobis-(2-nitrobenzoic acid) DTNB assay, followed by a (glutathione) GSH stability assay. We employed the CovLib in a (differential scanning fluorimetry) DSF-based screening against different targets: c-Jun N-terminal kinase 3 (JNK3), ubiquitin-specific protease 7 (USP7), and the tumor suppressor p53. Finally, the covalent binding was confirmed by intact protein mass spectrometry (MS).

Results

In general, the purchased fragments turned out to be sufficiently soluble. Additionally, they covered a broad spectrum of reactivity. All investigated α-cyanoacrylamides/acrylates and all structurally confirmed epoxides turned out to be less reactive compounds, possibly due to steric hindrance and reversibility (for α-cyanoacrylamides/acrylates). The SNAr and vinyl sulfone fragments are either highly reactive or stable. DSF measurements with the different targets JNK3, USP7, and p53 identified reactive fragment hits causing a shift in the melting temperatures of the proteins. MS confirmed the covalent binding mode of all these fragments to USP7 and p53, while additionally identifying the SNAr-type electrophile SN002 as a mildly reactive covalent hit for p53.

Conclusion

The screening and target evaluation of the CovLib revealed first interesting hits. The highly cysteine-reactive fragments VS004, SN001, SN006, and SN007 covalently modify several target proteins and showed distinct shifts in the melting temperatures up to +5.1 °C and -9.1 °C.

QM/MM simulations of EFGR with afatinib reveal the role of the β-dimethylaminomethyl substitution

Acrylamides are the most commonly used warheads of targeted covalent inhibitors (TCIs) directed at cysteines; however, the reaction mechanisms of acrylamides in proteins remain controversial, particularly for those involving protonated or unreactive cysteines. Using the combined semiempirical quantum mechanics (QM)/molecular mechanics (MM) free energy simulations, we investigated the reaction between afatinib, the first TCI drug for cancer treatment, and Cys797 in the EGFR kinase. Afatinib contains a β-dimethylaminomethyl (β-DMAM) substitution which has been shown to enhance the intrinsic reactivity and potency against EGFR for related inhibitors. Two hypothesized reaction mechanisms were tested. Our data suggest that Cys797 becomes deprotonated in the presence of afatinib and the reaction proceeds via a classical Michael addition mechanism, with Asp800 stabilizing the ion-pair reactant state β-DMAM+/C797- and the transition state of the nucleophilic attack. Our work elucidates an important structure-activity relationship of acrylamides in proteins.

Machine Learning Models to Interrogate Proteome-Wide Covalent Ligandabilities Directed at Cysteines

Machine learning (ML) identification of covalently ligandable sites may accelerate targeted covalent inhibitor design and help expand the druggable proteome space. Here, we report the rigorous development and validation of the tree-based models and convolutional neural networks (CNNs) trained on a newly curated database (LigCys3D) of over 1000 liganded cysteines in nearly 800 proteins represented by over 10,000 three-dimensional structures in the protein data bank. The unseen tests yielded 94 and 93% area under the receiver operating characteristic curves for the tree models and CNNs, respectively. Based on the AlphaFold2 predicted structures, the ML models recapitulated the newly liganded cysteines in the PDB with over 90% recall values. To assist the community of covalent drug discoveries, we report the predicted ligandable cysteines in 392 human kinases and their locations in the sequence-aligned kinase structure, including the PH and SH2 domains. Furthermore, we disseminate a searchable online database LigCys3D (https://ligcys.computchem.org/) and a web prediction server DeepCys (https://deepcys.computchem.org/), both of which will be continuously updated and improved by including newly published experimental data. The present work represents an important step toward the ML-led integration of big genome data and structure models to annotate the human proteome space for the next-generation covalent drug discoveries.

Machine Learning Models to Interrogate Proteomewide Covalent Ligandabilities Directed at Cysteines

Machine learning (ML) identification of covalently ligandable sites may accelerate targeted covalent inhibitor design and help expand the druggable proteome space. Here we report the rigorous development and validation of the tree-based models and convolutional neural networks (CNNs) trained on a newly curated database (LigCys3D) of over 1,000 liganded cysteines in nearly 800 proteins represented by over 10,000 three-dimensional structures in the protein data bank. The unseen tests yielded 94% and 93% AUCs (area under the receiver operating characteristic curve) for the tree models and CNNs, respectively. Based on the AlphaFold2 predicted structures, the ML models recapitulated the newly liganded cysteines in the PDB with over 90% recall values. To assist the community of covalent drug discoveries, we report the predicted ligandable cysteines in 392 human kinases and their locations in the sequence-aligned kinase structure including the PH and SH2 domains. Furthermore, we disseminate a searchable online database LigCys3D (https://ligcys.computchem.org/) and a web prediction server DeepCys (https://deepcys.computchem.org/), both of which will be continuously updated and improved by including newly published experimental data. The present work represents a first step towards the ML-led integration of big genome data and structure models to annotate the human proteome space for the next-generation covalent drug discoveries.

Quantum Descriptors for Predicting and Understanding the Structure-Activity Relationships of Michael Acceptor Warheads

Predictive modeling and understanding of chemical warhead reactivities have the potential to accelerate targeted covalent drug discovery. Recently, the carbanion formation free energies as well as other ground-state electronic properties from density functional theory (DFT) calculations have been proposed as predictors of glutathione reactivities of Michael acceptors; however, no clear consensus exists. By profiling the thiol-Michael reactions of a diverse set of singly- and doubly-activated olefins, including several model warheads related to afatinib, here we reexamined the question of whether low-cost electronic properties can be used as predictors of reaction barriers. The electronic properties related to the carbanion intermediate were found to be strong predictors, e.g., the change in the Cβ charge accompanying carbanion formation. The least expensive reactant-only properties, the electrophilicity index, and the Cβ charge also show strong rank correlations, suggesting their utility as quantum descriptors. A second objective of the work is to clarify the effect of the β-dimethylaminomethyl (DMAM) substitution, which is incorporated in the warheads of several FDA-approved covalent drugs. Our data suggest that the β-DMAM substitution is cationic at neutral pH in solution and promotes acrylamide's intrinsic reactivity by enhancing the charge accumulation at Cα upon carbanion formation. In contrast, the inductive effect of the β-trimethylaminomethyl substitution is diminished due to steric hindrance. Together, these results reconcile the current views of the intrinsic reactivities of acrylamides and contribute to large-scale predictive modeling and an understanding of the structure-activity relationships of Michael acceptors for rational TCI design.

Integrated Covalent Drug Design Workflow Using Site Identification by Ligand Competitive Saturation

Covalent drug design is an important component in drug discovery. Traditional drugs interact with their target in a reversible equilibrium, while irreversible covalent drugs increase the drug-target interaction duration by forming a covalent bond with targeted residues and thus may offer a more effective therapeutic approach. To facilitate the design of this class of ligands, computational methods can be used to help identify reactive nucleophilic residues, frequently cysteines, on a target protein for covalent binding, to test various warhead groups for their potential reactivities, and to predict noncovalent contributions to binding that can facilitate drug-target interactions that are important for binding specificity. To further aid covalent drug design, we extended a functional group mapping approach based on explicit solvent all-atom molecular simulations (SILCS: site identification by ligand competitive saturation) that intrinsically considers protein flexibility, functional group, and protein desolvation along with functional group-protein interactions. Through docking of a library of representative warhead fragments using SILCS-Monte Carlo (SILCS-MC), reactive cysteines can be correctly identified for proteins being tested. Furthermore, a machine learning model was trained to quantify the effectiveness of various warhead groups for proteins using metrics from SILCS-MC as well as experimental model compound warhead reactivity data. The ability to rank covalent molecular binders with similar warheads using SILCS ligand grid free energy (LGFE) ranking was also tested for several proteins. Based on these tools, an integrated SILCS-based workflow was developed, named SILCS-Covalent, which can both qualitatively and quantitatively inform covalent drug discovery.

Deciphering the intrinsic dynamics of unphosphorylated IRAK4 kinase bound to type I and type II inhibitors

Interleukin-1 receptor-associated kinase 4 (IRAK4) is a vital protein involved in Toll-like and interleukin-1 receptor signal transduction. Several studies have reported regarding the crystal structure, dynamic properties, and interactions with inhibitors of the phosphorylated form of IRAK4. However, no dynamic properties of inhibitor-bound unphosphorylated IRAK4 have been previously studied. Herein, we report the intrinsic dynamics of unphosphorylated IRAK4 (uIRAK4) bound to type I and type II inhibitors. The corresponding apo and inhibitor-bound forms of uIRAK4 were subjected to three independent simulations of 500 ns (total 1.5 μs) each, and their trajectories were analyzed. The results indicated that all three systems were relatively stable, except for the type II inhibitor-bound form of uIRAK4, which exhibited less compact folding and higher solvent surface area. The intra-hydrogen bonds corroborated the structural deformation of the type-II inhibitor-bound complex, which could be attributed to the long molecular structure of the type-II inhibitor. Moreover, the type II inhibitor bound to uIRAK4 showed higher binding free energy with uIRAK4 than the type I inhibitor. The free energy landscape analysis showed a reorientation of Phe330 side chain from the DFG motif at different metastable states for all the systems. The intra-residual distance between residues Lys213, Glu233, Tyr262, and Phe330 suggests a functional interplay when the inhibitors are bound to uIRAK4, thereby hinting at their crucial role in the inhibition mechanism. Ultimately, the intrinsic dynamics study observed between type I/II inhibitor-bound forms of uIRAK4 may assist in better understanding the enzyme and designing therapeutic compounds.

GPU-Accelerated All-Atom Particle-Mesh Ewald Continuous Constant pH Molecular Dynamics in Amber

Constant pH molecular dynamics (MD) simulations sample protonation states on the fly according to the conformational environment and user specified pH conditions; however, the current accuracy is limited due to the use of implicit-solvent models or a hybrid solvent scheme. Here, we report the first GPU-accelerated implementation, parametrization, and validation of the all-atom continuous constant pH MD (CpHMD) method with particle-mesh Ewald (PME) electrostatics in the Amber22 pmemd.cuda engine. The titration parameters for Asp, Glu, His, Cys, and Lys were derived for the CHARMM c22 and Amber ff14sb and ff19sb force fields. We then evaluated the PME-CpHMD method using the asynchronous pH replica-exchange titration simulations with the c22 force field for six benchmark proteins, including BBL, hen egg white lysozyme (HEWL), staphylococcal nuclease (SNase), thioredoxin, ribonuclease A (RNaseA), and human muscle creatine kinase (HMCK). The root-mean-square deviation from the experimental pK_a's of Asp, Glu, His, and Cys is 0.76 pH units, and the Pearson's correlation coefficient for the pK_a shifts with respect to model values is 0.80. We demonstrated that a finite-size correction or much enlarged simulation box size can remove a systematic error of the calculated pK_a's and improve agreement with experiment. Importantly, the simulations captured the relevant biology in several challenging cases, e.g., the titration order of the catalytic dyad Glu35/Asp52 in HEWL and the coupled residues Asp19/Asp21 in SNase, the large pK_a upshift of the deeply buried catalytic Asp26 in thioredoxin, and the large pK_a downshift of the deeply buried catalytic Cys283 in HMCK. We anticipate that PME-CpHMD will offer proper pH control to improve the accuracies of MD simulations and enable mechanistic studies of proton-coupled dynamical processes that are ubiquitous in biology but remain poorly understood due to the lack of experimental tools and limitation of current MD simulations.

Electrostatics in Computational Biophysics and Its Implications for Disease Effects

This review outlines the role of electrostatics in computational molecular biophysics and its implication in altering wild-type characteristics of biological macromolecules, and thus the contribution of electrostatics to disease mechanisms. The work is not intended to review existing computational approaches or to propose further developments. Instead, it summarizes the outcomes of relevant studies and provides a generalized classification of major mechanisms that involve electrostatic effects in both wild-type and mutant biological macromolecules. It emphasizes the complex role of electrostatics in molecular biophysics, such that the long range of electrostatic interactions causes them to dominate all other forces at distances larger than several Angstroms, while at the same time, the alteration of short-range wild-type electrostatic pairwise interactions can have pronounced effects as well. Because of this dual nature of electrostatic interactions, being dominant at long-range and being very specific at short-range, their implications for wild-type structure and function are quite pronounced. Therefore, any disruption of the complex electrostatic network of interactions may abolish wild-type functionality and could be the dominant factor contributing to pathogenicity. However, we also outline that due to the plasticity of biological macromolecules, the effect of amino acid mutation may be reduced, and thus a charge deletion or insertion may not necessarily be deleterious.

A Guide to the Continuous Constant pH Molecular Dynamics Methods in Amber and CHARMM [Article v1.0]

Like temperature and pressure, solution pH is an important environmental variable in biomolecular simulations. Virtually all proteins depend on pH to maintain their structure and function. In conventional molecular dynamics (MD) simulations of proteins, pH is implicitly accounted for by assigning and fixing protonation states of titratable sidechains. This is a significant limitation, as the assigned protonation states may be wrong and they may change during dynamics. In this tutorial, we guide the reader in learning and using the various continuous constant pH MD methods in Amber and CHARMM packages, which have been applied to predict pK a values and elucidate proton-coupled conformational dynamics of a variety of proteins including enzymes and membrane transporters.