Combining Mutational Signatures, Clonal Fitness, and Drug Affinity to Define Drug-Specific Resistance Mutations in Cancer

Summary The emergence of mutations that confer resistance to molecularly targeted therapeutics is dependent upon the effect of each mutation on drug affinity for the target protein, the clonal fitness of cells harboring the mutation, and the probability that each variant can be generated by DNA codon base mutation. We present a computational workflow that combines these three factors to identify mutations likely to arise upon drug treatment in a particular tumor type. The Osprey-based workflow is validated using a comprehensive dataset of ERK2 mutations and is applied to small-molecule drugs and/or therapeutic antibodies targeting KIT, EGFR, Abl, and ALK. We identify major clinically observed drug-resistant mutations for drug-target pairs and highlight the potential to prospectively identify probable drug resistance mutations.


INTRODUCTION
Although targeted cancer therapies, for example, against kinases, hormone receptors, or hormone-synthesizing enzymes, have shown clinical success, many patients develop resistance to treatment and subsequently relapse. Second-and third-generation drugs are being developed to target these resistant mutants; however, there is a significant time span between the detection of clinically validated resistance mutations and the availability of suitably targeted treatment options. Early identification of drug-specific mutations is therefore critical and the aim of this study.
Several mechanisms underlying resistance to targeted drugs have been described (Holohan et al., 2013), including mutations directly affecting the drug target. Such mutations may, for example, increase affinity for the endogenous co-factor ATP, thereby decreasing the relative affinity of an ATP-competitive drug. Mutations within the binding site may also alter drugprotein interactions and directly interfere with drug binding (Barouch-Bentov and Sauer, 2011).
Several computational studies have investigated the impact of protein mutations on drug efficacy (Kamasani et al., 2017;Mar-tínez-Jimé nez et al., 2017;Paul et al., 2016;Reeve et al., 2015), mainly in an antiviral or antibacterial context (Frey et al., 2010;Hosseini et al., 2016;Fowler et al., 2018). However, rarely do these studies prospectively identify drug resistance mutations (Frey et al., 2010;Paul et al., 2016;Reeve et al., 2015). In predicting mutations that render Staphylococcus aureus resistant to an antifolate antibiotic, Reeve et al. (2015) evaluated the likely effect of possible mutations on both binding of the inhibitor and on binding of the endogenous ligand -an important aspect since any mutation that significantly abrogates the native activity of the wild-type (WT) protein is unlikely to survive selective evolutionary pressure (Gil and Rodriguez, 2016;Sprouffske et al., 2012;Pandurangan et al., 2017). However, Reeve et al. do not consider the likelihood of whether each mutation can be formed in bacteria.
In cancer, the mutation landscape of a tumor can be characterized by the mutational signatures operating in a particular cancer type (Alexandrov et al., 2013). These signatures describe the probability of a specific base exchange within a defined trinucleotide context. Some of these signatures have been associated with known mutagenic processes, such as UV irradiation or aging, while the mechanism of others still remains elusive (Alexandrov et al., 2013). These mutagenic processes can generate a single clone harboring the disease-causing ''driver mutation,'' which ultimately leads to the development of cancer (Greaves and Maley, 2012). In addition, non-transforming somatic mutations, so-called passenger mutations, are randomly created. While not oncogenic per se, passenger mutations can provide the substrate for an evolutionary advantage throughout cancer progression, for example, under the selective pressure of a targeted molecular therapy, leading to drug resistance. Known drug resistance mutations have not only been detected in treatment-naive patients (Inukai et al., 2006;Roche-Lestienne et al., 2002), but also in healthy individuals (Gurden et al., 2015). This suggests that small pools of viable treatment-resistant clones can pre-exist in patients and that drug treatment puts a selection pressure on a heterogeneous cancer cell population that selects for resistant sub-clones.
Each drug interacts with its biological target in a unique way, and each protein target mutation will differentially affect diverse classes of drugs. As a consequence, each compound can be expected to exhibit a unique resistance mutation profile. Three factors contribute to the probability and functional impact of a residue change: (1) the probability that the protein mutation can be generated from a DNA mutational signature (signature-driven probability), (2) whether the mutation maintains protein function and clones harboring the mutation are still viable (fitness), and (3) whether the mutation confers lower drug affinity with respect to the endogenous ligand for the target protein (affinity). Martínez-Jimé nez et al. (2017) recently reported a workflow classi-fying potential drug resistance mutations based on Random Forest models and mutation signatures. However, the effect of mutations on the fitness of the clone was not taken into account. In addition, only single-point mutations (SPMs) were considered, despite the notable detection of double-point mutations (DPMs) in cancer patients (Table S1).
We report an in silico cascade that sequentially evaluates the probability of generating any mutant within 5 Å of a bound ligand, the clonal fitness of each mutation, and the effect of each mutation on drug affinity in order to systematically and objectively prioritize mutations that are highly likely to arise under drug treatment. Importantly, our workflow classifies the impact of a mutation on drug affinity relative to endogenous ligand and does not rely upon accurate calculation of binding free energies; it also ranks mutations according to their likelihood of being generated in particular cancer type. The workflow ( Figure 1 and described in detail below) is validated on a comprehensive benchmark dataset that describes the effect of nearly all possible extracellular signal-regulated kinase 2 (ERK2) missense mutations on sensitivity to the ERK2 inhibitor SCH772984 (Brenan et al., 2016). In addition, we apply the workflow to four well-established cancer targets and evaluate at least two US Food and Drug Administration-approved drugs (small molecules or biologics) per target.

Affinity
In the first step, the software package Osprey Gainza et al., 2013) was applied to evaluate the potential impact of protein active site mutations on both drug and endogenous ligand binding using protein crystallographic data.
Osprey uses ensembles of minimized side chain rotamers to calculate a K* score Gainza et al., 2013), which estimates ligand binding to the evaluated protein construct. Importantly, by employing ensembles of rotamers, Osprey takes the dynamic nature of protein amino acid side chains into account. K* scores were computed for the WT protein and all possible amino acid exchanges for both the drug and the endogenous ligand-protein complex with the exception of Pro, which cannot be processed by Osprey. We identify mutants that reduce affinity for the drug in comparison with the endogenous ligand and are therefore more likely to drive a viable resistant clone. Importantly, calculation of absolute drug affinity for WT and mutant proteins is not required. Please refer to the STAR Methods for additional details.

Fitness
Clones harboring a mutant oncoprotein that is unable to bind its endogenous ligand lose their function and are unlikely to persist. All such deleterious mutants with a K* score of 0.0 for the endogenous ligand were discarded.

Exclusion of Triple Point Mutations
Base codons for each WT residue were retrieved from the COSMIC database . Amino acid exchanges that would require mutation of all three codon bases were discarded. We hypothesized that three contiguous coding Potential mutations are evaluated based on their predicted effect on the affinity of both the drug and endogenous ligand (orange), the fitness of the resultant clone (blue), and the requirement for triple-point mutations to generate a mutant (lime green). Resistance hotspots are identified within the remaining set of resistant mutants; these resistance hotspots are protein residues where multiple amino acid changes are predicted to lead to resistance and which therefore have a high likelihood of functional relevance. Resistant mutations at these hotspots are prioritized based on the probability that they will be generated according to the known DNA mutational signatures operating in a particular cancer type. base mutations are highly unlikely and, to the best of our knowledge, no such patient case has been reported. While the majority of drug resistance mutations arise from SPMs, we also included DPMs (so-called tandem mutations), which represent 5% (7 of 132) of unique drug resistance mutations with defined genetic alterations in the COSMIC database (access date 5 January 2018) and occur with similar frequency to insertions (Table S1).
Resistance Hotspots and Signature-Driven Probability Functional protein mutants were further triaged according to the number of predicted mutations at a specific residue position. This approach identifies resistance hotspots, where multiple different viable mutants are predicted to cause resistance to the drug under study. While comparing the prevalence of mutations at hotspots in a particular drug-target pair can be used for ranking positions, this method cannot be applied to crosscompare the prevalence of mutations at hotspots between drug-target pairs to infer global inter-target susceptibility to resistance. The top three ranked resistance hotspots were then analyzed in more detail in order to further prioritize amino acid changes at these positions that have a high probability (relP) of being generated in particular tumor types according to the DNA mutational signatures reported by Alexandrov et al. (2013). Each possible resistance hotspot mutation was ranked according to its derived relP and the top-ranked amino acid changes were considered most relevant. Please refer to the STAR Methods for additional details.

ERK2: A Benchmarking Study
In 2016, Brenan and colleagues systematically studied the effect of a wide range of ERK2 genetic alterations on the cellular response to the ERK2 inhibitor SCH772984 (Brenan et al., 2016). They generated nearly every possible missense mutation using multiplexed site-directed mutagenesis. Pooled vectors were used to express ERK2 mutants in A375 human malignant melanoma cells; after drug exposure, the enrichment of resistant variants relative to their abundance in the initial cell population was determined by massively parallel sequencing.
This comprehensive dataset of possible resistant mutants to the ERK2 inhibitor SCH772984 provides a benchmark to evaluate the performance of our workflow. However, as a consequence of the aggressive mutant-generation protocol, this dataset includes protein mutants that require exchange of all three coding bases. Such mutants are excluded in our workflow (see above) and were also removed from the Brenan dataset.
Brenan and colleagues reported 46 experimentally observed resistance mutations within 5 Å of the bound ligand that were not mutations to Pro or triple-point mutations. For SCH772984, 559 mutants covering 31 residues were evaluated by our workflow. The number of true-positive (TP) (mutants predicted to confer resistance that were confirmed by experimental testing), true-negative (mutants predicted to be sensitive that did not confer resistance in experimental testing), false-positive (FP) (mutants predicted to confer resistance that were found to be sensitive in experimental testing), and false-negative (mutants predicted to be sensitive, but conferred resistance in the experimental testing) predictions was calculated ( Figure 2A).
Our workflow correctly classifies 84.6% of mutants and identifies almost 40% of the experimentally observed resistance mu-tations reported by Brenan et al. When considering only residue positions, 80% of experimentally observed mutated residues are discovered. The majority of experimentally evaluated ERK2 mutants remained sensitive to SCH772984 treatment and are correctly identified as sensitive by our approach, contributing to the high overall proportion (84.6%) of correct predictions. The low TP rate for SCH772984 (3.2% of mutants predicted to confer resistance were confirmed by experimental testing) is consistent with the experimental observation that only a small proportion of possible mutants (8.2%) cause drug resistance. The prioritized hotspot residues and their corresponding reported resistance mutations are provided in Table 1. For a detailed list of all predicted mutations please refer to Table S2.

Clinical Case Studies
We studied four protein kinases (KIT, epidermal growth factor receptor [EGFR], breakpoint cluster region-Abelson [Bcr-Abl] kinase, and anaplastic lymphoma kinase [ALK]). In each case protein-ligand crystal structures of approved drugs exemplify multiple generations of compound design; furthermore, as inhibitors of these targets provide the current standard of care for multiple cancer types, clinically observed resistance mutations are documented for most compounds and provide the potential for clinical validation of our computational workflow.

KIT
The first-generation inhibitor imatinib is resistant to the KIT gatekeeper mutation T670I (Antonescu et al., 2005;Tamborini et al., 2004). We investigated 648 possible mutations, all of which lie within 5 Å of either imatinib or ADP; 68 mutations remained after filtering for reduced drug binding, fitness of the clones, and triple-point mutants. For residue T670, eight mutations are predicted to cause resistance to imatinib, making it the first-ranked resistance hotspot; within this set, gatekeeper mutation T670I has the highest relP. Furthermore, the known resistance mutation T670E (Wardelmann et al., 2005) is ranked sixth by relP. Positions C809 and V668 are identified as second-and thirdranked resistance hotspots, respectively. C809R, the mutation with the highest relP, has been reported in a patient diagnosed with myelodysplastic syndrome-derived leukemia (Lorenzo et al., 2006), highlighting its biological relevance; however, similar to all predicted V668 mutations, this residue has not yet been associated with imatinib resistance.
The second-generation inhibitor sunitinib is reported to overcome resistance to the gatekeeper mutation T670I and, consistent with these reports, T670I/E mutations are not predicted to elicit sunitinib resistance by our workflow. Indeed, none of our predicted resistance mutations to sunitinib has yet been associated with clinically observed resistance. Only one resistance mutation has so far been reported for sunitinib in the COSMIC database , access date 23 November 2017); this mutation is not within the ATP-binding site.
Ponatinib has been shown to bind KIT, but is not approved for KIT-associated cancer types and clinical data are therefore lacking. Our method predicts resistance to ponatinib due to less-favorable interactions with the gatekeeper T607I mutation; however, in vitro data suggest that this mutant is sensitive to ponatinib treatment (Garner et al., 2014). Prioritized KIT resistance hotspot residues and their corresponding reported clinical resistance mutations are provided in Table 2. A full list of predicted resistance mutations and their relP values is available in Table S3.

EGFR
Applying our workflow to the first-generation EGFR inhibitor erlotinib revealed that G796 is the resistance hotspot with the most predicted mutations. Consistent with this prediction, the G796R mutation with the second-ranked relP, has been reported to weaken the affinity for erlotinib in comparison with the WT enzyme (Avizienyte et al., 2008). T790 is the second-ranked resistance hotspot with gatekeeper mutation T790M (Pao et al., 2005) having the second highest relP. L718, L788, and T854 share the third-ranked hotspot, with three predicted (B) Predicted resistance hotspots that are consistent with clinically observed resistant mutants for the representative case study EGFR (gray) and osimertinib (orange), PDB: 4ZAU (Yosaatmadja et al., 2015). Residues within 5 Å of the ligand are depicted in gray, predicted and clinically observed resistance hotspots are highlighted as crimson sticks and labeled, predicted hotspot residues that have not yet been observed in the clinic are shown as pink sticks; figure created with PyMOL (PyMOL Molecular Graphics System, Version 1.7, Schrö dinger, LLC). (C) Contribution of filtering steps to the identification of resistance mutations. The majority of mutants were discarded because they did not decrease drug affinity in comparison with binding of the endogenous ligand (orange). Mutations were further removed because of abrogated clonal fitness (blue) or because they required triple codon changes to be formed (green). The remaining pool of mutations (crimson) is predicted likely to confer resistance to drug treatment. mutations each. L788F has been detected in lung adenocarcinoma patient samples highlighting its clinical relevance (Liu et al., 2013); however, similar to all other third-rank mutations, L788F has not yet been associated with clinical resistance to erlotinib. Gefitinib, another first-generation EGFR inhibitor, elicits a similar predicted resistant mutation profile to erlotinib (Table 3).
The third-generation inhibitor osimertinib binds reversibly to EGFR prior to covalent bond formation (Yosaatmadja et al., 2015). Osimertinib was reported to overcome resistance to T790M and, consistent with this finding, the T790M mutation is not predicted to elicit resistance to osimertinib according to our protocol. Residues G796 (14 mutations), C797 (4 mutations), L718, G719, V726, and A743 (3 mutations each) are the topranked hotspot residues predicted to confer resistance. During the preparation of this manuscript, osimertinib resistance mutations affecting C797 (Ou et al., 2017b), G796 (Ou et al., 2017b;Zheng et al., 2017), and L718 (Ou et al., 2017a) (Figure 2B) were reported in clinical studies, building confidence in the prospective utility of our approach.
Further to small-molecule EGFR inhibitors, anti-EGFR monoclonal antibodies cetuximab and panitumumab are approved for the treatment of metastatic colorectal cancer. In contrast to small-molecule ATP-competitive kinase inhibitors, cetuximab and panitumumab interact with the extracellular domain of EGFR.
We evaluated the effect of extracellular domain mutations on antibody affinity and affinity of the endogenous ligands EGF and transforming growth factor a (TGF-a)--interactions that have been characterized by protein structural data. Both EGF and TGF-a are peptidic macromolecules with larger interaction surfaces compared with small molecules; thus evaluation of a higher number of possible mutations was necessary to ensure coverage of the extensive protein-protein interaction interface (1,153 and 1,207 mutations for cetuximab and panitumumab, respectively).
For both cetuximab and panitumumab, resistance mutation G441R (residue G465R in the mature protein sequence including the signaling peptide [mps]) has been detected in the clinic (Braig et al., 2015). G441 is the second-ranked resistance hotspot, with mutation G441R ranked first according to relP for both antibodies (Table 3). Confirmed clinical resistance mutations to cetuximab have also been identified at residues S468 (S468R [S492R in the mps]) (Montagut et al., 2012) and I467 (Arena et al., 2015) (I467M [I491M in the mps]); the fourth-and fifthranked resistance hotspots in our workflow (Table S4).
Our top-ranked predicted hotspot residues, S418 and G471 have not yet been associated with clinical resistance to cetuximab or panitumumab; however, seven endogenous ligands can activate EGFR (Schneider and Wolf, 2009) of which only EGF and TGF-a interactions are characterized by protein structural data. Residues S418 and G471, while not predicted to affect the interaction of EGFR with endogenous ligands EGF and TGF-a, may be important for interaction with one or more of the remaining five endogenous ligands such that mutation at these residues reduces clonal fitness and persistence. For example, we predict that several mutations affecting EGFR residues S440 and V417 tolerate TGF-a binding, but abrogate EGF binding, and are therefore filtered out in our workflow. Given the lack of structural data for the majority of endogenous EGFR ligands, and the number of potential mutations, the identification of clinically relevant resistance mutations is encouraging. All prioritized EGFR hotspot residues and corresponding reported clinical resistance mutations are summarized in Table 3. A list of all predicted mutations and relPs is provided in Table S4.

Abl
The hotspot residue with the highest number of predicted mutants for the Abl inhibitor imatinib is A380, for which clinically observed occurrences have yet to be reported. The secondrank predicted hotspot is shared by residues V256, A269, and T315. For V256 and T315, resistance mutations V256L Tiribelli et al., 2013), gatekeeper mutation T315I (Gorre et al., 2001), and T315V (Redaelli et al., 2012) have all been identified in the clinic ( Table 4). The contribution of mutational signatures to the landscape of mutations has not yet been reported for chronic myeloid leukemia (CML), and therefore the relP for hotspot mutations predicted for Abl residues could not be calculated. We obtained similar and consistent predictions for nilotinib, dasatinib, and bosutinib (Table 4). Additional resistance mutations affecting L248 are predicted for dasatinib consistent with in vitro data (Redaelli et al., 2012). While these mutations have been observed in the clinic (Redaelli et al., 2012), they have not yet been reported as resistant mutants for patients treated with dasatinib. Similarly, the effect of the predicted V256L resistant mutation on nilotinib activity has not been reported.
Axitinib has been reported to overcome resistance to the T315I gatekeeper mutation of Abl (Pemovska et al., 2015). Although not ranked among the top three resistance hotspots, axitinib is predicted to be resistant to T315I. Crystallographic studies of the T315I mutation in complex with axitinib revealed that this mutation causes large conformational changes of Abl compared with WT protein (Pemovska et al., 2015), which may not be adequately captured by our method. The L248R mutation, reported to confer resistance to axatinib in vitro (Pemovska et al., 2015), is prioritized by our workflow; however, clinical resistance data for axitinib targeting Abl is not yet available. A summary of the prioritized Abl hotspot residues and their corresponding reported clinical resistance mutations are provided in Table 4. A full list of predicted resistance mutations is provided in Table S5.

ALK
The hotspot residue with the highest frequency of predicted mutations for the first-generation ALK inhibitor crizotinib is G1269 (13 mutants). The clinically observed resistance mutation G1269A (Doebele et al., 2012;Gainor et al., 2016) and G1269S and G1269C mutations (Zhang et al., 2011), which confer resistance in vitro, rank at positions four, seven, and eight, respectively, based upon their relP values. For the second-ranked resistance hotspot G1202, 12 mutants were predicted and the known resistance mutation G1202R (Gainor et al., 2016;Katayama et al., 2012) has the second highest relP. Hotspots I1122, G1201, and D1203 rank third equal. While no G1201 mutations are associated with resistance, the D1203N mutation confers resistance in vitro (Gainor et al., 2016;Heuckmann et al., 2011) and has been associated with resistance in patients (Zhang et al., 2016). I1122V has been detected in a resistance screen against the second-generation inhibitor brigatinib, and has also T670I was predicted to confer resistance to ponatinib; however, T670I is reported to be sensitive to ponatinib (Garner et al., 2014). f NA, not applicable, tied resistance hotspot at rank 2. The gatekeeper mutation is underlined.
been confirmed to confer resistance to crizotinib (Ceccon et al., 2015). We generated similar results for the second-generation inhibitor ceritinib (Table 5), although, in this case, the G1269A mutation was not predicted to cause resistance, consistent with in vitro and clinical data (Gainor et al., 2016;Shaw et al., 2015). Entrectinib and lorlatinib are investigational ALK inhibitors currently in clinical trials. Increased half maximal inhibitory concentration values (decreased binding) in comparison with the WT protein have been demonstrated in in vitro experiments for several of the mutants prioritized by our workflow. In particular, G1269A  and G1202R  for entrectinib and G1269A (Gainor et al., 2016;Johnson et al., 2014;Shaw et al., 2015) and D1203N (Gainor et al., 2016) for lorlatinib. In addition, the G1123S and G1123D mutants were identified in an in vitro resistance screen for the ALK inhibitor TAE684 (Heuckmann et al., 2011) and in a patient resistant to ceritinib (Toyokawa et al., 2015). However, no data are currently available on entrectinib. The prioritized hotspot residues for ALK and corresponding reported clinical resistance mutations are provided in Table 5. A comprehensive list of all predicted resistance mutations is provided in Table S6.

DISCUSSION
We present a computational workflow to identify clinically relevant drug resistance mutations to targeted cancer therapies, both small molecule and biological. The workflow consists of consecutive filtering steps addressing three factors, which determine whether a mutation is likely to confer resistance in the clinic. Of these, the activity cutoff for protein-drug affinity compared with affinity for the endogenous ligand proved the most important and excluded 75.8% ± 1.1% (n = 18 [all 18 target-interaction partner case studies] ±SEM) of all potential mutants for each drug-target combination ( Figure 2C). This filtering step contains a fitness component; all mutants with predicted decreased affinity for the endogenous ligand (i.e., less fit clones), but with equivalent or increased affinity for the drug (i.e., more sensitive to drug treatment) are not progressed. In further considering the potential for disrupted binding of an endogenous ligand to mutated proteins, 8.6% ± 0.7% (n = 18, ±SEM) of mutants are predicted to completely abrogate binding of the endogenous ligand and were discarded. Removing triple-point mutations excluded a further 3.7% ± 0.2% (n = 18, ±SEM) of possible mutants. Interestingly, the majority of remaining mutations predicted to cause resistance (68.1% ± 1.7% [n = 18, ±SEM]) arise from DPM consistent with a similar proportion of DPMs observed in the Brenan dataset (Brenan et al., 2016). While DPMs commonly exhibit a low relP, they may become relevant when the predominant resistance clone harboring an SPM has been eradicated by cancer therapy. Taken together, 11.8% ± 1.0% (n = 18, ±SEM) of all evaluated mutations were considered to cause resistance ( Figure 2C).
The pool of potential resistant mutations was further analyzed to identify resistance hotspots where multiple different viable mutants are predicted to cause resistance to the drug under study. Mutants at these hotspots were prioritized based on their relP, which quantifies the relative probability of each specific amino acid mutation in the context of a defined cancer type. Applying this workflow, we correctly classify 84.5% of mutations  (Pao et al., 2005) L718 ( Table S4. a The number of mutants predicted to confer resistance (no. of pred mut) from the initial pool of possible mutants within 5 Å of the ligands (no. of all mut). b Resistance hotspots are identified and ranked according to the number of viable mutants (no. of mut) predicted for a residue. c The relP was calculated for all resistance hotspot mutations. Clinically observed resistance mutations and their rank according to relP (rank relP) are highlighted for each resistance hotspot. d NR, not reported--none of the predicted mutations are reported to confer resistance against the drug. The gatekeeper mutation is underlined.  (7) A269 (7) T315 (7) V256L Tiribelli et al., 2013) T315I (Gorre et al., 2001) NA e -Nilotinib 69 (540) V256 (8) A380 (8) NR NA -Y253 (6) A269 (6) T315 (6) G321 (6) T315I Redaelli et al., 2012;Weisberg et al., 2005) Dasatinib 58 (414) A380 (12) NR L248 (7) A269 (7) NR NA -Bosutinib 58 (414) V299 (7) T315 (7) G321 (7) V299L ( Table S5. a The number of mutants predicted to confer resistance (no. of pred mut) from the initial pool of possible mutants within 5 Å of the ligands (no. of all mut). b Resistance hotspots are identified and ranked according to the number of viable mutants (no. of mut) predicted for a residue. c Clinically observed resistance mutations are highlighted for each resistance hotspot. The relP could not be calculated as signatures for CML were not available. d NR, not reported--none of the predicted mutations were reported to confer resistance. In the case of axitinib, clinical resistance data on Abl are not yet available. e NA, not applicable, tied resistance hotspot at rank 1 or 2. The gatekeeper mutation is underlined.
in a comprehensive ERK2 mutation dataset and identify almost 40% of mutants conferring resistance to the ERK2 inhibitor SCH772984. Furthermore, we identify clinically observed drug resistance mutations within the top three predicted hotspot residues for first-generation compounds in all cases studied and for most second-and third-generation compounds where clinically observed resistance mutations have been reported. Importantly, in all cases studied, except ALK, the gatekeeper mutation is predicted among the top three resistance hotspots. In the ALK case study, two mutations that are also commonly observed in the clinic are highly ranked by our method. For the osimertinib-EGFR drug-target pair, we predicted resistance mutations before they were confirmed by clinical reports. Our workflow highlights mutations not yet observed in the clinic, or that may constitute FP predictions. The number of mutants progressed to further evaluation can be user defined. In the cases exemplified here, we prioritize three mutations with the highest relP for the top three resistance hotspots and significantly narrow the pool of potential resistance mutations for each drug from $350-1,200 possibilities to 9. While it may not be practical to further explore many hundreds of possible mutations, we suggest that a set of 9 prioritized mutants is more amenable to experimental testing within a drug discovery project or clinical setting. Thus, this workflow facilitates more focused monitoring of potential resistant mutations, as well as the design of next-generation compounds sensitive to likely resistant mutants.
Notably, our workflow evaluates the effect of a mutation on drug binding with respect to the WT protein. In some cases, exquisite drug affinity for the WT protein may mitigate a mutation-driven loss of potency. For example, our workflow correctly predicted the reduced affinity of lorlatinib for ALK G1269A; however, the exquisite potency of lorlatinib for WT ALK may render the relative loss of affinity for ALK mutants inconsequential (Gainor et al., 2016;Johnson et al., 2014;Shaw et al., 2015).
The workflow is dependent upon the availability and quality of both protein-ligand structural information and mutational signatures across diverse tumor types. For example, it is not clear which mutational signatures operate in CML, and, as a consequence, we could not determine the relP for mutants predicted to interfere with binding of Abl inhibitors in CML.
The approach presented here investigates the effects of protein mutation on drug binding mode as characterized by the input protein-ligand structure. How such mutations affect the overall conformation of the target protein and/or drug-target complex is not encompassed by our method. Extension of the approach by considering the effect of residue mutation on global protein conformation and by inclusion of structural models for highly homologous protein families would further expand the potential impact and is the aim of our future studies. For example, FP mutations that stabilize inactive protein conformations that are unlikely to persist could be excluded, and mutations that further stabilize an active protein conformation could be included. Furthermore, drug-resistant mutations distant from the binding site, including those which influence protein flexibility and conformation, which are beyond the scope of our current method, could also be evaluated.
This workflow includes three critical determinants of clinically relevant drug-resistant mutations. Importantly, the workflow can be further adapted to test more specific hypotheses. For example, Osprey parameters can be changed to include larger backbone movements and/or to include multiple ligand rotamers; furthermore, the number of prioritized mutants can also be user-defined.
In conclusion, we have developed and validated a computational method to prospectively identify clinically relevant drugresistance mutations. We suggest that this approach can have a significant impact on the design and development of targeted therapies by proactively signposting drug resistance hotspots. Prior knowledge of resistant mutants enables their timely detection in patients and the early development of effective treatment options against the resistant tumor cell population.

SIGNIFICANCE
Although molecularly targeted cancer therapies have shown great success in the clinic, drug resistance has emerged as the major challenge. Resistance mutations are commonly identified and characterized during clinical evaluation, often resulting in a reactive approach to tackling drug resistance. We report a computational method to prospectively identify drug resistance mutations during the design phase of potential therapeutics. This approach enables early signposting of likely resistance hotspots and supports more focused monitoring of potential emergent resistant clones as well as the timely development of alternative treatment options.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

AUTHOR CONTRIBUTIONS
T.K. conducted the experiments. J.B. and T.K. designed the experiments and wrote the paper.

DECLARATION OF INTERESTS
The authors declare no conflict of interest.

STAR+METHODS KEY RESOURCES TABLE CONTACT FOR REAGENT AND RESOURCE SHARING
Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact Julian Blagg (Julian.blagg@icr.ac.uk).

DATA AND SOFTWARE AVAILABILITY
Data are available upon request to the Lead Contact.

Structure Selection and Preparation
When available, crystal structures of wild type protein were used. Formation of the Lys-Glu salt bridge was the minimum requirement for kinase co-factor structures to be considered representative of an active conformation. Due to the paucity of high resolution ATP-bound kinase protein structures, we employed ADP or ATP analogue complexes for all kinase structures representative of endogenous co-factor binding. A list of employed PDB entries is provided in Table S7. All crystal structures were prepared using the Preparation wizard (Sastry et al., 2013) in Maestro version 9.8.016 (Schrö dinger Release, 2014). Bond orders were assigned and hydrogens were added. Zero-order bonds to metals and disulphide bonds were created and selenomethionines were converted to methionines. Missing side chains were added and all waters beyond 5 Å from heteroatoms, buffer compounds, and additional chains were deleted. Ionization states were generated using Epik. H-bonds were assigned, water conformations were sampled and structures minimized using default settings.
Residues within 12 Å of the ligand were selected using the selector tool in MOE 2015.1001 (Molecular Operating Environment (MOE), 2017) and only selected residues were included in the input structure. As described in the Osprey documentation, the protonation state of histidines was defined in the pdb input file, all HETATM identifiers were replaced by ATOM, and chain information was deleted.

Evaluation of Ligand Affinity
Osprey version 2.2beta Gainza et al., 2013), using the Amber94 force field, was employed to evaluate the impact of protein mutation on ligand binding. All non-peptide interaction partners were parameterized using Antechamber (AmberTools16 ) as described in the Osprey documentation and only the selected crystal structure protein conformation was used. Small molecule ligands were allowed to rotate and translate during the calculations.
Residues within 5 Å of the ligand were systematically mutated to all other possible amino acids, except if they were known to be crucial for catalytic activity; namely the AspPheGly motif, catalytic Asp, and Lys-Glu salt bridge in kinases. All three His protonation REAGENT