Adaptive simulations, towards interactive protein-ligand modeling

Abstract

Modeling the dynamic nature of protein-ligand binding with atomistic simulations is one of the main challenges in computational biophysics, with important implications in the drug design process. Although in the past few years hardware and software advances have significantly revamped the use of molecular simulations, we still lack a fast and accurate ab initio description of the binding mechanism in complex systems, available only for up-to-date techniques and requiring several hours or days of heavy computation. Such delay is one of the main limiting factors for a larger penetration of protein dynamics modeling in the pharmaceutical industry. Here we present a game-changing technology, opening up the way for fast reliable simulations of protein dynamics by combining an adaptive reinforcement learning procedure with Monte Carlo sampling in the frame of modern multi-core computational resources. We show remarkable performance in mapping the protein-ligand energy landscape, being able to reproduce the full binding mechanism in less than half an hour, or the active site induced fit in less than 5 minutes. We exemplify our method by studying diverse complex targets, including nuclear hormone receptors and GPCRs, demonstrating the potential of using the new adaptive technique in screening and lead optimization studies.

Figure 1

Protein-ligand complexes studied. (a) Trypsin with benzamidine as a ligand (TRP, PDB ID: 3PTB). (b) Progesterone nuclear hormone receptor with progesterone as a ligand (PR, PDB ID: 1A28). (c) Corticotropin-releasing factor GPCR with CP-376395 as a ligand (B-GPCR, PDB ID: 4KY5). (d) M3 muscarinic acetylcholine GPCR with tiotropium as a ligand (A-GPCR, PDB ID: 4DAJ). The initial structures for the protein-ligand exploration, with the ligand ~20 Å away from the binding site, are shown. The square inset in each panel depicts a 2D scheme of each ligand.

Figure 2

Energy landscape exploration of B-GPCR with 512 different explorers. (a,b) The RMSD variation along MC steps and the binding energy against the RMSD for the adaptive results. Each color code corresponds to a different epoch number, for a total of 12 adaptive iterations. (c,d) Analogous plots for the standard executions. Each color corresponds to a different trajectory (performed in a different computing core). Notice the change in scale in the X-axis between (a) and (c).

Figure 3

Binding times for all systems and MC techniques. (a) Number of steps for observing a binding event against the number of trajectories (processors) for the TRP system, using the standard PELE (in red) and the adaptive-PELE with the inversely proportional (in blue) and the ε-greedy guided strategies with binding energy (in green) and RMSD (in orange). Actual data (MC steps) with their standard deviation for three different sets of processors is shown at the bottom table inset for the standard PELE and the inversely proportional adaptive-PELE methods. (b–d) Analogous plots for PR, B-GPRC, and A-GPCR. A complete list of all data is shown in Supplementary Information.

Figure 4

PR binding mechanism. Two different views of the ligand entrance and the plasticity upon progesterone binding in PR. (a) Different ligand snapshots along the binding with two protein structures highlighting the initial closed (red cartoon) and intermediate open states (white cartoon). (b) A closer zoom at the entrance region with the ligand shown in the native bound structure; same color-coding as in the (a) panel but for the ligand (shown with atom element colors).

Figure 5

A-GPCR binding mechanism. (a) Different ligand snapshots showing the binding pathway from the initial structure (in red) to the bound pose (in blue), including Y148 and L225, which follow the same color-code. The white cartoon protein and the colored licorice ligand correspond to the bound crystal structure. (b) Side chain conformations for Y148 and L225, where the red licorice corresponds to the crystal structure. In grey lines, we show all the different conformations for those cluster centers along the adaptive process, and in colored licorice we show the resulting main conformations after a k-medoids clustering.

Figure 6

Induced-fit docking studies. (a) PR system: protein structure from PDB ID:1A28 and ligand structure from PDB ID:3KBA. (b) sHE system: protein structure from PDB ID:5AKE and ligand structure from PDB ID:5AM4. (c) sHE system: protein structure from PDB ID:5ALX and ligand structure from PDB ID:5AI5. In the upper panels we show the RMSD evolution along the simulation, in the middle ones the binding energy for the different RMSD values, and in the lower panels the native structure (atom-type colored), the lowest binding energy ligand structure (blue) and the starting ligand structure (red). Notice that in panel (b) the initial docking structure is slightly outside the active site (shown in the inset).