Introduction

Despite the groundbreaking impact of Alphafold2 (AF2) ^1,2 on the computational prediction of the ligand-free protein native apo structures, it appears that the determination of high quality crystal or cryo-EM ligands bound holo structures remains irreplaceable in the field of structure based drug design. When ligands bind, residues within the protein pockets may adjust their side-chain rotamer configurations to optimize contacts with ligands, a phenomenon known as the induced-fit effect. Furthermore, thermodynamic fluctuations induce protein dynamics and structural flexibility, leading to rearrangements of side-chains and even large-scale backbone movements that can reveal cryptic pockets.³ These metastable conformations and associated cryptic pockets can be stabilized upon binding to specific ligands, a phenomenon known as conformational selection. Moreover, due to the similarity of native structures among protein homologs, it has been widely believed that ligands targeting highly diverse metastable conformations should result in better selectivity.⁴ It is thus highly desirable to account for metastable protein conformations or states, instead of investigating native states only. Several AF2-based techniques, including reduced multiple sequence alignment (rMSA) AF2 (or MSA subsampling AF2),^5–7 AF2-cluster⁸ and AlphaFlow,⁹ have been devised to generate distinct decoy structures from native states. However, the suitability of these decoys for subsequent docking and virtual screening remains uncertain. Moreover, accurately assigning Boltzmann weights to decoys produced by those methods lacking a direct physical interpretation is challenging. Such a Boltzmann ranking is critical simply because of the explosion in number of decoys that can be hallucinated from AlphaFold2 or future generative AI methods, now including AlphaFold3.¹⁰

The AF2RAVE protocol integrates rMSA AF2 and the machine learning-based Reweighted Autoencoded Variational Bayes for Enhanced Sampling (RAVE) method^11–13 to systematically explore metastable states and accurately rank structures using Boltzmann weights.^14,15 Subsequently, traditional grid-based docking methods or recent generative diffusion models, like Diffdock ¹⁶ and DynamicBind¹⁷ can be employed to dock ligands with a few topranked structures in metastable states, enabling further virtual screening on large ligand libraries. In this work, we chose the docking method Glide XP^18–20and Induced-fit docking (IFD)^21,22 from the Schrödinger suite²³ as our primary docking method. By combining these two steps, we propose the AF2RAVE-Glide workflow (Fig 1) as an innovative approach for small molecule drug design, initiated from protein sequences. In this workflow, the conformational selection effect is addressed by probing various metastable states using AF2RAVE. Glide-IFD is then applied to account for the induced-fit effect, refine the holo-like pockets further, and predict the ligand-bound holo structures.

In this work, we demonstrate this AF2RAVE-Glide workflow retrospectively on three different protein kinases and their type I and type II inhibitors. Protein kinases are involved in the regulation of various cellular pathways by catalyzing the hydrolyzing of ATP and transferring the phosphate group to substrate peptides/proteins. Dysfunctions of various kinases are know to cause human pathologies and cancers. The human genome contains about 500 protein kinases which share highly conserved structures in their catalytic ATP-binding pocket, due to the selection pressure towards functional catalysis. This poses significant challenges in developing selective small-molecule ATP-competitive kinase inhibitors, as they must effectively target the intended kinase while avoiding off-target interactions and associated side effects. Research efforts aimed at achieving selectivity have led to the development of 4 distinct types of kinase inhibitors,²⁴ including two types of ATP-competitive kinase inhibitors: type I inhibitors bind to the ATP-binding site adopting the active catalytic state, while type II inhibitors target the binding pockets adjacent to the ATP-binding site adopting the inactive state. These two states primarily differ in the configurations of their activation loop (A-loop), which is a flexible loop of approximately 20 residues. In the active state, the A-loop is “extended” to create a cleft for substrate peptides to bind, while in the inactive state, it is collapsed or “folded” onto the protein surface, blocking substrate binding. Additionally, in the active state, the three residues ‘Asp-Phe-Gly’ (DFG motif) at the N-terminus of the A-loop bind to an ATP-binding Mg²⁺ ion, with the Asp side-chain pointing inward to the ATP-binding pocket, while in the inactive state, the Asp side-chain is flipped outward, and the DFG-motif adopts the DFG-out conformation.²⁵ Henceforth, we will refer to the type I inhibitor binding state as the active DFG-in state and designate the type II inhibitor binding state as the classical DFG-out state, as previously proposed in the literature^26–28 (demonstrated in Fig 2 A and B using Abl1 kinase as an example).

In this work, we investigated three kinases: (i) Abl1, which is targeted by the first clinically approved small-molecule kinase inhibitor imatinib for cancer therapy; (ii) DDR1, a structurally more flexible kinase which is identified as a promiscuous kinase targeted by chemically diverse inhibitors; ²⁹ and (iii) Src kinase, another crucially important member of the tyrosine kinase family. In the following section, we applied AF2RAVE to enrich hololike structures adopting the classical DFG-out state from the AF2-generated ensembles of DDR1 and Abl1 kinase, and further validated those holo-like structures by docking them with known type II kinase inhibitors. We initiated our investigation by conducting docking experiments involving both type I and type II inhibitors with the AF2 structures of these two kinases. We observed the incapacity of AF2 structures to effectively dock with type II inhibitors targeting the metastable classical DFG-out state. Subsequently, we employed rMSA AF2, an associative memory like process that generates diverse structure ensembles for both kinases^7,30, exploring their relatively improved but still limited potential to contain holo-like structures suitable for type II inhibitor binding. Ultimately, our study showcases AF2RAVE’s effectiveness in enhancing the generation and selection of hololike structures in metastable states by integrating AF2-based ensemble generation with physics-based methods.

Figure 1

Results and discussion

AF2 structures fail to dock with type II kinase inhibitors

The utility of AF2-generated structures for structure-based drug discovery and virtual screening campaigns has been a subject of controversy and skeptical enthusiasm.^31–34 For proteins lacking crystal and cryo-EM structures but with homologous structures available in PDB, such as CDK20 kinase, AF2 demonstrates effectiveness as a homology modeling method for generating initial structures suitable for subsequent virtual screening.³⁵ Other AF2-based homology modeling method can bias AF2-generated structures towards user-selected template structures with specific druggable conformations.³⁶ However, a significant drop in the hits enrichment factor during virtual screening has been reported when employing AF2 structures as rigid receptors for docking, compared to using holo PDB structures. This occurs even in cases where the binding pockets of AF2 structures differ only slightly at 2 to 3 residues, from those of holo PDBs. ^33,34 Given that ligands can induce slight relocation and side-chain rotation of pocket residues upon binding, it is important to note that AF2 does not account for this induced-fit effect, as it does not encode co-factors like ligands. Therefore, it appears necessary to perform ligand induced-fit modeling or relaxation on AF2 structures before engaging in any further structure-based drug design. Molecular dynamics (MD) simulations biased towards adjacent holo-template structure have proved effective in refining apo structures and improving their early enrichment performance. ³⁷ It has also been demonstrated that AF2 structures can achieve comparable accuracy to crystal holo structures in Free Energy Perturbation (FEP) calculations, by superposing AF2 structures with crystal structures, grafting the co-crystallized ligands onto the AF2 structure, and optimizing the AF2 structure/ligand complex to account for subtle induced-fit effects.³⁸ AF2 structures, decorated with ligands from template ligand grafting method or known-hits docking method, and further refined by Schrödinger IFD-MD protocol, exhibit promising performance in early enrichment³⁹ and the prediction of novel ligand/protein complex structures. ⁴⁰ Therefore, if large backbone motion is not required, it appears feasible to refine the apo pockets in AF2 structures into holo-like pockets for structure based drug design. However, when AF2 structures exhibit significant steric clashes with holo ligands, especially in cases ligands targeting metastable states, direct application of out-of-the-box AF2 structures in docking methods for virtual screening and early enrichment may pose more challenges. The AF2 predicted kinase structures predominantly exhibit the DFG-in state, with over 95% of human kinases predicted in this conformation.^41,42 Significant A-loop motion and backbone flipping of the DFG motif are necessary to transition from the AF2-predicted DFG-in state to holo-like states, for type II ligands targeting the classical DFG-out state. As a result, AF2 structures of Abl1/DDR1 kinases exhibit superior performance when docking with type I inhibitors (achieving a minimum ligand RMSD of 2.14 Å) compared to docking with type II inhibitors using Glide XP (Fig 2 C&D). Even with Glide’s Induced-fit Docking (Fig S16) and the highly side-chain clashes forgiving docking method DiffDock (Fig S20), AF2 structures struggle to dock with these metastable state-targeting ligands (type II inhibitors), with ligand RMSDs above 8 Å across all docking poses from all three docking methods.

Holo-like metastable structures may be present among decoys generated from rMSA AF2

AF2-based methods can achieve structural diversity by introducing dropouts in MSA inputs stochastically (rMSA AF2) or in a clustering manner (AF2-cluster). Additionally, models modified from the AF2 framework, such as the flow-match generative model AI- phaFlow,⁹ have been developed to explore the diversity of conformational space. Similar protocols to rMSA AF2 have demonstrated the potential to address the induced-fit effect by generating diverse structures at the binding pocket, ranging from closed apo pockets to opened holo-like pockets. ⁴³ Other investigations have also indicated that larger backbone motions, such as DFG-motif backbone flipping and A-loop movement in at least some protein kinases, can be captured by the rMSA AF2 ensemble, although the distributions of conformations deviate significantly from the correct Boltzmann distributions. ^6,15

Figure 2

In this section, we utilized rMSA AF2 to generate 1280 diverse structures for Abl1 kinase or DDR1 kinase: 640 for MSAs of depth 16 and 32. See Supplementary Material for information regarding Src kinase. Following a filtering step based on RMSD from corresponding AF2 structures, 1198 and 1147 structures remain for Abl1 and DDR1, respectively (Fig S1). As shown in Fig 3 A&B, for Abl1 kinase, only 4 out of 1198 structures have a folded A-loop, when using a distance cutoff of 15 Å between CB atoms of N98 and R162 in Abl1. However, for DDR1 kinase, 124 out of 1147 rMSA AF2 structures exhibit a folded A-loop, when using a salt bridge distance cutoff of 10 Å between the aligned residue pairs in DDR1 (E110 and R191). We then clustered the A-loop folded structures in the Dunbrack space. For Abl1, only two clusters were identified: one adopts the DFG-in state, and the other adopts the DFG-inter state, referring to an intermediate conformation during the backbone flipping from DFG-in to DFG-out, according to the Dunbrack definition. For DDR1, structures were divided into 5 clusters, with one cluster of size 15 being the closest to the classical DFG-out state.

We subsequently docked type II inhibitors, ponatinib and imatinib, to the most “classical DFG-out”-like cluster from rMSA AF2 ensembles (indicated as red circles in Fig 3 A&B), using IFD. Despite the A-loop being folded, the DFG-motif in structures from the Abl1 cluster is not strictly DFG-out. As expected, type II ligands fail to dock with those structures, with all IFD poses exhibiting ligand RMSDs above 9 Å (Fig S16 A, Fig 5 B). While AF2-based methods can indeed generate decoy structures that deviate from the native state, it remains uncertain whether these decoys correspond to metastable basins. Additionally, it’s unclear whether these decoys includes structures that can represent the specific metastable states required for the intended types of drug design.

In contrast to Abl1, rMSA AF2 ensemble for DDR1 contains holo-like structure for type II inhibitors. One structure (which we label the “holo-model”, shown as a red circle filled with green in Fig 3 B), out of the 15 in the DDR1 classical DFG-out cluster docked pona-tinib with a remarkably low RMSD of 0.89 Å using IFD (Fig 3 C). However, IFD poses from all the other structures exhibit large ligand RMSDs above 6 Å (Fig S16 B). Hence, an enrichment process to select holo-like structures from decoys becomes essential to ensure a practical number of pocket structures for ensemble docking and virtual screening.

We have also docked the “holo-model” structure with another type II ligand, imatinib. In this case, the steric clashes are more significant, rendering the refinement for the induced-fit effect more challenging compared to the ponatinib scenario. We thus introduced an extra trimming step (see SI text for details) for the DFG-Phe residue which exhibit significant clashes with the holo-imatinib in the “holo-model” structure (Fig S17 B). The IFD-trim protocol successfully achieves a minimum ligand RMSD of 1.04 Å for docking poses of the “holo-model” structure with imatinib (Fig 3 C, Fig S17 C&D).

AF2RAVE on DDR1 enriches holo-like classical DFG-out structures in rMSA AF2 decoys

To assess the utility of physics-based methods in selecting holo-like structures from AF2- generated ensembles, we utilized a physics-based protocol, AF2RAVE¹⁴ to explore the energy landscape of DDR1 kinase. We employed the identical set of collective variables (CVs) as in our previous study¹⁵ to perform regular space clustering for the DDR1 rMSA AF2 ensemble. Two structures were selected from the clustering centers for each combination of DFG type (in, inter, or out) and A-loop position (folded or extended), resulting in a total of 12 initial structures (Fig S2). Subsequently, 50 ns unbiased MD simulations were conducted starting from each initial structure. The CVs extracted from all MD trajectories were input into the SPIB model with a linear encoder, to learn the reaction coordinates indicative of the slow motions linked to transitions between metastable states. The learnt SPIB reaction coordinates possess clear physical interpretations, with the first coordinate signifying the A-loop position and the second indicating the DFG-type (Fig 4 A&B). This alignment with physical features is expected as we deliberately selected diverse and representative initial structures for unbiased MD simulations. Additionally, it’s noteworthy that these reaction coordinates are transferable across various kinases and can effectively discern different metastable states.

Figure 3

Interestingly, the 15 classical DFG-out structures within the DDR1 rMSA AF2 ensemble are situated in a region of the latent space that were not thoroughly explored by the 12 unbiased MD trajectories. To address this gap, we employed enhanced sampling to sample along the SPIB-approximated reaction coordinates and compute the free energy profile inside the classical DFG-out basin. Considering both the flipping of the DFG motif and the overall motion of the large flexible A-loop, it might be impractical to sample direct back-and-forth transitions between various states using metadynamics. Therefore, we opted for umbrella sampling for its simplicity. The reliability of umbrella sampling hinges on two issues, first whether the latent space adequately represents the conformational space and second, whether there is sufficient overlap between different windows for efficient reweighting. Addressing the first challenge remains an ongoing endeavor in the dimensionality reduction research field, and we anticipate that our SPIB latent space is sufficient enough for our current purpose. The second challenge can be managed through careful setup of umbrella sampling windows and bias strength. Given our current setup, sampling the extensive motion of A-loop relocation remains challenging, showing insufficient overlap between the A-loop folded and extended regions (Fig S3). Consequently, the quantitative reliability of the absolute ΔG values between different states is limited, allowing us only to qualitatively assess the relative thermodynamic stability of the DFG-in versus the DFG-out basins. Nevertheless, the qualitative relative stability observed from umbrella sampling aligns with previous studies^15,29 (Fig S4). Furthermore, the local potential of mean force (PMF) surrounding each basin should provide quantitatively reliable insights, given their thorough sampling through umbrella sampling. Reassuringly, when we ranked the classical DFG-out structures within the DDR1 rMSA AF2 ensemble, using the local PMF values in the latent space, the “holo-model” structure emerged among the top 2 structures with free energy relative to the minimum smaller than 1 kJ/mol, as illustrated in Fig 4 D.

Additionally, we used DiffDock to conduct docking experiments with type II ligands on the 15 classical DFG-out structures within the DDR1 rMSA AF2 ensemble. Despite poses from most structures surprisingly demonstrating very low ligand RMSD, typically below 2 Å (Fig S20), they exhibited significant steric clashes and low DiffDock confidence score due to Diffdock’s disregard for side-chain configurations and clashes. Notably, upon ranking structures based on their AF2RAVE PMF and examining the corresponding poses with the lowest ligand RMSD, we observed that only the poses derived from the top 2 structures selected by AF2RAVE had Diffdock confidence scores higher than −1.5, surpassing the default threshold of the DiffDock model (Table S1). This implies the general advantage of structures selected by physics-based methods across various docking methods.

In summary, Boltzmann ranking from AF2RAVE effectively distinguishes holo-like structures from other decoys. Beginning with the DDR1 rMSA AF2 ensemble, AF2RAVE substantially enhances the likelihood of identifying the holo-like structure from 1 out of 15 to at least 1 out of 2 when using a PMF cutoff of 1 kJ/mol. We noted that although the PMF values and Boltzmann ranks may fluctuate with the setup of umbrella sampling, the enrichment effect of holo-like structures remains (Fig S7). Besides, we applied an alternative protocol to compute the PMF profile within the DDR1 classical DFG-out basin, by running unbiased MD simulations starting from the 15 classical DFG-out decoys in the DDR1 rMSA AF2 ensemble. However, this unbiased protocol poses a risk of failing to sample rare events if barriers exit within the region of interest. In this case, the minibarriers inside the classical DFG-out basin are low enough to enable efficient sampling using 50 ns unbiased MD simulations. Consequently, the Boltzmann ranks derived from unbiased MD simulations also successfully enriched the single holo-like structure in the DDR1 rMSA AF2 ensemble to the top 5 of the structure set (Fig S11).

Transferable learning of holo-like structure for Abl1 and Src kinases with AF2RAVE from DDR1 templates

As mentioned earlier, our current Abl1 rMSA AF2 ensemble lacks any decoy structure in the classical DFG-out state. This is further illustrated in Fig 5 A, where the Abl1 rMSA AF2 ensemble is projected onto the same latent space learnt during the DDR1 AF2RAVE protocol. Hence, it’s necessary to prepare the Abl1 decoy set adopting the classical DFG-out state before we can rank them using Boltzmann weights derived from physics-based methods and conduct further docking for selected holo-like structures.

There are several approaches to generate the Abl1 decoy sets. One can conduct enhanced sampling on the latent space starting from Abl1 rMSA AF2 structures to reach the classical DFG-out basin. Subsequently, MD structures from this basin can serve as templates for asking AF2 to generate crystal-like structures in classical DFG-out state. However, for simplicity, we opted to use the 15 classical DFG-out structures from the DDR1 rMSA AF2 ensemble directly as templates and employed an AF2-based homology modeling protocol, referred to as AF2-template (tAF2 for short, detailed protocol can be found in the SI text), to generate a decoy set comprising 30 Abl1 structures, as illustrated by the green stars in Fig 5 A.

Compared to the AF2 and rMSA AF2 structures, the performance of IFD on type II inhibitors shows a significant improvement when using the 30 tAF2 decoy structures of Abl1. The lowest ligand RMSD achieved is 2.74 Å for imatinib and 0.78 Å for ponatinib (Fig 5 B). However, only 4 structures (labeled as “holo-model” structures hereafter) out of the 30 decoys are capable of docking with type II inhibitors with ligand RMSD < 3 Å, and all other structures produce IFD poses with ligand RMSD> 6 Å (Fig S18).

Figure 4

We then investigated whether physics-based methods could enrich the 4 “holo-model” structures from the 30 tAF2 decoys to a practical number of holo-like candidates. To explore the local classical DFG-out basin, we conducted 50 ns unbiased MD simulations starting from each tAF2 decoy structure and combined the trajectories to obtain the final Boltzmann distribution. Remarkably, the PMF for the Abl1 classical DFG-out basin calculated from the unbiased protocol (Fig 5 C) showed the enrichment of all 4 “holo-model” structures among the top 8 tAF2 structures with PMF value < 1 kJ/mol. The PMF values of the top 8 tAF2 structures and the lowest ligand RMSD from the corresponding structures are presented in Fig 5 D. We also employed umbrella sampling for the Abl1 kinase by transferring setups and initial structures using the AF2-template from DDR1 umbrella sampling. However, we noted a significant occurrence of αC-helix breakage in the Abl1 umbrella sampling trajectories compared to DDR1 (Fig S5). After excluding windows with broken αC-helix, the Abl1 PMF derived from umbrella sampling ultimately enriched all 4 “holo-model” structures to the top 6 among the 30 tAF2 structures (Fig S6). In Supplementary Information, we also provide results for Src kinase which are in the same quality as for Abl1 kinase.

Discussion

Through our retrospective analysis, we have thus demonstrated that the default AlphaFold2 models are ineffective for docking ligands targeting metastable protein kinase conformations. While AF2-based methods can be coaxed into generating diverse structures, they still struggle to produce reliable accuracy decoys for metastable conformations since the AF2 ensembles do not follow Boltzmann distribution. This failure is evident in the inability to generate an Abl1 AF2-ensemble containing holo-like structures for type II inhibitors. To further investigate whether this limitation is common among AF2-based methods, including the rMSA AF2 method we employed earlier, we tested another AF2-based approach, AF2-cluster.⁸ The AF2-cluster ensemble of Abl1 comprises more A-loop folded structures (21 out of 197) compared to our Abl1 rMSA AF2 ensemble (4 out of 1198). However, similar to rMSA AF2, there are still no decoys in the classical DFG-out state, and all A-loop folded structures are located far from the PMF basin of the classical DFG-out state in the latent space (Fig S14 C). This indicates that AF2-cluster, like rMSA AF2, also fails to generate Abl1 metastable states effectively. Interestingly, we observed comparable structural diversity in sampling the A-loop folded configurations within the AF2-cluster ensembles for DDR1 and Abl1 (Fig S14 C&D). While for the rMSA AF2 method, the promiscuous kinase DDR1 ensemble exhibit superior structural diversity compared to the Abl1 kinase. This enhanced diversity leads to the identification of one dockable structure for type II inhibitors among decoys in DDR1 rMSA AF2 ensemble. Through the application of a homology modeling method, AF2-template, we demonstrated that the classical DFG-out decoys in the DDR1 rMSA AF2 ensemble can be transferred to Abl1 kinase. Furthermore, we tested an additional kinase, Src, for which non-native state decoys are reported to be even more challenging to produce using AF2 subsampling methods than Abl1 kinase.⁶ We have also verified that rMSA AF2 and AF2-cluster struggle in producing distinct decoys from native structure of Src kinase (Fig S1 & Fig S12). Besides, neither the rMSA AF2 nor the AF2- cluster ensembles for Src kinase adequately sample the classical DFG-out basin in the latent space (Fig S13). Remarkably, AF2-template, as a homolgy modelling method, can easily produce the classical DFG-out structure of Src kinase, using the top 2 AF2RAVE-picked DDR1 classical DFG-out structures as templates (Fig S9). The IFD poses from the tAF2 structure of Src kinase shown a minimal ligand RMSD of 2.82 Å with imatinib (Fig S19).

Figure 5

With the rapid expansion of chemical space, virtual screening on libraries containing billions of diverse molecules becomes enticing for novel drug discovery.⁴⁴ Therefore, the enrichment of candidate holo-like structures emerges as a necessary step, offering significant benefits in terms of computational efficiency and feasibility. As summarized in Table 1, unlike the non-dockable AF2 structures for type II inhibitors, the diverse rMSA AF2 ensemble (referred to as rAF2 in Table 1 for brevity) shows potential in generating holo-like structures within a large set of decoys. However, it’s only upon AF2RAVE ranking and selection that the ratio of holo-like structures in selected structure set increases to a plausible value of 50%, facilitating further virtual screening on computational models of protein pockets.

Table 1

Conclusion

AlphaFold2 has arguably revolutionized protein structure prediction, but it remains to be constructively demonstrated if it can be reliably used for drug discovery purposes, especially involving non-native protein conformations. In this work we have demonstrated through retrospective studies on kinase inhibitors that a combination of AlphaFold2, statistical mechanics based enhanced sampling and induced fit docking can be deployed for such calculations. Specifically, we have utilized the AF2RAVE protocol by inputting the sequences of the DDR1 kinases, along with two additional pieces of prior information: a pairwise distance cutoff for evaluating A-loop positions and the Dunbrack definition of DFG-type. We then employed Glide Induced-fit Docking (IFD) to assess our AF2RAVE-generated computational models for type II kinase inhibitor binding pockets. This AF2RAVE-Glide workflow yielded holo-like structure candidates with a 50% successful docking rate for known type II inhibitors. Notably, the holo-like structures in metastable state and the latent space constructed from AF2RAVE of DDR1 are transferable to other kinases. This includes the challenging cases⁶ of Abl1 and Src kinases, wherein we showed that SPIB and sampling performed for DDR1 allowed generating classical DFG-out structures for both Abl1 and Src kinases. This severely reduces the computational cost for retraining SPIB to learn low-dimensional latent space for different kinases.

This demonstration of AF2RAVE-Glide on kinase inhibitors shows its promising application for discovering drugs targeting general proteins in addition to kinases, such as G-Protein Coupled Receptors (GPCRs), which are the targets for over one-third of Food and Drug Administration (FDA)-approved drugs.⁴⁵ For the design of novel drugs targeting general proteins, developing a protocol that does not require prior information about the system is left for future exploration. Besides, in this study we only investigate the classical DFG-out metastable state for kinases. For a comprehensive protocol, all top-ranked metastable states identified by AF2RAVE should be explored in subsequent docking experiments. Integration of algorithms capable of predicting ligand binding sites on protein surfaces, such as the Graph Attention Site Prediction (GrASP) model,⁴⁶ is then essential before utilizing AF2RAVE-selected structures in docking, thus expanding the workflow to AF2RAVE-GrASP-Glide. Additionally, the inclusion of free energy perturbation (FEP) calculations for front-runner ligands to evaluate the actual ligand binding affinity can further enhance this workflow.

The integration of AF2-based and physics-based methods presents a promising approach toward the development of a mature workflow for computer-aided drug design. AF2-based methods are capable of producing ensembles with structural diversity, which aids physics-based methods in better sampling and exploring the energy landscape of proteins. Additionally, AF2-generated structures can serve as crystal-like decoys, free from distortion that may occur in biased simulations. Physics-based methods play a crucial role in accurately assigning Boltzmann weights and ranking decoy structures to guide the enrichment of holo-like structures, essential for virtual screening on large libraries. This collaborative approach leverages the strengths of both methodologies, leading to enhanced efficiency and efficacy in drug discovery efforts.

We conclude this manuscript by making a final comment on the horizons opened by Generative AI methods, including those involving diffusion models^10,47 and any such future frameworks. These approaches make it possible to easily hypothesize regions of the conformational space underlying arbitrarily complex molecules of life, that then could serve as a starting point to launch more careful investigations. However, these predictions without careful in situ or a posteriori testing through advanced simulations or experiments, are only predictions and could just be hallucinations. Having the capability to quickly generate numerous - thousands or more - such structural hypotheses is what made AlphaFold2 so crucial to this current work. We believe that a deep integration of the hypothesis creation possibilities of Generative AI with careful Molecular Dynamics, ^48,49 experiments or other forms of rapid testing is one way these methods will truly facilitate new and reliable scientific discoveries.

Methods

Summary of systems tested and tools used

In this work we tested AF2RAVE-Glide protocol on active and inactive conformations of DDR1, Abl1 and Src kinases against type I inhibitors VX-680 and Dasatinib, and type II inhibitors Imatinib and Ponatinib. The top-ranked structures obtained from AF2RAVE-Glide were compared against publicly available crystal structures deposited in the Protein Data Bank (PDB), with PDB codes provided in the main text.

To generate diverse structural ensembles for DDR1, Abl1, and Src kinases, we primarily used the reduced MSA AF2 (rMSA AF2) method. For comparison, we also produced AF2-cluster ensembles for these three kinases.

For MD simulations, all the systems in this paper were parameterized with the Amber99SB*- ILDN force field^50,51 with the TIP3P water model⁵² and neutralized with Na⁺ ions or Cl⁻ ions. The simulations are performed at 300 K with the LangevinMiddleIntegrator⁵³ in OpenMM⁵⁴ with the step size of 2 fs. Particle Mesh Ewald⁵⁵ is used for calculating electrostatics and the lengths of bonds to hydrogen atoms are constrained using LINCS⁵⁶ throughout all simulations. Before performing MD simulations for analysis, energy minimization is conducted for all initial structures, followed by equilibrium runs under NVT and NPT for 500 ps and 1 ns respectively.

To account for the induced-fit effect, Induced-Fit Docking (IFD) from the Schrödinger suite is the primary docking method used in this work. For comparison, we also tested the Glide XP docking from the Schrödinger suite and DiffDock to dock the two known type II inhibitors (Imatinib and Ponatinib) against DDR1 AF2 structure and the 15 classical DFG-out structures in the DDR1 rMSA ensemble.

Generation of rMSA AF2 ensemble, AF2-cluster ensemble and AF2- template structures

In this work, Colabfold² was employed to generate all the AF2-based ensembles and structures. The multiple sequence alignments (MSAs) were produced using mmseq2. For rMSA AF2 ensembles, the MSA depth was reduced to either 16 or 32 for each kinase. For each depth, 128 random seeds were utilized, and each seed produced 5 structural models via Colabfold, resulting in a total of 1280 rMSA AF2 structures per kinase. The structure with the highest pLDDT score among all 1280 was identified as the AF2 structure (native structure). Subsequently, any unphysical structures with an RMSD greater than 7 Å from the AF2 structure were discarded, culminating in the final rMSA AF2 ensembles.

AF2-cluster for DDR1, Abl1 and SrcK were run with default setups as provided in the ColabFold notebook in the original paper of AF2-cluster.⁸

The AF2-based homology modeling protocol, AF2-template (tAF2) method is implemented using Colabfold. For a given query sequence, tAF2 structures are generated by Colabfold upon uploading desired template structure and deactivating the Evoformer module. For each template, 5 AF2-template models are generated, and the last 3 structures exhibiting lower pLDDT will be discarded. The AF2-template (tAF2) method is employed here to transfer structure sets between homologous systems. To generate tAF2 structures for Abl1 in the classical DFG-out state, each of the 15 classical DFG-out structures from the DDR1 rMSA AF2 ensemble is used as a template, resulting in 30 tAF2 Abl1 structures in total. For Src kinase, a single representative tAF2 structure is generated using the “holo-model” DDR1 structure as the template.

AF2RAVE protocol

Regular space clustering on the rMSA AF2 ensemble

We used the same 14 collective variables (CVs) for regular space clustering as in the previous AF2RAVE work on kinases. ¹⁵ These CVs are pariwise distances, selected to describe the kinase conformations around the ATP-binding pocket and the A-loop.

Unbiased MD and SPIB

50 ns unbiased MD simulation was run for each AF2RAVE initial structure of DDR1 from Fig S2. The standard deviations of the 14 CVs are calculated after concatenating all the 12 unbiased MD trajectories. 8 CVs with standard deviations larger than 0.25 of the maximum standard deviation remain as the input features of the SPIB model. We conducted a parameter screening of the SPIB time lag, ranging from 1 ns to 40 ns with intervals of 1 ns. Eventually, we selected a time lag of 16 ns based on the performance of SPIB coordinates in representing physical features, including DFG-type and A-loop position.

PMF calculations from umbrella sampling

2D umbrella sampling is conducted along the two learnt SPIB coordinates, employing 11x11 windows. The bias potential equilibrium points are uniformly distributed in the SPIB latent space, with σ₁ ranging from −0.8 to −0.1 and σ₂ ranging from 0 to 0.8. The strength of the bias potential is set to 1000 kJ/mol/nm². Each window originates from the structure closest to the window’s equilibrium point in Euclidean distance within the latent space among the 12 AF2RAVE initial structures and lasts for 100 ns (with the first 10 ns discarded in PMF calculation). For Abl1 kinase, we noticed a substantial proportion of the αC helix breaking in the umbrella sampling trajectories (Fig S5). This finding aligns with earlier enhanced sampling investigations on DDR1 using metadynamics, ¹⁵ where the authors imposed restraints to prevent αC helix breakage. In our study, we opted to exclude all umbrella sampling windows where the ratio of broken αC helix exceeded 20% for the Abl1 PMF calculation (Fig S6). The WHAM algorithm is applied to bin and reweight the biased trajectories and compute the final PMF.

PMF calculations from unbiased MD

50 ns unbiased MD simulation was run starting from each structure in the 15 classical DFG-out structures in DDR1 rMSA AF2 ensemble. Upon discarding first 10 ns, all the unbiased trajectories are simply concatenated to calculate the Boltzmann distribution and PMF for each bin around the classical DFG-out basin in the latent space. For Abl1, unbiased simulations start from 30 AF2-template structures to calculate the PMF around the classical DFG-out basin.

Boltzmann ranks assignment for structures in AF2-based ensembles

After calculating the PMF value for each bin in the latent space, we projected AF2- generated structures into latent space. We then directly assign the PMF values of the corresponding bins to these AF2-generated structures.

We must acknowledge the limitations of the way we assigned PMF values to AF-generated candidate holo structures. First of all, the free energy profiles are derived from MD simulations, and the PMF values directly correspond to the MD structures. Here, we assumed that the latent space adequately represents the conformational changes of protein pocket within specific metastable states. Additionally, while the enrichment of holo structures in the top Boltzmann-ranked structures persists, the absolute PMF values and Boltzmann ranks of proper holo structures may fluctuate with the umbrella sampling setups, as depicted in Fig S7. Theoretically, the number of umbrella sampling windows and the simulation length should be sufficiently large for PMF convergence. However, there is always a trade-off between PMF accuracy and computational costs, so we opted to stick with the current setups.

Docking details

All the input structure for our docking experiments were first relaxed in solution with a MD energy minimization step. This work investigates two type I inhibitors (VX-680 and dasatinib) and two type II inhibitors (imatinib and ponatinib) by docking. For Glide XP docking or IFD, we used Ligprep in Maestro to prepare the ligand inputs from SMILES files. For DiffDock, the ligand inputs were directly provided as SMILES files.

Glide XP Docking

Glide XP docking experiments in this work were run with default setups in the Maestro software. Glide XP docking was performed for all four ligands on the AF2 structures of DDR1 or Abl1, as well as on the 15 classical DFG-out conformations of DDR1 in the rMSA AF2 ensemble.

Induced-Fit Docking (IFD)

For Induced-fit Docking, we used the Glide XP for initial docking, followed by Prime relaxation and final Glide XP docking. Parameters remain default in Maestro IFD. IFD was performed only for the type II ligands on the AF2 structures of DDR1 or Abl1, on the classical DFG-out conformations of DDR1 or Abl1 in the rMSA AF2 ensembles, as well as tAF2 structures of Abl1 and SrcK.

Given ponatinib’s backbone features, notably its lengthy and slender carbon-carbon triple bond, it exhibits reduced sensitivity to steric clashes, resulting in successful docking with the DDR1 “holo-model” structure at a ligand RMSD of 0.89 Å. Conversely, the “holomodel” structure struggles to accurately dock with imatinib using IFD (Fig S15B). We then employed an extended-sampling IFD approach by initially trimming the DFG-Phe residue from the “holo-model” structure. The trimmed residue is temporarily mutated to alanine during initial docking and later restored in subsequent Prime relaxation and final docking steps. Here, in the “holo-model” structure, we manually chose the DFG-Phe residue which is significantly hindered by holo-imatinib (Fig S16B). For generic systems lacking groundtruth information, broader screening of single residue trimming for the protein pocket may be necessary for this extended-sampling IFD method. In essence, it’s a tradeoff between the quality of holo-like structure to dock with and the accuracy/complexity of the docking method.

Diffdock performace on DDR1 classical DFG-out in rMSA AF2 ensemble

DiffDock docking experiments in this work were run in the webserver with default setups in the version before 3/8/2024 (https://huggingface.co/spaces/simonduerr/diffdock). ¹⁶ DiffDock was performed only for type II ligands on the AF2 structures of DDR1, as well as on the 15 classical DFG-out conformations of DDR1 in the rMSA AF2 ensemble.

Acknowledgements

Research in this publication was supported by the National Institute of General Medical Sciences of the National Institutes of Health under Award Number R35GM142719. The content is solely the responsibility of the authors and does not represent the official views of the National Institutes of Health. We thank UMD HPC’s Zaratan and NSF ACCESS (project CHE180027P) for computational resources. A.A. was supported by NCI-UMD Partnership for Integrative Cancer Research. P.T. is an investigator at the University of Maryland-Institute for Health Computing, which is supported by funding from Montgomery County, Maryland and The University of Maryland Strategic Partnership: MPowering the State, a formal collaboration between the University of Maryland, College Park and the University of Maryland, Baltimore. We thank UMD HPC’s Zaratan and NSF ACCESS (project CHE180027P) for computational resources. We thank Bodhi Vani, Dedi Wang, Zachary Smith and Anjali Verma for helpful discussions.

Notes

The authors declare the following competing financial interest(s): P.T. is a consultant to Schrodinger, Inc. and is on their Scientific Advisory Board.

Supplementary Information

Figure S1

Figure S2

Figure S3

Figure S4

Figure S5

Figure S6

Figure S7

Figure S8

Figure S9

Figure S10

Figure S11

Figure S12

Figure S13

Figure S14

Figure S15

Figure S16

Figure S17

Figure S18

Figure S19

Figure S20

Table S1