Covalent Docking of Large Libraries for the Discovery of Chemical Probes
Chemical probes that form a covalent bond with a protein target often show enhanced selectivity, potency, and utility for biological studies. Despite these advantages, protein-reactive compounds are usually avoided in high-throughput screening campaigns. Here we describe a general method (DOCKovalent) for screening large virtual libraries of electrophilic small molecules. We apply this method prospectively to discover reversible covalent fragments that target distinct protein nucleophiles, including the catalytic serine of AmpC \xce\xb2-lactamase and noncatalytic cysteines in RSK2, MSK1, and JAK3 kinases. We identify submicromolar to low-nanomolar hits with high ligand efficiency, cellular activity and selectivity, including the first reported reversible covalent inhibitors of JAK3. Crystal structures of inhibitor complexes with AmpC and RSK2 confirm the docking predictions and guide further optimization. As covalent virtual screening may have broad utility for the rapid discovery of chemical probes, we have made the method freely available through an automated web server (http://covalent.docking.org).
Introduction<!>Overview of the method<!>Retrospective assessment of covalent docking<!>New boronic acid inhibitors of AmpC \xce\xb2-lactamase<!>New cyanoacrylamide inhibitors for RSK2 and MSK1 kinases<!>Selective, reversible covalent inhibitors of JAK3 kinase<!>Discussion<!>Ligand Generation<!>Availability<!>Virtual ligand libraries<!>Covalent Virtual Screening<!>Selection criteria<!>RMSD calculations<!>Enzymology<!>Reversibility<!>Selectivity<!>In vitro proteasome activity assay<!>Microbiology<!>Crystallography<!>RSK2 kinase assays<!>JAK3 kinase assays<!>JAK kinase selectivity<!>Assay demonstrating dissociation of covalent inhibitors from RSK2 CTD and JAK3<!>Cell-based assay and Western blotting<!>Compound sources<!>Structural Data Deposition
<p>Small-molecule modulators of protein activity are useful as tools for investigating biology and as leads for drug discovery. Development of genuinely useful probes typically involves iterative rounds of medicinal chemistry to optimize the potency and selectivity of initial active molecules. An effective strategy for enhancing both properties is via covalent bond formation with a nucleophilic residue that is specific to a target of interest and ideally absent from off-targets. Such covalent-acting chemical probes have increasingly been used in proteome-wide target identification 1, imaging 2, and for finding inhibitors with high specificity among related enzymes and enzyme isoforms 3,4. Covalent drugs 5 and natural products 6 are also well known. A challenge in developing covalent probes is identifying reactive functional groups ("warheads") that do not make the molecule so reactive as to be promiscuous. Less recognized is the challenge of screening a wide variety of scaffolds for optimal presentation of such reactive functionality.</p><p>The most widely used technique for novel ligand discovery is high-throughput screening (HTS), and one could potentially screen extant libraries for new small molecules that react covalently. However, protein-reactive compounds are rarely screened 7 and are typically avoided in HTS 8 or flagged as artifacts due to concerns about promiscuous activity 9. Whereas this is sensible for drug discovery, it removes potential starting points for covalent chemical probes 10.</p><p>Covalent ligands can target either catalytically essential nucleophiles such as those in serine and cysteine hydrolases, or noncatalytic nucleophiles, usually cysteine, found in small-molecule binding sites on proteins, including proteins without enzymatic activity (e.g., GPCRs and nuclear receptors) 11. Depending on the electrophile and nucleophile, they can bind reversibly or irreversibly. In all cases, specific noncovalent interactions contributed by the scaffold are critical for orienting the electrophile relative to the protein nucleophile, thereby increasing the rate (and stability, in the case of reversible covalent ligands) and selectivity of covalent bond formation. A key unsolved problem in the discovery of covalent probes is how to identify a protein-binding scaffold that optimally orients the electrophile, while minimizing the number of compounds that must be synthesized and tested.</p><p>In principle, structure-based docking screens 12,13 can address the gap left by HTS and its libraries. Given the structure of a protein target, docking programs computationally screen large compound libraries for molecules predicted to bind favorably within a defined binding site. The technique has been widely used for the discovery of reversible, non-covalent ligands 14,15. To date, there have been few docking screens for covalent ligands. Key obstacles include combining classical non-covalent scoring with covalent restraints and bond energies, and developing compound libraries suited to covalent modification of proteins. For the compound libraries, one would prefer not only commercially available electrophiles, but also the ability to design new, readily synthesizable molecules bearing a particular electrophile. Recently, there has been encouraging progress in developing covalent docking methods, but these have been restricted to retrospective recapitulation of covalent complexes 16–18, or to screens of a few hundred compounds 19,20; we are unaware of prospective, large-scale covalent docking screens to find novel inhibitors, or of any such screens against targets for which covalent ligands are unprecedented.</p><p>Here, we adapt the non-covalent docking program DOCK3.6 to large-scale, covalent virtual screening of electrophilic small molecules, including low-molecular weight electrophilic fragments (Fig. 1). Nine libraries of ligands bearing different electrophiles, amounting to over 650,000 commercially available or synthetically accessible small molecules, are developed for use with the method. We used the method (DOCKovalent) to prospectively screen compound libraries against three targets of therapeutic interest: AmpC β-lactamase, RSK2 kinase and JAK3 kinase. Multiple potent, reversible covalent inhibitors were found against all three targets. X-ray crystal structures of predicted ligands, and the occasional false negative, illuminated not only the method's ability to prospectively identify ligands and to predict their structures, but also its limitations. Several of the new covalent ligands were tested in cell culture experiments that established biological efficacy and target engagement. To ensure that the method may be used by a broad community, it has been made available on an easy to use web server (http://covalent.docking.org).</p><!><p>We begin by constructing large virtual libraries of electrophiles, either commercially available or synthetically accessible in 1–2 steps. We created libraries of well-studied electrophiles including: α,β-unsaturated carbonyls, aldehydes, boronic acids, α-cyanoacrylamides, alkyl halides, carbamates, α-ketoamides, and epoxides (Supplementary Results, Supplementary Fig. 1, Supplementary Table 1); other electrophilic chemotypes, such as vinyl sulfones, sulfonyl fluorides, 2-chloropyridines, and cyanopyrimidines, are also imaginable. All stereoisomers, protonation states and conformations of the covalent adduct are pre-generated for each ligand, enabling rapid docking of the library to any target (Fig. 1).</p><p>For each ligand, DOCKovalent exhaustively samples all poses and ligand conformations with respect to the covalent bond to the target nucleophile, constrained by ideal bond lengths and angles (Fig. 1 and Supplementary Fig. 2). The nucleophile is immobile during the sampling, and a separate screen should be run for each likely rotamer of the nucleophile. Each sampled conformation is scored using the physics-based scoring function in DOCK3.6 21, which evaluates the ligand's van der Waals and electrostatic interactions, and corrects for its desolvation. Using this scoring function, the entire library is ranked from most to least favorable. The top 1–3% of the ranked list is inspected for mis-docked ligands, which are common in large-scale docking screen, and molecules with incorrect ionization states, tautomers, or strained conformations are removed. The remaining molecules are prioritized for experimental testing based on their availability or synthetic accessibility, the presence of unprecedented chemotypes, and diversity of chemical structure.</p><!><p>We first tested method's ability to find known covalent ligands in five retrospective screens against four targets, and to recapitulate geometries for a previously published benchmark of covalent ligand complexes 16. DOCKovalent performed well in pose recapitulation (Supplementary Table 2,3) and in four of the five retrospective screens (Supplementary Note 1 and Supplementary Table 4). We thus turned to prosecuting prospective screens for new reversible covalent ligands for three enzymes.</p><!><p>AmpC β-lactamase is the leading cause of resistance to cephalosporin antibiotics in clinical settings 22, and several new β-lactamase inhibitors are in clinical trials 23. Boronic acids inhibit AmpC by forming a reversible covalent adduct with its active-site nucleophilic serine (Ser64). We first assessed the ability of our covalent docking method to recapitulate known boronic acid complexes with AmpC. In 15 of 23 cases, the ligand pose was accurately recovered to less than 2 Å RMSD (Supplementary Table 5 and Supplementary Fig. 3). Surprisingly, a relatively simple compound, m-aminophenylboronic acid (MAPB), failed this retrospective test. The boronic acid in the published MAPB/AmpC structure (PDB: 3BLS; Fig. 2a) adopted a pose that differed from other boronic acid complexes with the enzyme, and from our other docking poses. To investigate this discrepancy, we re-determined the crystal structure of AmpC in complex with MAPB. In this 1.65 Å resolution x-ray structure we observed unambiguous density for the ligand that corresponds to the docking prediction (0.46 Å RMSD, Fig. 2a).</p><p>Encouraged by these results, we used the method to covalently dock a library of 23,000 commercially available boronic acids against AmpC. Among the top-ranked 4.5% of the library, we sought boronic acids with scaffolds that had not been tested previously against AmpC. Five such compounds ranked between 11 and 646 of the docked library (top 0.04 to 2%) and were purchased and tested (compounds 1–5;
Fig. 2b). An additional, lower ranking compound (6) was purchased as a proxy for structurally related, high-ranking predictions that were commercially unavailable. Five of the six compounds inhibited AmpC with Ki values ranging from 40 nM to 3.55 μM (Table 1), and three of these exhibited submicromolar potency (ligand efficiency, LE 0.38–0.66, Table 1; see docking predictions in Fig. 2c–d and Supplementary Fig. 4). None of the five inhibitors resemble a known boronic acid inhibitor of AmpC (Tanimoto coefficients < 0.3 to AmpC boronic acids in ChEMBL, using ECFP4-based fingerprints).</p><p>A 1.74 Å crystal structure of compound 3, the most potent inhibitor from our initial set of six compounds (Ki = 40 nM), confirmed the docking pose prediction (1.38 Å RMSD, Fig. 2c). The boronic acid occupies the oxyanion hole formed by the backbone amides of Ala318 and Ser64 and hydrogen bonds with Tyr150. More importantly, noncovalent interactions between the scaffold and AmpC were well predicted. The pyrazole N2 accepts hydrogen bonds from Asn152 and Gln120, while the phenyl moiety stacks against Tyr221. The only substantial discrepancy between the docking prediction and the crystal structure is the position of the distal chlorine atom. This may reflect the presence of a conserved water network in the active site, which was not included in the calculation (Fig. 2c–e).</p><p>Several compounds related to pyrazole 3 were also highly ranked by docking. We therefore purchased seven additional pyrazole boronic acids (Compounds 7–13; Supplementary Fig. 5), one of which showed four-fold greater potency (compound 7, Ki = 10 nM; Fig. 2b). In a crystal structure we determined, compound 7 binds to AmpC in essentially the same manner as compound 3 (Fig. 2c–d). Its increased affinity may arise from a favorable interaction between the new pyrimidine ring and the conserved water network observed in both complexes, or to a stronger electrostatic interaction with the carbonyl of Gln120. Ultimately, low-nanomolar inhibitors were obtained by purchasing only 13 compounds.</p><p>We characterized the selectivity of the four most potent compounds (2, 3, 5, and 7) by testing them against three common serine proteases known to bind boronic acids: trypsin, elastase, and α-chymotrypsin 24, and against the yeast 20S proteasome. The new AmpC inhibitors typically showed >1000-fold selectivity vs. the serine proteases, and none inhibited the 20S proteasome greater than 20% at 100 μM (Supplementary Table 6 & Supplementary Fig. 6). An exception was compound 3, which inhibited α-chymotrypsin with a Ki of 300 nM. However, pyrimidine 7, the most potent AmpC inhibitor, showed 104-fold selectivity over α-chymotrypsin and >105-fold selectivity over trypsin and elastase.</p><p>A concern when screening electrophilic compounds is that the electrophile will be so reactive that most compounds in the library will bind the target. To control for this, we tested five boronic acids from the bottom of the ranked docking list (Compounds 14–18; Fig. 2b). We avoided trivial non-binders, selecting only those molecules for which the docking program found a non-clashing pose. Four of the five predicted non-binders showed less than 10% AmpC inhibition at 10 μM, consistent with prediction (Supplementary Table 7). Compound 14, however, did have measurable activity (Ki = 3.2 μM).</p><p>To investigate the origins of this docking false negative, we determined the crystal structure of 14 in complex with AmpC, which revealed unambiguous ligand density in a pose different from the predicted docking model (Fig. 2e). To accommodate the observed geometry, an active-site loop (L117-Q120) changes conformation, with Leu119 adopting a new rotamer and the loop moving by 0.7 Å (Cα RMSD, Fig. 2e). This binding mode is incompatible with the AmpC structure used for docking and highlights a caveat of our approach: to enable fast screening of large libraries, we treat the receptor as fixed. The new loop conformation is unique across 23 AmpC structures (Fig. 2e; Supplementary Fig. 7).</p><p>We next tested whether the new boronic acid inhibitors could reverse antibiotic resistance in bacteria that express AmpC. We determined the minimum inhibitory concentration (MIC) of cefotaxime, alone or in combination with inhibitors, against eight clinical isolates resistant to third-generation cephalosporins (Table 1). Consistent with the enzymatic assays, compound 7 was the most potent at reversing antibiotic resistance, lowering the MIC for six strains to ≤2 μg mL, the empirical threshold for resistance defined by the Clinical and Laboratory Standards Institute 25. None of the compounds had significant antibiotic activity in the absence of cefotaxime (Supplementary Table 8).</p><!><p>The C-terminal kinase domains (CTD) of p90 ribosomal protein S6 kinase-2 (RSK2) and the closely related paralog, mitogen- and stress-activated kinase-1 (MSK1), contain a noncatalytic active-site cysteine shared by only 11 of the 518 human protein kinases. Starting with an established kinase inhibitor scaffold, we previously designed irreversible 3 and reversible 26,27 covalent inhibitors that target this cysteine (Cys436 in RSK2). To achieve reversible covalent inhibition, we exploited the atypical reactivity of cyanoacrylamide Michael acceptors, which react rapidly and reversibly with cysteine thiols at physiological pH. Both RSK2 and MSK1 are attractive therapeutic targets implicated in tumor metastasis 28,29, neurodegeneration 30, and atherosclerosis 31, among other pathological conditions. We sought novel RSK2/MSK1 inhibitors by covalent docking screens of thousands of cyanoacrylamide fragments.</p><p>As an initial blind test, we used the method to predict the poses of two cyanoacrylamide fragments bound to RSK2, prior to determining the crystal structures. The predicted binding modes anticipated the experimental structures to 1.93 Å and 1.56 Å RMSD (Fig. 3a and 3b). Retrospective docking of two larger cyanoacrylamides also recapitulated their crystal structures (0.66 Å and 1.52 Å RMSD; Supplementary Fig. 8a–b). In each prediction, the scaffold portion of the molecule, which forms critical non-covalent interactions with RSK2, closely matched the x-ray structures (0.91–1.36 Å RMSD).</p><p>Encouraged by these results, we used covalent docking to screen for novel cyanoacrylamide inhibitors. Cyanoacrylamide fragments are rare in commercial collections (602 out of 474,770 of the "fragments in-stock" in ZINC 32). However, β-substituted cyanoacrylamides can be synthesized in one step by condensing aldehydes with cyanoacetamide (Supplementary Fig. 9). We therefore assembled ~12,000 aldehyde fragments from ZINC 32 (molecular weight < 250 Da). These were converted in silico to generate a virtual library of cyanoacrylamide fragments.</p><p>We docked this library against Cys436 of RSK2. After manually inspecting the top-ranked compounds for novelty, diversity, and accessibility, we pursued eight virtual cyanoacrylamide fragments ranked between 96 and 391 (top 3%; Compounds 19–26;
Fig. 3c). The corresponding aldehydes were purchased and converted to the cyanoacrylamides, which were tested against wild-type RSK2 and the T493M gatekeeper mutant (Table 2). We have previously used this mutant as a biochemical surrogate for MSK1, as MSK1 CTD kinase activity has yet to be reconstituted in vitro
27. Five of the eight high-ranking cyanoacrylamides inhibited RSK2 with IC50 values less than 10 μM in the presence of competing 0.1 mM ATP and 10 mM reduced glutathione. Pyridine 21 and pyrazole 24 (Fig. 3d–e) were the most potent against wild-type RSK2 and inhibited the T493M mutant with submicromolar potency (IC50 430 nM and 370 nM, respectively).</p><p>We determined the co-crystal structure of compound 24 bound to T493M RSK2 at 3.0 Å resolution. Even at this modest resolution, the electron density allowed unambiguous modeling of the phenylpyrazole fragment and the covalent bond to Cys436 (Supplementary Fig. 10). The crystallographic structure superposed well with the docking prediction, with an RMSD of 0.86 Å over the phenylpyrazole fragment (Fig. 3d). The pyrazole forms two hydrogen bonds with the hinge region (N1 to the carbonyl of Glu494 and N2 to the amide NH of Met496). This ring also packs edgewise against the methionine gatekeeper, and the additional van der Waals contacts may explain the enhanced potency for the T493M mutant.</p><p>We tested 21 and 24 for activity in mammalian cells stimulated with phorbol myristate acetate (PMA), which activates kinase cascades upstream of RSK2 and MSK1. Both compounds inhibited the activating autophosphorylation of wild-type MSK1 (EC50 < 5μM; Fig. 3f). While these compounds inhibit WT RSK2 less potently in cells, quantitation of the normalized pS380 signal reveals dose-dependent inhibition (EC50 ~ 20 uM). A cysteine to valine substitution in the C-terminal kinase domain of MSK1 (C458V) conferred complete resistance to both inhibitors, consistent with on-target efficacy. Hence, the unoptimized cyanoacrylamide fragments identified by covalent docking inactivate the target kinases in cells without affecting upstream kinases such as RAF, MEK, and ERK.</p><p>The docking pose of 21 (Fig. 3e) suggested the possibility of improving its potency by fusing a pyrrole ring, an H-bond donor, to either a pyridine H-bond acceptor, as in 21, or a pyrimidine ring. Exemplified by pyrrolopyrimidine 27 (Fig. 3g), this would bury more hydrophobic surface area and form an additional hydrogen bond with the hinge region, as suggested by covalent docking to RSK2 (Fig. 3h). Compound 27 was not included in the original screen, as the corresponding aldehyde is not commercially available. On synthesis and testing, compound 27 inhibited WT RSK2 kinase in vitro with an IC50 of 42 nM, over 25-fold better than 21 (Fig. 3g). Correspondingly, 27 was substantially more potent than 21 in cells, blocking MSK1 autophosphorylation with an EC50 < 1 μM (Fig. 3i).</p><!><p>Members of the Janus kinase family, comprised of JAK1, JAK2, JAK3, and TYK2, are essential for signaling downstream of many cytokine receptors 33. JAK3 is expressed predominantly in immune cells and is a potential therapeutic target for autoimmune diseases like rheumatoid arthritis (RA) 34. A pan-JAK inhibitor, tofacitinib 35, was recently approved for RA, but it suffers from adverse effects such as elevated liver enzymes and LDL cholesterol 36. Selective JAK3 inhibitors may avoid such toxicities, and moreover, could help illuminate JAK3-specific roles in cytokine signaling. To date, development of selective JAK3 inhibitors has been hampered by the high sequence identity among JAK-family kinases 37. JAK3 contains a solvent-exposed cysteine residue just outside the ATP binding site (Cys909), which is not found in JAK1, JAK2, or TYK2, and is present in only nine other human kinases. We used DOCKovalent in an effort to find the first reversible covalent inhibitors of JAK3, which might be expected to have specificity over closely related JAK kinases that lack Cys909.</p><p>The vector from Cys909 to the hinge differs greatly from the previously targeted Cys436 of RSK2. A preliminary screen of the virtual cyanoacrylamide fragment library developed initially for RSK2 suggested that greater diversity and perhaps larger fragments would be required to engage both Cys909 and the hinge of JAK3. Inspired by the simple two-step synthesis of 27, we designed a combinatorial virtual library based on two synthetic transformations: a Suzuki-Miyaura cross-coupling reaction between an aryl or heteroaryl bromide and an aldehyde-containing boronic acid, followed by a Knoevenagel condensation of the aldehyde with cyanoacetamide. We selected 50 commercially available boronic acids and 4,400 aryl bromides, which were converted to their corresponding products in silico. This approach afforded a diverse virtual library of 220,000 hetero-biaryl cyanoacrylamides (Supplementary Table 1, Supplementary Fig. 11), which was screened against eight JAK3 crystal structures.</p><p>We purchased eight arylbromides and three boronic acids common among the top 0.2% of the docked library (Supplementary Table 9). From these building blocks, we synthesized 15 inhibitors (Compounds 28–42; Fig. 4a). Compound 42 was prepared from a commercially available aldehyde. Each compound was initially tested against JAK3 at two concentrations. Nine of the 15 compounds inhibited the kinase by > 50% at 5 μM, and five maintained > 50% inhibition at 1 μM (Fig. 4b).</p><p>We focused on the two most potent compounds 31 and 33 (for docking poses, see Supplementary Fig. 12). To assess their selectivity, we measured full dose-response curves against JAK3 and against the other three JAK-family kinases lacking Cys909. Cyanoacrylamides 31 and 33 inhibited JAK3 with IC50 values of 49 nM and 93 nM, respectively, but neither inhibitor affected any other JAK kinase at concentrations up to 10 μM (Fig. 4c). This marked selectivity for JAK3 may be largely attributed to covalent bond formation with the target cysteine. By virtue of this covalent targeting of a noncatalytic cysteine, 31 and 33 are among the most selective JAK3 inhibitors reported to date 37. We further assessed the selectivity of 31 against nine additional human kinases that have a cysteine at the equivalent position to JAK3. 31 potently inhibited three of the nine kinases (IC50= 22 nM, 44 nM and 221 nM for BLK, ERB-B4 and ITK, respectively) but had at least 30-fold specificity for JAK3 over the remaining six kinases (IC50 > 1 μM; Supplementary Fig. 13).</p><!><p>Covalent probes play a crucial role in chemical biology 1–4, yet electrophilic molecules that might serve as initial hits for developing such probes have largely been expunged from empirical screening libraries. Accordingly, we sought to enable large-scale screening for covalent probe molecules via structure-based docking. Although the method has limitations, its utility is supported by the discovery of new chemotypes in three prospective docking screens. For AmpC, the new inhibitors bear little topological resemblance with previously known inhibitors — with ECFP4-based Tanimoto coefficients < 0.3 to the most similar previous inhibitor — despite intense study of boronic acids against this enzyme over the last 30 years. For RSK2/MSK1, docking identified the active phenylpyridine 21, even though isomeric phenylpyridines were previously found to be inactive 27. The predicted pose of 21 supported the design of the unique pyrrolopyrimidine variant 27. Finally, the 4-phenylindazole scaffold of the newly discovered JAK3 inhibitor 31 has little precedent among kinase inhibitors; a related scaffold found in linifanib binds in a completely different orientation to VEGFR2 (PDB: 1YWN).</p><p>For all targets, the hit rates were high — five of six molecules predicted and tested for β-lactamase (Table 1), five of eight for RSK2 kinase (Table 2) and nine of fifteen for JAK3 kinase (Fig. 4). These hit rates did not reflect broad promiscuity on the part of the electrophiles. The boronic acids active against β-lactamase were typically inactive against related serine proteases (Supplementary Table 6) and against the proteasome (Supplementary Fig. 6), and most poorly ranked boronic acids were inactive against the enzyme (Supplementary Table 7). The one exception to this was shown by crystallography to reflect an unexpected and previously unobserved conformational change. Similarly, MSK1 inhibitors were selective in cells, being inactive against the Cys458→Val substitution of MSK1, as well as the upstream kinases, MEK and ERK, which also have noncatalytic cysteines in their active sites. Finally, the most potent JAK3 inhibitors were virtually inactive against other JAK kinases, underscoring the advantages of covalently targeting poorly conserved noncatalytic cysteines. Moreover, compound 31 showed selectivity over six of the nine human kinases bearing a homologous cysteine. Given its low molecular weight (288 Da) and simple architecture, accessible in two steps from inexpensive building blocks, the selectivity of 31 among these cysteine-containing kinases is remarkably high, comparing favorably with the approved covalent drug, ibrutinib 38, which targets the same cysteine in BTK.</p><p>Whereas most of the initial hits do not have the potency typically associated with optimized chemical probes, the best examples from the AmpC and RSK/MSK screens were sufficiently potent and selective to be active in cell-based assays. Hits from all screens were remarkably ligand efficient, with LE values as high as 0.66. Moreover, the initial hits could be improved, as shown in two cases. Boronic acid 7 has four-fold higher potency relative to the initial docking hit, 3. Indeed, this inhibitor appears to have the highest ligand efficiency of any serine β-lactamase inhibitor 24. As for RSK2/MSK1, compound 27 was improved 25-fold versus its parent screening hit, had sub-micromolar cellular activity, and is one of the most ligand-efficient RSK2 inhibitors known. The origins of the selectivity and affinities of these ligands are captured by their docking poses, which superpose well with the subsequently determined x-ray structures. The high fidelity of the docking to the x-ray structures suggests that the covalently docked poses alone may guide ligand optimization in cases where obtaining experimental structures is challenging, as demonstrated by 27.</p><p>An innovation introduced in this study, one with few precedents in docking as a field, is the development and experimental testing of truly virtual, readily synthesizable, molecular libraries. The few commercially available cyanoacrylamides necessitated the design of larger and more diverse libraries, which is likely the case for other classes of reactive molecules as well.</p><p>One of our concerns at the outset was that, in mixing covalent and non-covalent scoring terms, the former would overwhelm the latter 39. To avoid this, we did not explicitly score the covalent bond energy of the adducts, but rather imposed distance and angle constraints on bond formation and then scored the docked poses by our standard physics-based scoring function, ignoring the covalent bond. An advantage is that the scoring function is not dominated by the covalent bond to the electrophile, and instead favors structural complementarity and specificity provided by noncovalent interactions with the scaffold. The approach is especially well suited to the reversible, and hence thermodynamically driven, inhibitors discovered here (our best compounds were explicitly shown to be rapidly reversible by dilution; Supplementary Fig. 14). Still, our docking method can be extended to model high-energy intermediate states of the covalent bond-forming reaction 40, as reflected in the high enrichment in retrospective docking screens against FAAH and AChE (Supplementary Table 4), where the carbamates were modeled as high-energy intermediates (Supplementary Fig. 1). The method also may be suited to irreversible Michael acceptors, such as acrylamide-based inhibitors of the EGFR kinase, where retrospective enrichment was also high (Supplementary Table 4).</p><p>The failure to explicitly model the energy of covalent bond formation is one of several methodological gaps that merit mention. Without considering the covalent bond energy, we cannot directly compare molecules bearing different electrophiles. Thus, currently each library of different electrophiles is docked and ranked separately. Incorporating quantum-level approaches to predict the reactivity of specific scaffold/electrophile combinations may be a direction forward 41. Even within the current constraint-based approach, sampling can be improved by relaxing around ligand dihedral angles, and perhaps by accounting for bonded energies. Lastly, as demonstrated by compound 14 with AmpC, allowing for receptor flexibility can be essential for accurate modeling of receptor/ligand interactions. This remains an area of active research 12,42.</p><p>These caveats should not obscure the potential benefits offered by the method. We hope it will begin to address what has been a substantial gap in our ability to screen widely for useful covalent chemical probes, and we have made it and its associated libraries freely available through a general use, public access web-portal (http://covalent.docking.org).</p><!><p>Ligand flexibility is sampled by generating ligand conformations prior to docking. Given a SMILES string of a ligand with a specific electrophile, we use the OEChem library (OEChem TK 1.7.4; Openeye Scientific Software: Santa Fe, NM. http://www.eyesopen.com) to convert the ligand to its final, reacted form (Supplementary Fig. 1). The receptor's nucleophilic atom involved in the covalent bond is represented by a dummy atom (silicon, for technical reasons). Following the generation of the 'reacted' electrophile, the ligand's 3D structures and stereoisomers are built by Corina 43 (Molecular Networks, Erlangen, Germany), and then protonated and tautomerized by EPIK (Schrodinger software, Catsville NY). Partial atomic charges and solvation energies are calculated for each of these structures with AMSOL 44. The electrophile serves as a starting rigid body fragment, and conformations are generated using Omega 45 (Omega parameters: EnergyWindow=30.0; MaxConfs=10,000; RMSThreshold=0.5). The collection of pre-generated ligand conformations in the reacted state is saved to a DOCK-readable flexibase format file.</p><p>Sampling of ligand poses within the protein binding-site is restricted to exhaustive ligand placement with respect to the covalent bond (Supplementary Fig. 2). The covalent attachment point is sampled in steps of 20° around the terminal dihedral of the nucleophilic side chain. Based on the electrophile geometry determined during ligand generation, and user provided parameters, the vectors of the covalent bond from the ligand and receptor sides are aligned and the ligand is rotated around this vector in 20° steps. For each placement, all of the pre-generated ligand conformations are scored and the score for the best pose is saved. This process is repeated for different values of the covalent bond length and angles, centered on ideal values (Supplementary Fig. 2). The magnitude of deviation from the ideal values, as well as the step sizes, are user specified.</p><p>Scoring is performed as previously described, using pre-calculated van der-Waals, electrostatic, and ligand solvent-excluded desolvation grids, correcting for ligand desolvation 21. Receptor structures were prepared using an automated procedure as described in 46 using DELPHI 47 for electrostatics. The ligand's electrophilic atom participating in the bond is omitted from the overall ligand score.</p><!><p>As noted, the method is accessible through a public web-server (http://covalent.docking.org), and for download as part of the next DOCK3.x release (http://dock.compbio.ucsf.edu/).</p><!><p>For the curation of ligand libraries, the electrophiles were represented by SMARTS regular expressions (Supplementary Table 1). The Full, Lead-like, and Fragment-like subsets of the ZINC database of commercially available molecules were filtered using these patterns to identify electrophile-bearing molecules 32. Known inhibitors for the five retrospective virtual screens were collected from ChEMBL14 48 for AChE, FAAH (carbamates with <1uM activity), and NS3 (α-ketoamides with <100nM activity). For EGFR ligands, we used the same ligand set as in 16, and FAAH boronic acid inhibitors (Ki < 10 nM) were collected from 49. The best receptor structure was selected based on its ability to enrich known ligands.</p><p>Known AmpC binders were defined as having < 1uM activity in ChEMBL16 48. Compound similarity was calculated using ECFP4-based Tanimoto coefficients 50 as implemented in Pipeline-Pilot version 6.1 (SciTegic Inc. San Diego, CA). Boron is not parameterized in several of our ligand generation programs, to overcome this in the construction of our boronic acids library, we searched ZINC's commercial catalogs for boronic acids, replacing the boron with a carbon during ligand generation. This doesn't affect ranking, as the energy of the ligand's bonding atom is omitted from the docking score.</p><p>The one step cyanoacrylamide based library is based on filtering ZINC's fragment-like subset for fragments containing an aldehyde. The two step, Suzuki-Miyaura based cyanoacrylamide library is made according to (Supplementary Fig. 11) by combining any of 50 commercially available aldehyde containing boronic acids with 4,397 aryl bromide fragments from ZINC's fragments-in-stock subset, containing one of the following motifs: 'ncN', [NX3][CX3](=[OX1])[#6], or a pyridine, pyrimidine, pyrrole, pyrazine, pyrazole, or triazole.</p><!><p>For the AmpC retrospective pose recapitulation benchmark (Supplementary Table 5) and prospective screen, PDB: 4E3N was used as the receptor. The flexible side chain of Gln120, was truncated at Cβ. The backbone of residues Ala318 and Ser64 and side chains of Asn152 and Tyr150 were polarized to emphasize the electrostatic interactions with the boronic acid. Covalent bond sampling parameters were set to d=1.5±0.1Å a=116.0°±5° and b=109.5°±5°. (Step size of 0.05Å and 1°; Supplementary Fig. 2). Compounds 1, 3, 4, and 14–18 were selected for testing based on this screen. Following the determination of the crystal structure of AmpC in complex with 3, a second screen was performed using the new structure (PDB: 4LV1; in this screen, the bond length was fixed to d=1.6Å, and Gln120 was not truncated). Compounds 2, 5 and 6 were selected for testing. Run time for the entire screen was 227 CPU hours (ran in parallel on an 800-CPU computer cluster; Wall time < 1 hour). The RSK2 screen used PDB: 4D9T for the receptor structure. Given that nearly all kinase inhibitors form one or more hydrogen bonds with the hinge region, the backbone amides of the kinase hinge residues (Glu494 and Met496) were polarized to emphasize hinge-binding hydrogen bonds. A technique we have long used. Covalent bond sampling parameters were set to d=2.0Å a=109.5°±10° and b=109.5°±10° (Step size of 2.5°; Supplementary Fig. 2). This bond length was chosen based on the two available RSK2 cyanoacrylamide complexes (PDBs: 4D9U,4D9T). While longer than a typical thioether by about 0.2Å, this helps ensure that the scores are dominated by the non-covalent terms, minimizing hard van der Waals repulsion in the region of the new bond. We have not attempted to optimize this term. Calculation time for the docking screen was 103 CPU hours, elapsed wall time was less than 1 hour owing to use of a cluster. JAK3 docking was performed using eight available crystal structures PDB: 1YVJ, 3LXK, 3LXL, 3PJC, 4HVD, 4HVG, 4HVH, and 4HVI. For each structure we used both the native cysteine rotamer as well as an alternative rotamer (χ1=−60°). Docking to the native rotamer produced few plausible results, likely due to steric clashes of the electrophile with nearby residues, and only one candidate was chosen based on this rotamer (Supplementary Table 9). The backbone amides of the hinge residues (Glu903 and Leu905) were polarized similarly to RSK2. Covalent bond sampling parameters were set to d=1.8Å a=109.5°±10° and b=109.5°±10° (Step size of 2.5°; Supplementary Fig. 2). Calculation time for the entire docking screen was 11,740 CPU hours, elapsed wall time ~14 hours.</p><!><p>As noted, following docking, the top 500 (kinases) −1000 (AmpC) molecules were manually inspected for exclusion criteria based on considerations that are orthogonal to the docking scoring function such as structural diversity, commercial availability and or synthetic accessibility, perceived ability to derivatize and improve hits, correct representation of the molecule, and internal strain (ligand internal energy is not part of the scoring function). We remove redundant, highly similar molecules and compounds with potentially reactive or unstable functional groups. Additionally, for AmpC we mostly selected poses in which the boronic acid was predicted to occupy the oxyanion hole. For the kinase inhibitors we only selected poses predicted to form at least one hydrogen bond with the hinge. For JAK3 inhibitors we usually restricted the Cβcysteine-Sγcysteine-Cβacrylamide-Cαacrylamide dihedral angle to be >90° to avoid internal strain.</p><!><p>RMSD values were calculated using the Hungarian matching algorithm as implemented in DOCK6 51. For comparison purposes, RMSDs for the β-lactam benchmark were calculated using software generously provided by Xuchang Ouyang 16. For AmpC, receptors were superimposed based on chain A before calculations. For RSK2, receptors were superimposed based on C436, M496, and C560.</p><!><p>AmpC Enzyme inhibition was measured from initial rates using curve fitting in the native Agilent software. Compounds were initially dissolved in DMSO at 100 mM and individually diluted from such stocks. AmpC activity on CENTA (Km =15 μM) was monitored by the change of absorbance at 405 nm 24. AmpC was expressed and purified as described 24, CENTA was purchased from Tydock Pharma (Modena, Italy). IC50 values were obtained by fitting percent inhibition to a sigmoidal dose-response equation using GraphPad Prism (GraphPad Software Inc., Supplementary Fig. 16). Ki values (average of two biological experiments) were determined using the Cheng-Prusoff equation assuming competitive inhibition. Reactions were performed at room temperature in 50 mM sodium-cacodylate, pH 6.5 in the presence of 0.01% Triton-X-100 in 1 mL cuvettes with 50–100 μM CENTA and initiated by addition of AmpC.</p><!><p>Compounds 3 (1 μM) and 7 (276 nM) were assayed for AmpC (1nM) inhibition with no incubation as described above or after 5 minutes incubation with 10nM AmpC. In the latter, the reaction was started by 10x dilution of AmpC and the incubated inhibitor into a reaction buffer with substrate. Reversibility experiments were performed in 50 mM potassium phosphate buffer, pH 7.0</p><!><p>Compounds 2, 3, 5, and 7 were tested against porcine pancreas elastase (Sigma E-0258), porcine pancreas trypsin (Sigma T-0134), and bovine pancreas α-chymotrypsin (Sigma C-7762). The following substrates were used (Bachem Biosciences): Suc-Ala-Ala-Pro-Ala-pNA (L1775), Suc-Ala-Ala-Pro-Arg-pNA (L1720), and Suc-Ala-Ala-Pro-Phe-pNA (L1400) respectively. Serine protease activities were assayed at a concentration of 0.01 mg/mL enzyme in 50 mM Tris buffer, pH 7.0, with 0.01% v/v Triton X-100. The reaction was initiated by the addition of 200 μM substrate and monitored at 405 nm. IC50 values were calculated from single inhibitor concentration measurements (usually 100/1000 μM), and Ki values were estimated for each inhibitor and protein pair using reported Km values: 190 μM for elastase 52, 37 μM for trypsin 53 and 50 μM for α-chymotrypsin 54.</p><!><p>Each reaction contained the substrate Suc-LLVY-AMC (R&D Systems) at a final concentration of 150 μM; assay buffer (25 mM HEPES pH 7.4, 100 mM KCl, 20 mM MgCl2, and 10% glycerol); 1% DMSO, 100 μM compound, or 5 nM PA26; and 1 nM of the yeast 20S proteasome (except for the control reaction). Following the addition of the proteasome, fluorescence was read every 20 seconds, for 20 minutes. Experiments were conducted in duplicates at 23°C. Activity rates were calculated based on the last 10 minutes of a reaction. The yeast 20S proteasome and 26S proteasome activator complex were generous gifts from the laboratory of Philip Coffino. 100 μl reactions were performed in 96-well plates using a SpectraMax M5 Microplate Reader.</p><!><p>Susceptibility testing was performed and interpreted following the guidelines of the Clinical and Laboratory Standards Institute 25. The compounds were dissolved in DMSO, and dilutions were made into Muller Hinton medium, keeping DMSO < 5%. Inhibitors were tested for synergy with the third-generation β-lactam cefotaxime against clinical bacteria. The ratio of β-lactam to inhibitor was 1:4. Each value reported reflects the average of three independent experiments. The bacteria exhibited high levels of resistance to cefotaxime because of the expression of class C β-lactamases (AmpC) or class A extended-spectrum β-lactamases (ESBLs). Three Escherichia coli strains and one strain each of Citrobacter freundii, Enterobacter aerogenes and Enterobacter cloacae showed an AmpC-overproduction phenotype. Finally, two of the E. coli strains produced the plasmid-mediated class A ESBLs, TEM-3 and CTX-M-14 (Table 1).</p><!><p>Co-crystals of AmpC in complex with the inhibitors MAPB, 3, 7, and 14 (PDB codes: 4LV0, 4LV1, 4LV2, 4LV3) were grown by the hanging drop vapor diffusion method equilibrated over 1.7 M potassium phosphate buffer (pH 8.6–8.8). Protein (4 mg/mL) was mixed with 1 mM inhibitors and incubated for 20–30 minutes. Drops were set up by mixing 2 μL of protein-inhibitor with 2 μL of well solution; 1 μL of microseeding solution was added to promote crystal growth. Crystals appeared after 2–7 days of equilibration at 20°C. Before data collection, crystals were immersed in a cryoprotectant solution composed of 25% sucrose, 1.7 M potassium phosphate, pH 8.7, for about 30s and were flash-cooled in liquid nitrogen. The cryoprotectant solution also contained the respective inhibitor at 1 mM concentration.</p><p>Diffraction was measured at beamline 8.3.1 of the Advance Light Source (ALS, Lawrence Berkeley National Laboratory). Reflections were indexed, integrated and scaled using the XDS package 55 in the space group C2 with two molecules in the asymmetric unit. Structure refinement was carried out using Phenix 56. Coot 57 was used for model building and eLBOW 58 was used to generate coordinates and ligand restraints; The Ser64 Oxygen-Boron bond length was set to 1.45 Å.</p><p>The initial phasing model was based on apo-AmpC (PDB: 1KE4), with water molecules and ions removed. B-factors were refined isotropically and the protein was subjected to TLS refinement; three TLS groups were determined for each chain using the TLSMD server 59.</p><p>Similarly, co-crystals of RSK2 T493M in complex with compound 24 (PDB: 4M8T) were grown by hanging drop diffusion as described previously 27, and diffraction data were collected at the ALS on beamline 8.2.2. Diffraction images were indexed and integrated using XDS. Molecular replacement was performed using apo-RSK2 CTD as a starting model (PDB: 2QR8) using Phaser, and TLS refinement was carried out using Phenix. See Supplementary Table 10 for crystallographic statistics.</p><!><p>Wild-type and T493M RSK2 kinase activity were assayed as reported previously 26. Briefly, ERK2-activated RSK2 CTD (5 nM) was incubated with varying concentrations of each inhibitor for 30 minutes in the presence of 100 μM ATP and 10 mM GSH. Each reaction was initiated by the addition of 167 μM substrate peptide (RRQLFRGFSFVAK) and 0.3μCi/μL γ-32P-ATP in a final volume of 25 μL for an additional 30 minutes. Reactions were spotted on phosphocellulose membranes, washed once with 10% AcOH, twice with 0.1% H3PO4, and once with MeOH prior to drying. Blots were exposed to a phosphor storage plate (Amersham Biosciences), imaged with a Typhoon scanner (GE Healthcare), and quantitated using the SPOT program 60. IC50 values were calculated using a sigmoidal dose response fitting in the Prism program (Graphpad) and are reported as the average of two biological experiments. Full IC50 curves are presented in Supplementary Fig. 17.</p><!><p>JAK3 kinase activity was assayed using recombinant JAK3 (Invitrogen, catalog # PV5774). JAK3 (3.1 nM) was incubated with varying concentrations of each inhibitor for 30 minutes in the presence of 11.5 μM ATP. Each reaction was initiated by the addition of 17.9 μM substrate peptide and 0.3 μCi/μL γ-32P-ATP in a final volume of 25 μL for an additional 60 minutes. Reactions were spotted on phosphocellulose membranes, washed once with 10% AcOH, twice with 0.1% H3PO4, and once with MeOH prior to drying. Blots were exposed to a phosphor storage plate (Amersham Biosciences), imaged with a Typhoon scanner (GE Healthcare), and quantitated using the SPOT program60.</p><!><p>Kinase selectivity dose-response experiments for 31 and 33 were performed by Nanosyn (Santa-Clara, CA). Test compounds were diluted in 100% DMSO using 3-fold dilution steps. Final compound concentration in assay ranged from 10μM to 0.056 nM. Compounds were tested in a single well for each dilution, and the final concentration of DMSO in all assays was kept at 1%. All assays were performed with a substrate concentration of 1 μM and Km ATP concentration. Enzyme concentrations ranged 0.1–2 nM, and incubation times were 2–4 hours. See Supplementary Table 11 for exact concentrations. IC50 values were calculated using a sigmoidal dose response fitting in the Prism program (Graphpad) and are reported as the average of two biological experiments.</p><!><p>Compound 24 (5 μM) and WT RSK2 CTD (100 nM) were incubated in the presence of 10 mM GSH and 100 μM ATP for 1 hour in RSK kinase buffer (20 mM HEPES buffer, pH 8.0 with 10 mM MgCl2, 0.2 mg/mL BSA and 2 mM Tris(carboxyethyl)phosphine (TCEP)). The solution was then diluted 20-fold in duplicate into the RSK kinase buffer or RSK kinase buffer containing compound 24 (final concentration 5 μM). After 1 hour, kinase activity was then assayed as previously described.</p><p>Compound 31 (2 μM) and JAK3 (62.6 nM) were incubated in Invitrogen's recommended JAK3 kinase buffer (50 mM HEPES, pH 7.5, 0.01% BRIJ-35, 10 mM MgCl2 and 1mM EGTA) in the presence of 11.5 μM ATP, 2.5 mM DTT and 0.2 mg/ml BSA for 1 hour. The solution was then diluted 20-fold in duplicate into the JAK3 kinase buffer or JAK3 kinase buffer containing compound 31 (final concentration 2 μM) in the presence of 100 μM ATP, 2.5 mM DTT and 0.2 mg/ml BSA for 1 hour. After an additional hour, kinase activity was assayed as previously described.</p><!><p>Confluent COS-7 cells were transfected overnight with HA-tagged full-length RSK2, MSK1, and MSK1 C458V, as reported previously 27. Transfected cells were seeded into 6-well plates at 500,000 cells/well in 2 mL DMEM supplemented with 10% FBS (Axenia), 100 units/mL penicillin and 100 μg/mL streptomycin (Gibco) and allowed to adhere for 4 hours. Cells were serum starved for 18 hours, then treated with inhibitor or DMSO for an additional 2 hours. Cells were stimulated with PMA (100 ng/mL) for 30 minutes, washed with cold PBS, and frozen. Cells were thawed into 60 μL of lysis buffer (50 mM HEPES, pH 7.4, 150 mM NaCl, 0.1% Triton X100 supplemented with Roche phosphatase and protease inhibitor cocktails). Lysates were clarified by centrifugation at 14,000 rpm, normalized by Bradford assay, denatured in SDS and separated by 7.5% acrylamide SDS-PAGE. Gels were transferred to nitrocellulose, blocked with Odyssey LiCOR blocking buffer for one hour, and probed with 1:1000 HA, 1:1000 pS380 RSK2 or 1:1000 pS376 MSK1 antibody dilutions. After thorough washing with TBST, blots were incubated with 1:10,000 dilutions of fluorescent secondary antibodies (Odyssey) for one hour, washed with TBST, and scanned on an Odyssey LiCOR instrument. Raw gel images can be found in Supplementary Fig. 15. phospho-S380 RSK2 (Cat #9335S), phospho-S376 MSK1 (Cat #9591S) antibodies were purchased from Cell Signaling. The HA antibody was 12CA5 from Roche, Cat #11 666 606 001. COS-7 cells were obtained from the ATCC.</p><!><p>Compounds were sourced from the following vendors: Combi-Blocks: 1–11, 13–15,, Matrix Scientific: 12, Alfa-Aesar: 16, TCI: 17, PepTech: 18. These compounds were sourced at 95% or greater purity as described by the vendors. Compounds 19–42 were synthesized in-house, see Supplementary Note 2 for synthetic chemistry details and characterization.</p><!><p>All crystal structures reported here were deposited in the Protein Data Bank (PDB) with accession codes 4LV0, 4LV1, 4LV2, 4LV3 and 4M8T.</p>