Introduction

Intrinsically disordered proteins (IDPs) are now recognized as important drug targets ^1-4. For structured protein targets, a crucial step is identifying the drug-binding pocket. Although an IDP can be locked into a specific conformation by a drug molecule in rare cases ⁵, the prevailing scenario is that the protein remains disordered upon drug binding. Consequently, the IDP-drug complex typically samples a vast conformational space, and the drug molecule only exhibits preferences, rather than exclusiveness, for interacting with subsets of residues. Such drug-interacting residues (DIRs), akin to binding pockets in structured proteins, are key to optimizing compounds and elucidating the mechanism of action.

NMR chemical shift perturbation (CSP) is the best experimental method for identifying DIRs and has been applied to many IDPs ^6-23. Aromatic residues are frequently among DIRs ^{7,10,12,17,18,20,22}. This is understandable as drug molecules typically are rich in aromatic rings (Figure 1 and Table S1), which can form π-π interactions with aromatic residues. As recently pointed out by Heller et al. ²³, CSPs of IDPs elicited by drug binding can be small enough to fall within the spectral resolution of NMR spectroscopy, therefore making it difficult to unequivocally identify DIRs. The small magnitude of CSPs may arise because the drug does not have a strong preference for interacting with any residues. Another scenario may be that drug binding induces a conformational shift such as secondary structure formation or even partial folding ^11,16,20, so CSPs spread from the directly interacting residues to the rest of the IDP, adding to the difficulty in identifying DIRs.

Figure 1.

All-atom molecular dynamics (MD) simulations have presented atomic details in many IDP-drug systems ^{17,18,20,22,24-29}. These simulation studies have highlighted the frequent engagement of drug molecules with aromatic residues, particularly in π-π interactions ^{17,18,20,22,26-29}. MD simulations also revealed the emergence of a compact sub-population of the N-terminal disordered region of the tumor repressor protein p53 upon binding the drug epigallocatechin-3-gallate (EGCG) ¹⁸ and reduced solvent exposure of hydrophobic residues in the 42-residue amyloid-β peptide (Aβ42) upon binding the drug 10074-G5 ²⁵. However, MD simulations of IDPs still suffer from the perennial issues of inaccurate force fields and insufficient sampling, which are exacerbated when drug molecules are included. For example, two replicate simulations may identify different DIRs ³⁰.

Here we present a sequence-based method, DIRseq, as a complement to NMR CSP and MD simulations. This method was motivated by our observation that DIRs seem to overlap with residues exhibiting elevated transverse relaxation rates (R₂) ^31,32, as exemplified by the C-terminal 20 residues of α-synuclein ²⁰. Elevated R₂ in IDPs is caused by either local inter-residue interactions or residual secondary structure ³³. Based on this understanding, we developed a sequence-based method, SeqDYN, to predict R₂ of IDPs ³¹. As suggested previously ³², the propensities of residues to form intramolecular interactions and therefore elevate R₂ should be similar to those for forming intermolecular interactions with drug molecules. DIRseq is therefore an adaptation of SeqDYN. DIRseq predictions match well with DIRs identified by NMR chemical shift perturbation and other methods, including residues L₂₂WK₂₄ and Q₅₂WFT₅₅ in p53 and C-terminal residues in α-synuclein. DIRseq is available as a web server at https://zhougroup-uic.github.io/DIRseq/ and has many applications.

Results

Retooling of SeqDYN into DIRseq

SeqDYN ³¹ predicts the R₂ value of residue n by accumulating a contributing factor f(i; n) from each residue i:

where N is the total number of residues in the IDP and ϒ is a uniform scale parameter. The contributing factor f(i; n) has an amplitude q(i) that is determined by the amino-acid type of residue i and attenuates with increasing sequence distance, s = |i − n|, from residue n:

The 20 q parameters (one for each amino-acid type) and b were trained on the measured R₂ values for 45 IDPs. The original q values (Figure S1a) show that aromatic (Trp, Tyr, His, and Phe), Arg, and long aliphatic (Ile and Leu) amino acids are interaction-prone and tend to elevate R₂. We transformed SeqDYN into DIRseq by implementing four changes.

First, we reassigned four q parameters. We lowered the q values of the three long aliphatic amino acids, Leu, Ile, and Met, to the q value of a short aliphatic amino acid, Val, because long aliphatic amino acids primarily participate in hydrophobic interactions, which may be less important for stabilizing the binding of a small molecule in sites largely exposed to water. At the same time, we raised the q value of Asp to be the same as that of Glu, to increase the role of Asp’s electrostatic interactions. The modified set of q parameters is displayed in Figure S1b.

Second, the original b value, 0.0316, corresponds to a correlation half-length, b^−1/2, of 5.6 residues. Given the small size of drug molecules, we increased the b value to 0.3, corresponding to a correlation half-length of 1.8 residues. That is, we expect a drug molecule to interact with 3 to 4 residues at a time. Third, the SeqDYN-predicted R₂ profile, capturing experimental observations ³⁴, falls off at both termini, because no residues beyond the termini are present to provide a contributing factor to R₂. However, whereas intra-IDP interactions experience such a terminal effect, IDP-drug interactions do not. To eliminate the terminal effect, we padded the original IDP sequence by a stretch of 12 Gln residues at each terminus. Gln was selected because its q value is at the middle of the 20 q values (Figure S1b).

Lastly, we converted high and low R₂ values into high and low DIR propensities, respectively, using a sigmoid function

where the midpoint R_2th and width R_2wd of the transition region are determined by the mean (m) and standard deviation (SD) of the R₂ values over the entire sequence. Specifically,

We chose s₁ to be 1.5 and s₂ to be 14.0. On the DIRseq web server, users can either keep these default values for b, s₁, and s₂, or enter values of their choice. The values of p range from 0 to 100; p = 50 when R₂ is at the “threshold” value R_2th. Residues with p ≥ 50 are predicted to be DIR, but an isolated residue with p ≥ 50 is discarded, as it can be reasonably expected that interactions with at least two consecutive (or nearby) residues are needed to generate sufficient drug binding stability.

Detailed assessment of DIRseq on four IDPs

DIRs in four IDPs or intrinsically disordered regions (IDRs) have been thoroughly characterized by NMR CSP. In Figure 1 (left images), we show these IDPs in a single conformation, with DIRseq propensities rendered in a color spectrum (from low in yellow to high in red); predicted DIRs (i.e., those with p ≥ 50) are shown with sidechains in stick. Drug molecules that bind to these IDPs are also displayed in Figure 1 (right images). We display NMR CSPs as blue bars and DIRseq propensities as a red curve in Figure 2. Unless otherwise noted, we use m + 1.5 SD as the threshold for identifying DIRs, for both CSP and DIRseq; this threshold is indicated by a horizontal dashed line in Figure 2. Experimentally identified DIRs are further indicated by cyan shading. Below we present a detailed comparison between CSP and DIRseq.

Figure 2.

The kinase inhibitory domain (residues 22-105) of the cell cycle regulator p27 harbors three aromatic-centered motifs that interact with a group of compounds represented by SJ403 (Figure 1a), as found by Iconaru et al. ¹⁰. These motifs, W₆₀N₆₁, E₇₅WQ₇₇, and Y₈₈Y₈₉ (Figure 2a, cyan shading), are correctly predicted by DIRseq (Figure 2a, portions of the red curve above the 50% threshold). Specifically, the predicted DIRs are R₅₈KWN₆₁, Y₇₄EWQ₇₇, and Y₈₈Y₈₉ (Figure 1a, left image). Moreover, the DIRseq propensity profile matches well the CSP profile along the sequence (Figure 2a). The two parameters show a very strong correlation, with a Pearson correlation coefficient (r) of 0.82 (Figure S2a).

Iconaru et al. also studied SJ403 binding to a related cell cycle regulator p21. Their CSP data identified the first two corresponding motifs, W₄₉N₅₀ and F₆₃A₆₄, but not the third as DIRs (Figure 2b, cyan shading). DIRseq correctly predicts both the first two motifs above the 50% threshold and the third motif below the threshold (Figure 2b, red curve; Figure 1b, left image). As noted by Iconaru et al., the p21 counterpart to the p27 third motif F₈₇YY₈₉ is L₇₆YL₇₈; the two Leu residues in p21 are no match for F₈₇ and Y₈₉ in p27. For the DIRseq method, Leu has a lower q value than Tyr and Phe (Figure S1b). CSP and DIRseq propensity show a strong correlation (r = 0.66; Figure S2b).

NMR CSP revealed two motifs, W₂₃K₂₄ and D₄₉IEQWFT₅₅, in the N-terminal region of p53 as interacting with EGCG ¹⁸ (Figure 2c, cyan shading). MD simulations supported these drug-binding sites. DIRseq again correctly predicts these two motifs (Figure 2c, red curve). Specifically, the predicted DIR motifs are L₂₂WK₂₄ and E₅₁QWFT₅₅ (Figure 1c, left image). Here also CSP and DIRseq propensity show a strong correlation (r = 0.67; Figure S2c).

Using the m + 1.5 SD threshold, CSPs of Robustelli et al. ²⁰ select residues A₁₂₄YEMPS₁₂₉ and Y₁₃₃QDYE₁₃₇ of α-synuclein as interacting with the drug Fasudil (Figure 2d, cyan shading). These two C-terminal motifs are also identified by DIRseq, with residues A₁₂₄YE₁₂₆ and Y₁₃₃QDYEP₁₃₈ above the 50% threshold (Figure 2d, red curve; Figure 1d, left image). There is a moderate correlation (r = 0.51; Figure S2d) between CSP and DIRseq propensity. Robustelli et al. further characterized the IDP-drug interactions by all-atom simulations, using both full-length α-synuclein and a C-terminal fragment (residues 121-140). Further simulations of the C-terminal fragment binding with additional compounds led to a compound known as Ligand 47 with a somewhat higher affinity than Fasudil. The CSP profiles elicited by Fasudil and Ligand 47 are similar, though with larger amplitudes around Y₁₂₅ and Y₁₃₆ by the latter; their correlation coefficient is 0.78 (Figure S3a). Because DIRseq does not consider any information about the drug molecule, this correlation coefficient between CSPs elicited by two different compounds might be viewed as an upper bound of what can be achieved for the correlation between CSP and DIRseq propensity.

DIRseq as a complement to CSP for assigning DIRs

After assessing the achievable accuracy of DIRseq, we now consider an application: its combination with CSP to make robust assignments of DIRs. We note that the above CSP-identified drug-interacting motifs are all anchored on one or more aromatic residues, and this feature likely contributes to the good performance of DIRseq. When a clear “aromatic signal” is not present, CSP or DIRseq alone may not be able to conclusively identify DIRs. However, a consensus identification by CSP and DIRseq may be reliable. MD simulations have played such a complementary role to CSP in several studies ^18,20,22,28.

De Mol et al. ¹¹ reported the CSPs of an IDR called AF-1* (residues 142-448) in the androgen receptor elicited by EPI-001 and its stereoisomers including EPI-002; small but reproducible CSPs were found in three subregions, R1 (residues 341-371), R2 (residues 391-414), and R3 (residues 426-446). Conversely, AF-1* caused changes in the ¹H NMR spectrum of EPI-001 but the peptides corresponding to R1-R3 did not, suggesting that, rather than separately interacting with the three subregions, EPI-001 induces a partial folding of the three subregions. MD simulations captured the partial folding of the R2-R3 fragment (residues 391-446) induced by EPI-002, and underscored the importance of aromatic residues in drug binding ²⁸. Basu et al. ²² reported the CSPs of the transactivation unit 5 (Tau-5*, residues 336-448) in AF-1* by EPI-001 and a more potent variant, 1aa (Figure 3a and Figure S3b). The contact probabilities of the R2-R3 fragment with 1aa in MD simulations reaffirmed the importance of aromatic residues.

Figure 3.

De Mol et al. ¹¹ and Basu et al. ²² were careful not to name any DIRs based on CSPs. In Figure 3a, we compare the 1aa-elicited CSPs and DIRseq propensities of Tau-5*. In agreement with both NMR studies, DIRseq identifies DIRs in the middle of each of R1-R3: R₃₆₀DYY₃₆₃, A₃₉₆WAA₃₉₉, and W₄₃₃H₄₃₄. In addition, the latter two motifs showed the highest drug-contact probabilities in separate MD simulations of the R2-R3 fragment ^22,28. For the R2 and R3 subregions, CSPs above the m + 1.5 SD threshold are observed at residues downstream of the DIRseq identifications, so we propose to expand the drug-interacting motifs to A₃₉₆WAAAAAQ₄₀₃ and W₄₃₃HTLF₄₃₇. The three putative drug-interacting motifs are indicated as cyan shading in Figure 3a.

Heller et al. ²³ used ¹⁹F transverse relaxation measurements to determine the binding affinity of the disordered domains 2 and 3 of the hepatitis C virus NS5A protein (NS5A-D2D3, residues 247-466) for 5-fluorindole. They also measured ¹H-¹⁵N CSPs at two ligand concentrations but described them as “nearly undetectable” (Figure 3b). We speculate that the small CSPs may be due to the small size of the ligand, making it difficult to interact with multiple residues simultaneously and thus achieve sufficient binding stability. In any event, CSPs above the m + 1.5 SD threshold were largely isolated (and thus appeared to be random), and there was not much overlap between residues having these above-the-threshold CSPs at the two ligand concentrations. The one exception is a motif around W₃₁₂, which had above-the-threshold CSPs at both ligand concentrations, and nearby residues A₃₀₈ and A₃₁₃ also had above-the-threshold CSPs at one of the ligand concentrations. DIRseq predicts the motif P₃₁₀AWARPD₃₁₆ as DIRs. We propose the expanded motif, A₃₀₈LPAWARPD₃₁₆, as residues interacting with 5-fluorindole (Figure 3b, cyan shading). DIRseq also predicts two more motifs, E₃₂₃SWRRPDY₃₃₀ and R₃₅₂RRR₃₅₅, as DIRs, which remain to be tested.

The fibrillation of acid-denatured β₂ microglobulin is inhibited by rifamycin SV ³⁵. An aromatic-rich motif, W₆₀SFYLLYYTEF₇₀, was implicated in the nucleation of fibrillation and also involved in ligand binding, as a triple mutation, F62A/Y63A/Y67A, significantly weakened binding. Low intensities of NMR peaks from residues 58-79 (possibly due to the formation of residual structures) prevented the measurement of CSPs. DIRseq predicts D₅₉WSFY₆₃ as DIRs. Combined with the mutational data of Woods et al. ³⁵, we propose the expanded motif, D₅₉WSFYLLYY₆₇, as the major drug-binding site (Figure 3c, cyan shading). DIRseq also predicts an additional motif, K₉₄WDR₉₇, as DIRs.

The aggregation of human islet amyloid polypeptide (hIAPP; 37 residues) is inhibited by the small molecule YX-I-1 ²¹. CSPs elicited by this molecule were small (Figure 4a). In addition, CSPs of short IDPs may not exhibit strong disparities because amino acids may be too well mixed along the sequence or drug binding may induce a conformational shift. We thus reduce the threshold for identifying DIRs to m + 1.0 SD when the number of residues is ≤ 50. With this threshold, three residues are identified by CSP as DIRs: R₁₁ and V₁₇H₁₈. In comparison, DIRseq identifies T₉QRLA₁₃ and F₁₅ as DIRs, which partially overlap with the CSP identifications. Combining the two types of data, we propose T₉QRLANFLVH₁₈ as the primary drug-binding site (Figure 4a, cyan shading). We note that this motif is also prone to α-helix formation ³⁶.

Figure 4.

Ono et al. ⁸ acquired ¹H-¹⁵N heteronuclear single quantum coherence spectra of the 42-residue amyloid-β (Aβ42) in the absence and presence of the oligomerization-blocking compound myricetin. CSPs were not calculated, but chemical shift movements were most pronounced at R₅, V₁₂H₁₃, K₁₆LVFFAED₂₃, and I₃₁I₃₂ (vertical lines in Figure 4b). In addition, NMR cross peaks suffered broadening upon ligand binding at four of these residues: R₅, V₁₂, K₁₆, and V₁₈ (dark vertical lines in Figure 4b), implicating elevated probabilities of ligand interactions. DIRseq predicts E₃FRH₆ and H₁₃H₁₄ as DIRs. Combined with the NMR data, we propose E₃FRH₆ and V₁₂HHQKLV₁₈ as the primary ligand-binding sites (Figure 4b, cyan shading). Aβ42 CSPs elicited by several other compounds were also widely distributed over the sequence, such that Iwaya et al. ₁₆ “failed to identify” DIRs, implicating a conformational shift.

Hammoudeh et al. ⁶ identified three distinct drug binding sites in a C-terminal IDR (residues 363-412) of the oncoprotein c-Myc. These authors measured the CSPs of three overlapping fragments (residues 363-381, 370-409, and 402-412) each in the presence of a single compound and of the full IDR in the presence of all three compounds. Using 20 parts per billion (ppb) as the threshold, Hammoudeh et al. named residues R₃₆₆RNELKRSFF₃₇₅, F₃₇₅ALRDQIPELE₃₈₅, and Y₄₀₂ILSVQAE₄₀₉ as the binding sites for 10074-G5, 10074-A4, and 10058-F4, respectively. We present their CSP data for the full IDR in Figure 4c. Using the m + 1.0 SD threshold, only six residues are identified as DIRs: R₃₆₇, L₃₇₀, F₃₇₄F₃₇₅, A₃₉₉, and S₄₀₅. DIRseq predicts two motifs, E₃₆₃RQRR₃₆₇ and K₃₇₁RSFFA₃₇₆, that overlap with the first two sites identified by CSP; a third motif, K₃₉₇KATAY₄₀₂, that corresponds to the third CSP-identified site has moderate DIR propensities. Combining CSP and DIRseq, we revise the three drug-binding sites to be E₃₆₃RQRRNELKRSFF₃₇₅, S₃₇₃FFALRDQI₃₈₁, and K₃₉₇KATAYILS₄₀₅ (Figure 4c, cyan, olive, and yellow shading, respectively). In addition to 10074-G5, 10074-A4, and 10058-F4, many compounds bind to these three sites 24,37–39.

Discussion

We have presented the first sequence-based method, DIRseq, to predict drug-interacting residues of IDPs. Assessment against NMR CSP demonstrates the accuracy of DIRseq. Drug-binding motifs are anchored on one or more aromatic residues, for forming π-π interactions with drug-like molecules that are rich in aromatic groups. The success of DIRseq comes without any specific information on the drug molecules, suggesting that IDPs may have a relatively simple sequence code for drug binding.

The notion that DIRs may be agnostic to the molecular details of drug compounds is supported by the fact that the same drug can bind to different IDPs. For example, in addition to p53 ¹⁸ (Figure 2c), the polyphenol EGCG also binds to many other IDPs, including α-synuclein ⁴⁰, hIAPP ^41,42, Aβ40 ¹³ and Aβ42 ⁴⁰, tau ⁴³, and merozoite surface protein 2 ⁴⁴. Likewise, 10074-G5 binds to c-Myc ⁶ (Figure 4c) but also to Aβ42 ²⁵. On the other hand, c-Myc represents a case where different compounds bind to distinct sites on a single IDP ⁶. A related example is presented by p27, where SJ403 typifies a group of compounds that share the same three binding sites (Figure 2a). Another group of compounds, typified by SJ710, binds only to the third site. Chemically, the presence of nitrogen atoms in the rings of SJ403 enhances its aromaticity and thus strengthens π-π interactions; in addition, the electronegative groups of SJ403 project into different directions, making it less restricted when forming electrostatic interactions (Figure S4). These features may explain why SJ403 can bind to all three sites whereas SJ710 can bind only to the third site, F₈₇YYR₉₀, where three consecutive aromatic residues followed by a basic residue ensure that SJ710 can form both π-π and electrostatic interactions.

We have illustrated the combination of DIRseq with NMR CSP to make robust identifications of drug-binding sites in IDPs. Indeed, CSP, MD simulations, and DIRseq are three orthogonal approaches that have great potential in complementing each other, not only for identifying drug-binding sites but also for elucidating the roles of amino acids, their sequence context, and different types of noncovalent interactions in forming such sites. DIRseq offers fast speed and a simple, direct link between sequence motif and propensity for drug binding. As more data from CSP and MD simulations become available, DIRseq can be tuned to make even better predictions. For example, it may be feasible to develop DIRseq versions specific to different pharmacophores.

Another application of DIRseq is to define IDP fragments for in-depth study by MD simulations, as shorter constructs both enable the use of a smaller simulation box and reduce the size of the conformational space. For example, based on CSP data and initial MD simulations of fill-length of -synuclein, Robustelli et al. ²⁰ chose a 20-residue C-terminal fragment for simulations of binding with additional compounds, leading to the identification of Ligand 47 as a stronger binder than the original Fasudil. Similarly, based on CSP data from longer constructs of the androgen receptor, both Zhu et al. ²⁸ and Basu et al. ²² chose the 56-residue R2-R3 fragment for MD simulations of drug binding. DIRseq can now play a similar role in selecting fragments for MD simulations when CSP data are unavailable. Longer constructs may also present challenges such as resonance assignments to NMR experiments, so well-chosen fragments guided by DIRseq can also benefit NMR studies.

Lastly, virtual screening has been conducted against conformational ensembles of IDPs ^17,45; drug-binding sites predicted by DIRseq can be used to guide such screening. As a simple illustration, we present poses of EGCG generated by screening against the two DIRseq-predicted binding sites in p53 in Figure 5. As IDPs sample a vast conformational space, knowledge of the binding site can drastically reduce the computational cost. The subset of conformations that generate high docking scores for a given drug at the known site can also provide insight into the mechanism of drug action.

Figure 5.

Acknowledgements

This work was supported by National Institutes of Health Grant GM118091.

Additional information

Computational Methods

The sequences of the IDP studied here and the drugs that bind to them are listed in Table S1. All DIRseq predictions were obtained using the web server at https://zhougroup-uic.github.io/DIRseq/. Conformations of IDPs were generated using the TraDES method ⁴⁶.

Docking of EGCG onto p53 was performed via the SwissDock web server at https://www.swissdock.ch/ ⁴⁷ utilizing the Autodock Vina docking engine ⁴⁸. The SMILES string for EGCG from PubChem (https://pubchem.ncbi.nlm.nih.gov/; CID 65064) and several conformations of p53 were used as input. A cubic region (13 to 20 Å in side length) around the center of each DIR was selected for docking.

Funding

National Institute of General Medical Sciences (GM118091)

Additional files

Supplementary Table and Figures