Improved Modeling of Peptide-Protein Binding Through Global Docking and Accelerated Molecular Dynamics Simulations

Introduction

Peptides mediate up to 40% of known protein-protein interactions in higher eukaryotes. Peptide binding plays a key role in cellular signaling, protein trafficking, immune response, and oncology (Petsalaki and Russell, 2008; Das et al., 2013). In addition, peptides have served as promising drug candidates with high specificity and relatively low toxicity (Ahrens et al., 2012; Fosgerau and Hoffmann, 2015; Kahler et al., 2018; Lee et al., 2019). The number of peptide-based drugs being marketed is increasing in recent years (Ahrens et al., 2012; Fosgerau and Hoffmann, 2015; Kahler et al., 2018; Lee et al., 2019). Therefore, understanding the molecular mechanism of peptide-protein interactions is important in both basic biology and applied medical research.

Rational design of peptide-derived drugs usually requires structural characterization of the peptide-protein complexes. X-ray crystallography and nuclear magnetic resonance (NMR) have been utilized to determine high-resolution structures of peptide-protein complexes. These structures are often deposited into the Protein Data Bank (PDB) and also collected in specific databases focused on peptide-protein complex structures, including the PeptiDB (London et al., 2010), PepX (Vanhee et al., 2010), and PepBind (Das et al., 2013). Particularly, PeptiDB is a set of 103 non-redundant protein-peptide structures extracted from the PDB. The peptides are mostly 5–15 residues long (London et al., 2010). PepX contains 1,431 non-redundant X-ray structures clustered based on the binding interfaces and backbone variations. There are 505 unique peptide-protein interfaces, including those for the major histocompatibility complex (MHC) (14%), thrombins (12%), α-ligand binding domains (8%), protein kinase A (5%), proteases and SH3 domains (Vanhee et al., 2010). The PepBind contains a comprehensive dataset of 3,100 available peptide-protein structures from the PDB, irrespective of the structure determination methods and similarity in their protein backbone. More than 40% of the structures in PepBind are involved in cell regulatory pathways, nearly 20% in the immune system and ~30% with protease or other hydrolase activities (Das et al., 2013). These databases have greatly facilitated structure-based modeling and drug design of peptide-protein interactions. However, the number of currently resolved structures is only a small fraction of the peptide-protein complexes, as limited by the difficulties and high cost of X-ray and NMR experiments.

Computational methods have been developed for predicting the peptide-protein complex structures. In this regard, modeling of peptide binding to proteins has been shown to be distinct from that of extensively studied protein-ligand binding and protein-protein interactions. Notably, small-molecule ligands are able to bind deeply buried sites in proteins, but peptides normally bind to the protein surface, especially in the largest pockets. On the other hand, protein partners usually have well-defined 3D structures before forming protein-protein complexes, despite possible conformational changes during association. In contrast, most peptides do not have stable structures before forming complexes with proteins (Petsalaki and Russell, 2008). The biggest and immediate challenge for modeling of peptide-protein binding is that peptide structures are not known a priori. Furthermore, peptide-mediated interactions are often transient. The affinity of peptide-protein interactions is typically weaker than that of protein-protein interactions, because of the smaller interface between peptides and their protein partners. Therefore, new and robust computational approaches are developed to address the above challenges in the modeling of peptide-protein binding.

Molecular docking has proven useful in predictions of peptide-protein complex conformations (Ciemny et al., 2018). The commonly used approaches include template-based docking such as GalaxyPepDock (Lee et al., 2015), local docking of peptides to pre-defined binding sites such as Rosetta FlexPepDock (Raveh et al., 2011), HADDOCK (Trellet et al., 2013), and MDockPep (Xu et al., 2018), and global docking of free peptide binding to proteins such as CABS-dock (Kurcinski et al., 2015), PIPER-FlexPepDock (Alam et al., 2017), and PeptiDock (Porter et al., 2017). The template-based docking is highly efficient, but often limited to the availability of templates (Lee et al., 2015). Local docking is able to generate good quality models that meet the Critical Assessment of PRediction of Interactions (CAPRI) criteria (Janin et al., 2003). However, it requires a priori knowledge of the peptide binding site on the protein surface. In comparison, global peptide docking provides sampling of peptide binding over the entire protein surface without the need for pre-defined binding sites, but it is challenging to account for the system flexibility. In this regard, ClusPro PeptiDock has been developed for docking of motifs (short sequences) of peptides, which are found to sample only a small ensemble of different conformations (Alam et al., 2017). Structural ensemble of a peptide motif is built by retrieving motif structures from PDB that are very similar to the peptide's bound conformation. A Fast-Fourier Transform (FFT) based docking is then used to quickly perform global rigid body docking of these fragments to the protein. PeptiDock is thus able to alleviate the peptide flexibility problem through ensemble docking of the peptide motifs. Nevertheless, it remains challenging to account for the high flexibility of the peptides. Overall, peptide docking often generates poor predictions that require further refinement to obtain CAPRI-quality models.

Molecular dynamics (MD) is a powerful technique that enables all-atom simulations of biomolecules. MD simulations are able to fully account for the flexibility of peptides and proteins during their binding (Knapp et al., 2015; Wan et al., 2015; Salmaso et al., 2017; Yadahalli et al., 2017; Kahler et al., 2018). MD has been used to refine binding poses of peptides in proteins in the pepATTRACT (De Vries et al., 2017) and AnchorDock (Ben-Shimon and Niv, 2015) docking protocols. However, it is challenging to sufficiently sample peptide-protein interactions through conventional MD (cMD) simulations, due to the slow dynamics and limited simulation timescales. Computational approaches that combine many cMD simulations provide improved sampling of peptide-protein interactions, including supervised MD (Salmaso et al., 2017) and weighted ensemble (Zwier et al., 2016). Notably, weighted ensemble of a total amount of ~120 μs MD simulations has been obtained to investigate binding of an intrinsically disordered p53 peptide to the MDM2 Protein (Zwier et al., 2016). The simulation predicted binding rate constant agrees very-well with the experiments. However, expensive computational resources would be needed for applications of cMD simulations in large-scale predictions of peptide-protein complex structures.

On the other hand, enhanced sampling MD methods have been developed to improve biomolecular simulations (Christen and Van Gunsteren, 2008; Gao et al., 2008; Liwo et al., 2008; Dellago and Bolhuis, 2009; Abrams and Bussi, 2014; Spiwok et al., 2015; Miao and Mccammon, 2016). Multi-ensemble Markov models (Paul et al., 2017), which combine cMD with Hamiltonian replica exchange enhanced sampling simulations, have been used to characterize peptide-protein binding and calculate kinetic rates of a nano-molar peptide inhibitor PMI to the MDM2 oncoprotein fragment (Paul et al., 2017). While cMD is able to simulate fast events such as peptide binding, enhanced sampling simulations can capture rare events such as peptide unbinding. The steered MD (Cuendet et al., 2011), temperature-accelerated MD (Lamothe and Malliavin, 2018) and MELD (Modeling by Employing Limited Data) using temperature and Hamiltonian replica exchange MD (Morrone et al., 2017) have also been applied to study peptide-protein binding. In comparison, more enhanced sampling methods have been applied in studies of protein-ligand binding and protein-protein interactions, including the umbrella sampling (Torrie and Valleau, 1977; Kastner, 2011; Rose et al., 2014), metadynamics (Laio and Parrinello, 2002; Alessandro and Francesco, 2008; Saleh et al., 2017a,b,c), adaptive biasing force (Darve and Pohorille, 2001; Darve et al., 2008), steered MD (Cuendet and Michielin, 2008; Gonzalez et al., 2011), replica exchange MD (Sugita and Okamoto, 1999; Okamoto, 2004), accelerated MD (aMD) (Hamelberg et al., 2004; Miao et al., 2015), and Gaussian accelerated MD (GaMD) (Miao et al., 2015; Miao and Mccammon, 2017, 2018; Pang et al., 2017). Overall, enhanced sampling simulations of peptide binding to proteins have been under explored. Peptide-protein binding shows distinct characteristics as described above and requires the development of improved enhanced sampling approaches.

Here, we present a prototype of a novel computational approach that combines global peptide docking using PeptiDock and all-atom enhanced sampling simulations using GaMD to model peptide-protein binding. Three model peptides have been selected from the PeptiDB database of non-redundant peptide-protein complex structures (London et al., 2010). They include peptide motifs “PAMPAR” (Peptide 1), “TIYAQV” (Peptide 2) and “RRRHPS” (Peptide 3), which bind to the SH3 domain, X-linked lymphoproliferative syndrome (XLP) protein SAP and human PIM1 kinase, respectively. Starting with the lowest RMSD conformation selected from top 10 models of PeptiDock, GaMD significantly refines the peptide-protein complex structures. Furthermore, the simulations provided important insights into the mechanism of peptide binding to target proteins at an atomistic level. Thus, PeptiDock+GaMD is a promising approach for exploring peptide-protein interactions.

Methods

A Computational Approach Combining PeptiDock and GaMD

A new computational approach was designed to predict peptide-protein complex structures by combining peptide docking with PeptiDock and all-atom enhanced sampling simulation with GaMD (Figure S1). Initial peptide-protein complex structures were obtained using the ClusPro PeptiDock server. The first step in the PeptiDock protocol is fragment search: the PDB database is searched for fragments containing the target peptide motif. The templates are clustered and an FFT-based rigid docking is applied to the cluster centroids. Top-scoring poses are clustered again and the centroids of the largest clusters are chosen as the final results (Porter et al., 2017). For the purpose of this study—to show the viability of the protocol—only one pose within top 10 models of PeptiDock, known to be near native, was selected for further refinement using GaMD simulations.

System Setup

Three model peptides were selected from the PeptiDB database of non-redundant peptide-protein complex structures (London et al., 2010). They included peptide motifs “PAMPAR” (Peptide 1), “TIYAQV” (Peptide 2) and “RRRHPS” (Peptide 3), which bind to the SH3 domain, XLP protein and human PIM1 kinase, respectively. The free X-ray structures of target proteins is 1OOT, 1D1Z and 2J2I, respectively. The corresponding bound structures are 1SSH, 1D4T (Poy et al., 1999) and 2C3I (Pogacic et al., 2007), respectively. The free X-ray structures of the target proteins were used in the peptide docking and GaMD simulation. Both capped/neutral and uncapped/zwitterion terminus models were investigated in the GaMD simulations. In the neutral terminus model, the N- and C-termini were capped with ACE and NHE, respectively.

Peptide Docking

The standard ClusPro PeptiDock protocol was used for all three systems. In the first step, receptor structures were specified: 1OOT chain A (Peptide 1), 1D1Z chain A (Peptide 2) and residues 125-305 of 2J2I chain B (Peptide 3). The next step was specifying motifs—the templates for searching fragments in PDB database. The motif was specified as subsequence of the peptide with one or more wildcard symbols. Wildcards could be of two forms: “X,” denoting any amino acid substitution, and “[…],” denoting substitution by any amino acid from the list. e.g., “[FT]” means that either Phe or Tyr can take this place. It is recommended to adjust the motif to yield between 100 and 1,000 hits, while preserving the essential features for binding. For the studied systems, the following motifs were used for fragment search: “PXMPXR” for Peptide 1 [107 hits, see Ref. Hou et al., 2012], “TI[YF]XX[VI]” for Peptide 2 [686 hits, see Ref. Poy et al., 1999] and “RXRHXS” for Peptide 3 [198 hits, see Ref. Bullock et al., 2005]. Since PDB contains bound structures of the studied systems, a number of PDB entries were explicitly excluded from template search, as listed in Table S4. The next steps were performed automatically by the server (Porter et al., 2017), being the same for all systems. The extracted fragments were changed to the target peptide sequence using backbone-dependent rotamer library (Dunbrack and Karplus, 1993). The extracted fragments (hits) were clustered using the greedy algorithm according to their pairwise root-mean-square deviation (RMSD), with 0.5 Å cluster radius. The centroids of top 25 clusters were docked to the receptor using rigid-body FFT docking (Kozakov et al., 2006), exhaustively sampling all possible mutual orientations of the receptor and ligand, and ranking them using a special scoring function with a mixture of physics-based and knowledge-based terms (Kozakov et al., 2006; Chuang et al., 2008). The top-scoring poses of each fragment were pooled together and clustered based on their pairwise RMSDs, with 3.5 Å cluster radius. The clusters were ranked according to their sizes (Kozakov et al., 2005). The centroids of ten largest clusters were subjected to energy minimization with a CHARMM19-based force field using the ABNR algorithm. To demonstrate the protocol, only the lowest RMSD conformation obtained from top 10 PeptiDock models of each peptide was selected for refinement using GaMD simulations. The ranks of docking poses with the lowest peptide backbone RMSDs used were 9, 5, and 10 for Peptides 1, 2 and 3, respectively. It is important to note that each of the top-10 docking models will be refined and scored in a full version of the protocol in further studies.

GaMD Enhanced Sampling Simulations

GaMD was applied to refine the peptide-protein complex structures. Complexes were solvated in explicit water using tleap from the AMBER 18 package (Case et al., 2018). The Na⁺ and Cl⁻ ions were added to neutralize the system charge. The AMBER ff14SB force field parameters (Maier et al., 2015) and TIP3P model (Jorgensen et al., 1983) were used for the proteins/peptides and water molecules, respectively. Each system was minimized using steepest descent for 50,000 steps and conjugate gradient for another 50,000 steps. After minimization, the system was heated from 0 to 310 K in 1 ns simulation by applying 1 kcal/(mol•Å²) harmonic position restraints to the protein and peptide heavy atoms with a constant number, volume and temperature (NVT) ensemble. Each system was further equilibrated using a constant number, pressure and temperature (NPT) ensemble at 1 atm and 310 K for 1 ns with same restraints as in the NVT run. Another 2 ns cMD simulations were performed to collect potential energy statistics (including the maximum, minimum, average, and standard deviation). Then 18 ns GaMD equilibration after applying the boost potential was performed. Finally, four independent 300 ns GaMD production simulations with randomized initial atomic velocities were performed on each peptide system. Simulation frames were saved every 0.2 ps for analysis. Snapshots of all four GaMD production simulations (1,200 ns in total) were combined for clustering to identify peptide binding conformations, for which the hierarchical agglomerative algorithm in CPPTRAJ (Roe and Cheatham, 2013) was applied. The cutoff was set to 3.5 Å for the peptide backbone RMSD to form a cluster. The PyReweighting toolkit (Miao et al., 2014) was applied to reweight four GaMD simulations combined and recover the original free energy or potential of mean force (PMF) profiles of each peptide-protein system. The RMSDs of the peptide and protein backbone were used as reaction coordinates. Detailed descriptions of GaMD theory and energetic reweighting were shown in Supplementary Material.

Results

Prediction of Peptide Binding Conformations Through Docking and GaMD Simulations

There were no significant conformational changes in the protein during binding of Peptides 1 and 3 (Figures 1A,C). In comparison, binding of Peptide 2 induced a large structural rearrangement of the loop involving residues 67–74 in the protein (Figure 1B). In addition, Peptide 3 is highly charged as its first three N-terminal residues in the sequence are all arginine. These features of Peptides 2 and 3 raised the difficulty in accurate prediction of their peptide-protein complex structures. Peptide docking with PeptiDock showed different levels of accuracy: RMSDs of the peptide backbone compared with the bound X-ray structures were 3.3, 3.5, and 4.8 Å for the three peptides, respectively (Figures 1A–C and Table 1). The first two were of acceptable quality predictions according to the CAPRI peptide docking criteria (Janin et al., 2003), and the third one was slightly above acceptability cutoff. It should be noted that our flexible protein-peptide docking protocol PIPER-FlexPepDock (Alam et al., 2017) mentioned above is successful in obtaining high-quality model only in the case of Peptide 1, whereas the other two cases are challenging due to either significant receptor flexibility (Peptide 2) or remoteness of rigid-body docking poses to the native conformation (Peptide 3).

FIGURE 1

Figure 1. Docking poses (red) of three peptide motifs obtained using PeptiDock are compared with X-ray structures (green): (A) Peptide 1 “PAMPAR”, (B) Peptide 2 “TIYAQV,” and (C) Peptide 3 “RRRHPS”; Binding poses (red) of three model peptides obtained using the “PeptiDock+GaMD” are compared with X-ray structures (green): (D) Peptide 1, (E) Peptide 2, and (F) Peptide 3.

TABLE 1

Table 1. Comparison of 10 top-ranked clusters of three model peptides using the PeptiDock+GaMD approach.

Next, GaMD simulations were performed to refine the docking models. Analysis of simulation trajectories showed that the GaMD simulations were able to effectively refine the peptide binding pose. For Peptides 1 and 2, RMSDs of the peptide backbone relative to the X-ray structures decreased to <1 Å during the GaMD simulations (Figures 2A,B). Peptide 1 bound tightly to the protein target site throughout the four GaMD simulations. Peptide 2 reached the native conformation within ~10, ~90, ~120, and ~170 ns in the four GaMD simulations and stayed tightly bound during the remainder of the simulations. In comparison, Peptide 3 exhibited higher fluctuations and sampled the near-native conformation transiently during the GaMD simulations (Figure 2C). Nevertheless, the minimum RMSDs of peptide backbone compared with X-ray structures were identified to be 0.20, 0.22, and 0.73 Å for the three peptides, respectively (Figures 2A–C).

FIGURE 2

Figure 2. Time courses of peptide backbone RMSD obtained from four 300 ns GaMD simulations on (A) Peptide 1, (B) Peptide 2, and (C) Peptide 3. 2D potential of mean force (PMF) regarding the peptide backbone RMSD and protein backbone RMSD for (D) Peptide 1, (E) Peptide 2, and (F) Peptide 3. The black stars indicate the initial binding poses obtained using PeptiDock.

Furthermore, GaMD simulation snapshots of the peptide conformations were clustered using the backbone RMSDs relative to the X-ray structures. This procedure was similar to analysis of the peptide docking poses. The 10 top-ranked clusters of peptide conformations with the lowest free energies were obtained. The 1st top-ranked cluster exhibited peptide backbone RMSDs of 0.94 and 0.61 Å for Peptides 1 and 2, respectively (Figures 1D–E and Table 1). For Peptide 3, the 3rd top-ranked cluster showed the smallest peptide backbone RMSD of 2.72 Å (Figure 1F and Table 1). According to the CAPRI criteria (Janin et al., 2003), structural predictions for Peptides 1 and 2 were of sub-angstrom high quality and medium quality for Peptide 3. Therefore, GaMD simulations significantly refined docking conformations of the three peptide-protein complex structures. The simulation predicted bound conformations of the peptides were in excellent agreement with experimental X-ray structures with 0.6–2.7 Å in the peptide backbone RMSDs. In comparison, docking poses of the three peptides obtained from PeptiDock showed RMSDs of 3.3–4.8 Å (Table 1).

Peptide Binding Mechanism Revealed From GaMD

Free energy profiles were calculated from the GaMD simulations using the protein and peptide backbone RMSDs relative to the bound X-ray structures as reaction coordinates. For Peptide 1, only one low-energy minimum was identified near the native bound state (Figure 2D). This was consistent with the clustering result that the peptide backbone RMSD of the 1st top-ranked cluster was only 0.9 Å.

For peptide 2, two low-energy minima were identified, corresponding to peptide backbone RMSDs of 0.5 and 4.2 Å, respectively (Figure 2E). As described above, the binding of Peptide 2 induced a significant conformational change in the protein loop of residues 67–74 (Figure 1B). Thus, the loop backbone RMSD and peptide backbone RMSD relative to the bound X-ray structure were also used as reaction coordinates to compute another two-dimensional free energy profile (Figure 3A). The protein loop was highly flexible, sampling a large conformational space. The loop backbone RMSD ranged from ~0.2–~8.0 Å. This loop sampled two low-energy conformations, including the “Open” (bound) (RMSD < 1 Å) and “Closed” (free) states (RMSD ~3–6 Å) (Figure 3). Compared to the “Open” state, the “Closed” loop moved closer to the core domain of protein (Figure 3B). GaMD simulations successfully captured the conformational change of this loop. The peptide and protein loop accommodated each other to form the final bound conformation (Figure 3), suggesting an “induced fit” mechanism.

FIGURE 3

Figure 3. (A) 2D PMF calculated for binding of Peptide 2 regarding RMSDs of the peptide backbone and protein loop (residues 67–74) relative to the X-ray structure (PDB: 1D4T). (B) Representative conformation of “Closed” state (blue) in compared with initial conformation from “PeptiDock” (red) and X-ray structure (green).

For Peptide 3, GaMD sampled a broad low-energy well, centered at the ~4.3 and ~1.0 Å RMSDs for the peptide and protein backbone relative to the bound X-ray structure (Figure 2F). Overall, this peptide-protein complex underwent high fluctuations, visiting a large conformational space. Nevertheless, GaMD simulations sampled the native binding pose of Peptide 3, for which the peptide backbone RMSD decreased to ~1 Å at ~60 ns and 160 ns during one of the GaMD production runs (Sim1) (Figure 2C). In contrast to binding of Peptide 2 that involved induced fit of the protein receptor, binding of Peptides 1 and 3 did not induce significant conformational change of the receptors.

Effects of the Terminal Residue Charges on Peptide Binding

In addition to the neutral terminus model as described above, we simulated another model of the three peptides with zwitterion terminal residues that were charged. Compared with the neutral terminus models, larger fluctuations were observed in the zwitterion terminus models of the three peptides (Figures S2–S4). For Peptides 2 and 3, their backbone RMSDs could reach large values of ~40 and ~20 Å, respectively. These results suggested that the peptides could dissociate from the initial near-native bound pose obtained from docking. Furthermore, 10 top-ranked clusters of peptide conformations with the lowest free energies were also calculated through structural clustering and energetic reweighting (see Methods for details). For Peptide 1, the 1st top-ranked cluster exhibited the smallest backbone RMSD of 1.22 Å relative to the X-ray structure (Figure S5A and Table S1). The 2nd top-ranked clusters exhibited the smallest backbone RMSDs of 0.62 and 3.88 Å for Peptides 2 and 3, respectively (Figures S5B,C and Tables S2–S3). In summary, peptides with zwitterion terminal residues underwent higher fluctuations and the simulation predicted bound conformations deviated more from the native X-ray structures compared with the neutral terminal models.

Improved Sampling Efficiency of GaMD Compared With Conventional MD

In addition to GaMD simulations, another set of cMD simulations of the same lengths were performed for comparison in their sampling efficiency to refine peptide binding conformations. The peptides were simulated with neutral terminal residues. Compared with GaMD, cMD needed typically longer simulation time to refine the binding mode of Peptide 1 (Figure S6A). The cMD mostly failed to refine binding poses of Peptides 2 and 3, for which RMSD decrease was not observed in 3 out of 4 cMD simulations of Peptide 2 (Figure S6B) and all 4 cMD simulations of Peptide 3 (Figure S6C). The 1st top-ranked cluster exhibited peptide backbone RMSD of 0.96 Å for Peptide 1 (Table S1). For Peptide 2, the 2nd top-ranked cluster showed the smallest peptide backbone RMSD of 2.79 Å, suggesting a medium-quality model similar to the docking pose was obtained (Table S2). For Peptide 3, the 6th top-ranked cluster showed the smallest peptide backbone RMSD of 4.68 Å, being closely similar to the PeptiDock result (Table S3). Therefore, cMD was significantly less efficient in refining docking poses of peptides compared with GaMD.

Discussion

We have demonstrated that GaMD can successfully refine PeptiDock docking poses, and thus established the possibility of PeptiDok+GaMD combination to predict peptide-protein complex structures and explore the peptide binding mechanism. Three peptides with different difficulty levels were selected as model systems. Peptide 1 was the easiest one as the peptide is rigid and there was no conformational change in the protein during peptide binding. Both Peptides 2 and 3 were challenging for predicting bound conformations accurately. The binding of Peptide 2 involved a significant structural rearrangement of the residue 67–74 loop in the protein. Peptide 3 with dense residue charges proved difficult for both docking and GaMD simulations. Nevertheless, the GaMD refinement achieved high quality models for both Peptides 1 and 2, and medium quality prediction for Peptide 3. This approach showed promise to be widely applicable for other peptide-protein binding systems.

It is difficult for the current docking programs to account for large conformational changes of proteins during peptide binding (Ciemny et al., 2018). Even in the flexible docking calculation, only movements of protein side chains are often taken into account. This raised a challenge in the modeling of Peptide 2. On the other hand, cMD simulations could account for flexibility of the peptide and protein and had been applied to refine docking poses of peptides in proteins (Ben-Shimon and Niv, 2015; De Vries et al., 2017). However, cMD could suffer from insufficient sampling and limited simulation timescales. Indeed, cMD is significantly less efficient in refining docking poses of the peptides compared with GaMD, especially for Peptides 2 and 3. Thus, the GaMD enhanced sampling method has been used in this study. Remarkably, GaMD effectively captured the loop movement of Peptide 2 (Figure 3) and greatly refined the peptide docking poses (Figure 1E). In addition, high-performance GaMD simulations were performed using AMBER 18 on the GPUs. With NVIDIA Pascal P100 GPU cards, each of the 300 ns GaMD simulations took about 38.1, 43.5, and 53.2 h for Peptides 1, 2 and 3, respectively.

In summary, PeptiDock+GaMD has been demonstrated on predicting the peptide-protein complex structures and revealing important insights into the mechanism of peptide binding to proteins, using three distinct peptides as model systems. In the future, all top-10 models of the ClusPro PeptiDock will be refined with GaMD and a larger number of protein-peptide systems will be evaluated systematically. Furthermore, the effects of different force fields (e.g., CHARMM36m) and solvent models (e.g., TIP4P, implicit solvent, etc.) (Kuzmanic et al., 2019) are to be further investigated. Since excellent performance was obtained using the CHARMM19-based force field in the previous study of protein-peptide docking with ClusPro PeptiDock (Porter et al., 2017), we continued to use it as implemented in the ClusPro PeptiDock server for docking calculations in the present study. For refinement of the docking poses with GaMD, because AMBER18 was applied for running the simulations, the widely used AMBER ff14SB force field was selected instead. Nevertheless, it might be better to use newer and the same force field for different stages of the modeling protocol, which will be tested in future studies. Development of novel protocols to increase the accuracy of peptide-protein structural prediction will facilitate peptide drug design. Advances in the computational methods and computing power are expected to help us to address these challenges.

Data Availability Statement

All datasets generated and analyzed for this study are included in the article/Supplementary Material.

Author Contributions

YM and DK designed research. JW and AA performed research. JW, AA, DK, and YM analyzed data and wrote the paper.

Funding

This work was supported in part by the National Institutes of Health (R01GM132572), National Science Foundation (AF 1816314, DBI 1759277), Binational Science Foundation Grant (2015207), American Heart Association (Award 17SDG33370094), and the startup funding in the College of Liberal Arts and Sciences at the University of Kansas (KU). This work used the KU Center for Research Computing Cluster and supercomputing resources with the allocation award TG-MCB180049 through the Extreme Science and Engineering Discovery Environment (XSEDE), which was supported by National Science Foundation grant number ACI-1548562, and project M2874 through the National Energy Research Scientific Computing Center (NERSC), which is a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmolb.2019.00112/full#supplementary-material

References

Abrams, C., and Bussi, G. (2014). Enhanced sampling in molecular dynamics using metadynamics, replica-exchange, and temperature-acceleration. Entropy 16, 163–199. doi: 10.3390/e16010163

CrossRef Full Text | Google Scholar

Ahrens, V. M., Bellmann-Sickert, K., and Beck-Sickinger, A. G. (2012). Peptides and peptide conjugates: therapeutics on the upward path. Future Med. Chem. 4, 1567–1586. doi: 10.4155/fmc.12.76

PubMed Abstract | CrossRef Full Text | Google Scholar

Alam, N., Goldstein, O., Xia, B., Porter, K. A., Kozakov, D., and Schueler-Furman, O. (2017). High-resolution global peptide-protein docking using fragments-based PIPER-FlexPepDock. PLoS Comput. Biol. 13:e1005905. doi: 10.1371/journal.pcbi.1005905

PubMed Abstract | CrossRef Full Text | Google Scholar

Alessandro, L., and Francesco, L. G. (2008). Metadynamics: a method to simulate rare events and reconstruct the free energy in biophysics, chemistry and material science. Rep. Prog. Phys. 71:126601. doi: 10.1088/0034-4885/71/12/126601

CrossRef Full Text | Google Scholar

Ben-Shimon, A., and Niv, M. Y. (2015). AnchorDock: blind and flexible anchor-driven peptide docking. Structure 23, 929–940. doi: 10.1016/j.str.2015.03.010

PubMed Abstract | CrossRef Full Text | Google Scholar

Bullock, A. N., Debreczeni, J., Amos, A. L., Knapp, S., and Turk, B. E. (2005). Structure and substrate specificity of the Pim-1 kinase. J. Biol. Chem. 280, 41675–41682. doi: 10.1074/jbc.M510711200

PubMed Abstract | CrossRef Full Text | Google Scholar

Case, D. A., Cerutti, D. S., Cheatham, T. E., Darden, T. A., Red, T. J., Giese, H., et al. (2018). AMBER 2018. San Francisco, CA: University of California.

Google Scholar

Christen, M., and Van Gunsteren, W. F. (2008). On searching in, sampling of, and dynamically moving through conformational space of biomolecular systems: a review. J. Comput. Chem. 29, 157–166. doi: 10.1002/jcc.20725

PubMed Abstract | CrossRef Full Text | Google Scholar

Chuang, G. Y., Kozakov, D., Brenke, R., Comeau, S. R., and Vajda, S. (2008). DARS (Decoys As the Reference State) potentials for protein-protein docking. Biophys. J. 95, 4217–4227. doi: 10.1529/biophysj.108.135814

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciemny, M., Kurcinski, M., Kamel, K., Kolinski, A., Alam, N., Schueler-Furman, O., et al. (2018). Protein-peptide docking: opportunities and challenges. Drug Discov. Today 23, 1530–1537. doi: 10.1016/j.drudis.2018.05.006

PubMed Abstract | CrossRef Full Text | Google Scholar

Cuendet, M. A., and Michielin, O. (2008). Protein-protein interaction investigated by steered molecular dynamics: the TCR-pMHC complex. Biophys. J. 95, 3575–3590. doi: 10.1529/biophysj.108.131383

PubMed Abstract | CrossRef Full Text | Google Scholar

Cuendet, M. A., Zoete, V., and Michielin, O. (2011). How T cell receptors interact with peptide-MHCs: a multiple steered molecular dynamics study. Proteins 79, 3007–3024. doi: 10.1002/prot.23104

PubMed Abstract | CrossRef Full Text | Google Scholar

Darve, E., and Pohorille, A. (2001). Calculating free energies using average force. J. Chem. Phys. 115, 9169–9183. doi: 10.1063/1.1410978

CrossRef Full Text | Google Scholar

Darve, E., Rodriguez-Gomez, D., and Pohorille, A. (2008). Adaptive biasing force method for scalar and vector free energy calculations. J. Chem. Phys. 128:144120. doi: 10.1063/1.2829861

PubMed Abstract | CrossRef Full Text | Google Scholar

Das, A. A., Sharma, O. P., Kumar, M. S., Krishna, R., and Mathur, P. P. (2013). PepBind: a comprehensive database and computational tool for analysis of protein-peptide interactions. Genom. Proteom. Bioinf. 11, 241–246. doi: 10.1016/j.gpb.2013.03.002

PubMed Abstract | CrossRef Full Text | Google Scholar

De Vries, S. J., Rey, J., Schindler, C. E. M., Zacharias, M., and Tuffery, P. (2017). The pepATTRACT web server for blind, large-scale peptide-protein docking. Nucleic Acids Res. 45, W361–W364. doi: 10.1093/nar/gkx335

PubMed Abstract | CrossRef Full Text | Google Scholar

Dellago, C., and Bolhuis, P. G. (2009). “Transition Path Sampling and Other Advanced Simulation Techniques for Rare Events,” in Advanced Computer Simulation Approaches for Soft Matter Sciences III, eds C. Holm and K. Kremer (Berlin, Heidelberg: Springer), 167–233. doi: 10.1007/978-3-540-87706-6_3

CrossRef Full Text | Google Scholar

Dunbrack, R. L. Jr., and Karplus, M. (1993). Backbone-dependent rotamer library for proteins. Application to side-chain prediction. J. Mol. Biol. 230, 543–574. doi: 10.1006/jmbi.1993.1170

PubMed Abstract | CrossRef Full Text | Google Scholar

Fosgerau, K., and Hoffmann, T. (2015). Peptide therapeutics: current status and future directions. Drug Discov. Today 20, 122–128. doi: 10.1016/j.drudis.2014.10.003

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y. Q., Yang, L. J., Fan, Y. B., and Shao, Q. (2008). Thermodynamics and kinetics simulations of multi-time-scale processes for complex systems. Int. Rev. Phys. Chem. 27, 201–227. doi: 10.1080/01442350801920334

CrossRef Full Text | Google Scholar

Gonzalez, A., Perez-Acle, T., Pardo, L., and Deupi, X. (2011). Molecular basis of ligand dissociation in beta-adrenergic receptors. PLoS ONE 6:e23815. doi: 10.1371/journal.pone.0023815

PubMed Abstract | CrossRef Full Text | Google Scholar

Hamelberg, D., Mongan, J., and Mccammon, J. A. (2004). Accelerated molecular dynamics: a promising and efficient simulation method for biomolecules. J. Chem. Phys. 120, 11919–11929. doi: 10.1063/1.1755656

PubMed Abstract | CrossRef Full Text | Google Scholar

Hou, T., Li, N., Li, Y., and Wang, W. (2012). Characterization of domain-peptide interaction interface: prediction of SH3 domain-mediated protein-protein interaction network in yeast by generic structure-based models. J. Proteome Res. 11, 2982–2995. doi: 10.1021/pr3000688

PubMed Abstract | CrossRef Full Text | Google Scholar

Janin, J., Henrick, K., Moult, J., Eyck, L. T., Sternberg, M. J., Vajda, S., et al. (2003). CAPRI: a critical Assessment of PRedicted interactions. Proteins 52, 2–9. doi: 10.1002/prot.10381

PubMed Abstract | CrossRef Full Text | Google Scholar

Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W., and Klein, M. L. (1983). Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 79, 926–935. doi: 10.1063/1.445869

CrossRef Full Text | Google Scholar

Kahler, U., Fuchs, J. E., Goettig, P., and Liedl, K. R. (2018). An unexpected switch in peptide binding mode: from simulation to substrate specificity. J. Biomol. Struct. Dyn. 36, 4072–4084. doi: 10.1080/07391102.2017.1407674

PubMed Abstract | CrossRef Full Text | Google Scholar

Kastner, J. (2011). Umbrella sampling. WIREs Mol. Sci. 1, 932–942. doi: 10.1002/wcms.66

CrossRef Full Text | Google Scholar

Knapp, B., Demharter, S., Esmaielbeiki, R., and Deane, C. M. (2015). Current status and future challenges in T-cell receptor/peptide/MHC molecular dynamics simulations. Brief Bioinform. 16, 1035–1044. doi: 10.1093/bib/bbv005

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozakov, D., Brenke, R., Comeau, S. R., and Vajda, S. (2006). PIPER: an FFT-based protein docking program with pairwise potentials. Proteins 65, 392–406. doi: 10.1002/prot.21117

PubMed Abstract | CrossRef Full Text | Google Scholar

Kozakov, D., Clodfelter, K. H., Vajda, S., and Camacho, C. J. (2005). Optimal clustering for detecting near-native conformations in protein docking. Biophys. J. 89, 867–875. doi: 10.1529/biophysj.104.058768

PubMed Abstract | CrossRef Full Text | Google Scholar

Kurcinski, M., Jamroz, M., Blaszczyk, M., Kolinski, A., and Kmiecik, S. (2015). CABS-dock web server for the flexible docking of peptides to proteins without prior knowledge of the binding site. Nucleic Acids Res. 43, W419–W424. doi: 10.1093/nar/gkv456

PubMed Abstract | CrossRef Full Text | Google Scholar

Kuzmanic, A., Pritchard, R. B., Hansen, D. F., and Gervasio, F. L. (2019). Importance of the force field choice in capturing functionally relevant dynamics in the von willebrand factor. J. Phys. Chem. Lett. 10, 1928–1934. doi: 10.1021/acs.jpclett.9b00517

PubMed Abstract | CrossRef Full Text | Google Scholar

Laio, A., and Parrinello, M. (2002). Escaping free-energy minima. Proc. Natl. Acad. Sci. U.S.A 99, 12562–12566. doi: 10.1073/pnas.202427399

PubMed Abstract | CrossRef Full Text | Google Scholar

Lamothe, G., and Malliavin, T. E. (2018). re-TAMD: exploring interactions between H3 peptide and YEATS domain using enhanced sampling. BMC Struct. Biol. 18:4. doi: 10.1186/s12900-018-0083-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, A. C., Harris, J. L., Khanna, K. K., and Hong, J. H. (2019). A comprehensive review on current advances in peptide drug development and design. Int. J. Mol. Sci. 20:2383. doi: 10.3390/ijms20102383

PubMed Abstract | CrossRef Full Text | Google Scholar

Lee, H., Heo, L., Lee, M. S., and Seok, C. (2015). GalaxyPepDock: a protein-peptide docking tool based on interaction similarity and energy optimization. Nucleic Acids Res. 43, W431–W435. doi: 10.1093/nar/gkv495

PubMed Abstract | CrossRef Full Text | Google Scholar

Liwo, A., Czaplewski, C., Oldziej, S., and Scheraga, H. A. (2008). Computational techniques for efficient conformational sampling of proteins. Curr. Opin. Struct. Biol. 18, 134–139. doi: 10.1016/j.sbi.2007.12.001

PubMed Abstract | CrossRef Full Text | Google Scholar

London, N., Movshovitz-Attias, D., and Schueler-Furman, O. (2010). The structural basis of peptide-protein binding strategies. Structure 18, 188–199. doi: 10.1016/j.str.2009.11.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Maier, J. A., Martinez, C., Kasavajhala, K., Wickstrom, L., Hauser, K. E., and Simmerling, C. (2015). ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB. J. Chem. Theory Comput. 11, 3696–3713. doi: 10.1021/acs.jctc.5b00255

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, Y., Feher, V. A., and Mccammon, J. A. (2015). Gaussian accelerated molecular dynamics: unconstrained enhanced sampling and free energy calculation. J. Chem. Theory Comput. 11, 3584–3595. doi: 10.1021/acs.jctc.5b00436

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, Y., and Mccammon, J. A. (2016). Unconstrained enhanced sampling for free energy calculations of biomolecules: a review. Mol. Simul. 42, 1046–1055. doi: 10.1080/08927022.2015.1121541

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, Y., and Mccammon, J. A. (2017). “Chapter six - gaussian accelerated molecular dynamics: theory, implementation, and applications,” in Annu. Rep. Comput. Chem., ed D. A. Dixon (San Diego, CA: Elsevier), 231–278. doi: 10.1016/bs.arcc.2017.06.005

CrossRef Full Text | Google Scholar

Miao, Y., and Mccammon, J. A. (2018). Mechanism of the G-protein mimetic nanobody binding to a muscarinic G-protein-coupled receptor. Proc. Natl. Acad. Sci. U.S.A 115, 3036–3041. doi: 10.1073/pnas.1800756115

PubMed Abstract | CrossRef Full Text | Google Scholar

Miao, Y., Sinko, W., Pierce, L., Bucher, D., Walker, R. C., and Mccammon, J. A. (2014). Improved reweighting of accelerated molecular dynamics simulations for free energy calculation. J. Chem. Theory Comput. 10, 2677–2689. doi: 10.1021/ct500090q

PubMed Abstract | CrossRef Full Text | Google Scholar

Morrone, J. A., Perez, A., Deng, Q., Ha, S. N., Holloway, M. K., Sawyer, T. K., et al. (2017). Molecular simulations identify binding poses and approximate affinities of stapled alpha-helical peptides to MDM2 and MDMX. J. Chem. Theory Comput. 13, 863–869. doi: 10.1021/acs.jctc.6b00978

PubMed Abstract | CrossRef Full Text | Google Scholar

Okamoto, Y. (2004). Generalized-ensemble algorithms: enhanced sampling techniques for Monte Carlo and molecular dynamics simulations. J. Mol. Graph. Model. 22, 425–439. doi: 10.1016/j.jmgm.2003.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Pang, Y. T., Miao, Y., Wang, Y., and Mccammon, J. A. (2017). Gaussian accelerated molecular dynamics in NAMD. J. Chem. Theory Comput. 13, 9–19. doi: 10.1021/acs.jctc.6b00931

PubMed Abstract | CrossRef Full Text | Google Scholar

Paul, F., Wehmeyer, C., Abualrous, E. T., Wu, H., Crabtree, M. D., Schoneberg, J., et al. (2017). Protein-peptide association kinetics beyond the seconds timescale from atomistic simulations. Nat. Commun. 8:1095. doi: 10.1038/s41467-017-01163-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Petsalaki, E., and Russell, R. B. (2008). Peptide-mediated interactions in biological systems: new discoveries and applications. Curr. Opin. Biotechnol. 19, 344–350. doi: 10.1016/j.copbio.2008.06.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Pogacic, V., Bullock, A. N., Fedorov, O., Filippakopoulos, P., Gasser, C., Biondi, A., et al. (2007). Structural analysis identifies imidazo[1,2-b]pyridazines as PIM kinase inhibitors with in vitro antileukemic activity. Cancer Res. 67, 6916–6924. doi: 10.1158/0008-5472.CAN-07-0320

PubMed Abstract | CrossRef Full Text | Google Scholar

Porter, K. A., Xia, B., Beglov, D., Bohnuud, T., Alam, N., Schueler-Furman, O., et al. (2017). ClusPro PeptiDock: efficient global docking of peptide recognition motifs using FFT. Bioinformatics 33, 3299–3301. doi: 10.1093/bioinformatics/btx216

PubMed Abstract | CrossRef Full Text | Google Scholar

Poy, F., Yaffe, M. B., Sayos, J., Saxena, K., Morra, M., Sumegi, J., et al. (1999). Crystal structures of the XLP protein SAP reveal a class of SH2 domains with extended, phosphotyrosine-independent sequence recognition. Mol. Cell 4, 555–561. doi: 10.1016/S1097-2765(00)80206-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Raveh, B., London, N., Zimmerman, L., and Schueler-Furman, O. (2011). Rosetta FlexPepDock ab-initio: simultaneous folding, docking and refinement of peptides onto their receptors. PLoS ONE 6:e18934. doi: 10.1371/journal.pone.0018934

PubMed Abstract | CrossRef Full Text | Google Scholar

Roe, D. R., and Cheatham, T. E. 3rd (2013). PTRAJ and CPPTRAJ: software for processing and analysis of molecular dynamics trajectory data. J. Chem. Theory Comput. 9, 3084–3095. doi: 10.1021/ct400341p

PubMed Abstract | CrossRef Full Text | Google Scholar

Rose, A. S., Elgeti, M., Zachariae, U., Grubmuller, H., Hofmann, K. P., Scheerer, P., et al. (2014). Position of transmembrane helix 6 determines receptor G protein coupling specificity. J. Am. Chem. Soc. 136, 11244–11247. doi: 10.1021/ja5055109

PubMed Abstract | CrossRef Full Text | Google Scholar

Saleh, N., Ibrahim, P., and Clark, T. (2017a). Differences between G-protein-stabilized agonist-GPCR complexes and their nanobody-stabilized equivalents. Angew. Chem. Int. Ed. Engl. 56, 9008–9012. doi: 10.1002/anie.201702468

PubMed Abstract | CrossRef Full Text | Google Scholar

Saleh, N., Ibrahim, P., Saladino, G., Gervasio, F. L., and Clark, T. (2017b). An efficient metadynamics-based protocol to model the binding affinity and the transition state ensemble of G-protein-coupled receptor ligands. J. Chem. Inf. Model. 57, 1210–1217. doi: 10.1021/acs.jcim.6b00772

PubMed Abstract | CrossRef Full Text | Google Scholar

Saleh, N., Saladino, G., Gervasio, F. L., and Clark, T. (2017c). Investigating allosteric effects on the functional dynamics of beta2-adrenergic ternary complexes with enhanced-sampling simulations. Chem. Sci. 8, 4019–4026. doi: 10.1039/C6SC04647A

PubMed Abstract | CrossRef Full Text | Google Scholar

Salmaso, V., Sturlese, M., Cuzzolin, A., and Moro, S. (2017). Exploring protein-peptide recognition pathways using a supervised molecular dynamics approach. Structure 25, 655–662.e652. doi: 10.1016/j.str.2017.02.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Spiwok, V., Sucur, Z., and Hosek, P. (2015). Enhanced sampling techniques in biomolecular simulations. Biotechnol. Adv. 33, 1130–1140. doi: 10.1016/j.biotechadv.2014.11.011

PubMed Abstract | CrossRef Full Text | Google Scholar

Sugita, Y., and Okamoto, Y. (1999). Replica-exchange molecular dynamics method for protein folding. Chem. Phys. Lett. 314, 141–151. doi: 10.1016/S0009-2614(99)01123-9

CrossRef Full Text | Google Scholar

Torrie, G. M., and Valleau, J. P. (1977). Nonphysical sampling distributions in Monte Carlo free-energy estimation: Umbrella sampling. J. Comput. Phys. 23, 187–199. doi: 10.1016/0021-9991(77)90121-8

CrossRef Full Text | Google Scholar

Trellet, M., Melquiond, A. S., and Bonvin, A. M. (2013). A unified conformational selection and induced fit approach to protein-peptide docking. PLoS ONE 8:e58769. doi: 10.1371/journal.pone.0058769

PubMed Abstract | CrossRef Full Text | Google Scholar

Vanhee, P., Reumers, J., Stricher, F., Baeten, L., Serrano, L., Schymkowitz, J., et al. (2010). PepX: a structural database of non-redundant protein-peptide complexes. Nucleic Acids Res. 38, D545–D551. doi: 10.1093/nar/gkp893

PubMed Abstract | CrossRef Full Text | Google Scholar

Wan, S., Knapp, B., Wright, D. W., Deane, C. M., and Coveney, P. V. (2015). Rapid, precise, and reproducible prediction of peptide-MHC binding affinities from molecular dynamics that correlate well with experiment. J. Chem. Theory Comput. 11, 3346–3356. doi: 10.1021/acs.jctc.5b00179

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, X., Yan, C., and Zou, X. (2018). MDockPeP: an ab-initio protein-peptide docking server. J. Comput. Chem. 39, 2409–2413. doi: 10.1002/jcc.25555

PubMed Abstract | CrossRef Full Text | Google Scholar

Yadahalli, S., Li, J., Lane, D. P., Gosavi, S., and Verma, C. S. (2017). Characterizing the conformational landscape of MDM2-binding p53 peptides using molecular dynamics simulations. Sci. Rep. 7:15600. doi: 10.1038/s41598-017-15930-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Zwier, M. C., Pratt, A. J., Adelman, J. L., Kaus, J. W., Zuckerman, D. M., and Chong, L. T. (2016). Efficient atomistic simulation of pathways and calculation of rate constants for a protein-peptide binding process: application to the MDM2 protein and an intrinsically disordered p53 peptide. J. Phys. Chem. Lett. 7, 3440–3345. doi: 10.1021/acs.jpclett.6b01502

PubMed Abstract | CrossRef Full Text | Google Scholar