Cyclic Peptide Inhibitors of the Tsg101 UEV Protein Interactions Refined through Global Docking and Gaussian Accelerated Molecular Dynamics Simulations

Tsg101 UEV domain proteins are potential targets for virus infection therapy, especially for HIV and Ebola viruses. Peptides are key in curbing virus transmission, and cyclic peptides have a greater survival time than their linear peptides. To date, the accurate prediction of cyclic peptide-protein receptors binding conformations still is challenging because of high peptide flexibility. Here, a useful approach combined the global peptide docking, Gaussian accelerated molecular dynamics (GaMD), two-dimensional (2D) potential of mean force (PMF), normal molecular dynamics (cMD), and solvated interaction energy (SIE) techniques. Then we used this approach to investigate the binding conformations of UEV domain proteins with three cyclic peptides inhibitors. We reported the possible cyclic peptide-UEV domain protein binding conformations via 2D PMF free energy profiles and SIE free energy calculations. The residues Trp145, Tyr147, and Trp148 of the native cyclic peptide (CP1) indeed play essential roles in the cyclic peptides-UEV domain proteins interactions. Our findings might increase the accuracy of cyclic peptide-protein conformational prediction, which may facilitate cyclic peptide inhibitor design. Our approach is expected to further aid in addressing the challenges in cyclic peptide inhibitor design.


Introduction
Human Tsg101 proteins are important receptor proteins required for budding of the enveloped viruses (HIV, RSV, HSV-1/2, Ebola and parainfluenza) [1][2][3][4]; thus, protecting the release of virions from an infected cell is crucial [5,6]. Tsg101 proteins are recruited to the virus budding sites by binding to the PT/SAP motifs located within the HIV Gag p6 region and the Ebola Vp40 proteins, and then Tsg101 proteins and other cellular factors complete the budding process [2]. Inactivation of the Tsg101-expressing gene [5] through mutation or deletion results in the budding process failure of HIV

Native
Pro-Glu-Pro-Thr-Ala-Pro-Pro-Glu-Glu Peptide-protein docking remains a great challenge for most computational docking programs because of the peptide flexibility [13,14]. Small chemical molecules usually can strongly interact with small binding motifs and grasp the motifs in protein receptors. Except GPCR receptors, the peptides normally interact with the largest pockets and only bind to the protein receptors surface. Therefore, most peptides and protein receptors can't form stable complex structures [15]. There are template-based, local and global methods in regard to the peptide-protein docking simulations [14]. Since the template-based and local peptide docking methods require the 3D structural binding information, these approaches must be tested and verified by qualifying with the Critical Assessment of Prediction of Interactions criteria (CAPPRI) [16,17]. These approaches usually have a high performance but are often limited by the availability of experimental 3D structure information. The global peptide docking is that free peptides bind to protein receptors without experimental 3D structure information, but it is challenging to account for the protein receptors and peptides conformational flexibility. The ClusPro docking is a fast and global peptide docking method [18,19]. Unfortunately, this method still cannot provide more accurate 3D structural predictions of the peptide binding [20]. To address this issue, the approach combined ClusPro and Gaussian accelerated MD (GaMD) methods have been successfully tested [13]. The combined approach can refine the complex structures and significantly reduced the linear peptide backbone RMSDs to 0.6-2.7 Å. We also followed our strategy to simulate the refinement of the Tsg101 UEV protein and the linear peptide (amino sequence: PEPTAPPEE). Then our refinement was compared with the x-ray structure (PDB ID: 3obu). Our results are shown in Figures S1-3. In this study, we used molecular simulation techniques to investigate the interactions of UEV domain protein with these cyclic peptides (N to C terminal, Table 1) with the approach combined the ClusPro and GaMD methods. For verifying our predictions, the solvated interaction energy (SIE) free energy calculations were performed, and our predicted free energies were comparied with the experimental binding free energies (Table 1) [21]. Our simulations provide further insight into the binding interactions between UEV domain protein and cyclic peptides.

Materials and Methods
We followed Miao et al.'s approach, and the ClusPro docking, normal MD, GaMD, and binding free energy calculations were performed in the cyclic peptide-UEV domain protein receptor calculations [13]. This strategy consists of the following steps: Step 1: Building and optimizing the cyclic peptides by using the EVGAZZ software.
Step 2: Docking the cyclic peptides to UEV domain protein receptor by using ClusPro software.
Step 3: Selecting possible cyclic peptide-UEV domain protein receptor complex structures from step 2.
Step 7: Performing analysis of reaction coordinates by using AmberTools cpptraj software to obtain 2D potential of mean force profiles (2D PMF) with PyReweighting toolkit.
Step 8: The complex structures with the lowest potential of mean force (PMF) values were performed in 10-ns cMD simulations for SIE free energy calculations and binding mode analysis [22].

Cyclic Peptide Building
The three cyclic peptides (Table 1) were built on VEGAZZ. UEV domain protein crystal structures (PDB ID: 3obu) were selected as receptors for the cyclic peptide docking simulations with the ClusPro Dock protocol [22].

Cyclic Peptide-UEV Domain Protein Docking
The ClusPro Dock web-server were applied in the cyclic peptide-UEV domain receptors docking simulations [18,23]. Next, the possible cyclic peptide-UEV domain protein receptor complex structures were selected after sorting the weighted score values. The selected complex structures were then selected for subsequent equilibrium and production GaMD simulations.

The Detail Information of GaMD Simulations
From the ClusPro peptide docking simulations, the complex structures with the lowest potential of mean force (PMF) values were selected and then inserted into TIP3P water molecules box (box size: approximately 7.14 × 9.21 × 7.80 nm 3 ). These initial structures were chosen as refence structures for analyzing the C.M. distances and backbone RMSDs in our PMF profiles. These initial cyclic peptide-UEV domain protein structures were then simulated using the AMBER 18 software with the all-hydrogen amino acid AMBER FF99SB force field [24][25][26]. All initial cMD simulations were performed in the NPT ensembles with a constant temperature of 310 K, unless stated otherwise, using the Verlet integrator with an integration time step of 0.002 ps and SHAKE constraints [27] for all covalent bonds involving hydrogen atoms. For the long-range electrostatics, electrostatic interactions between molecules were performed using the particle mesh Ewald (PME) method [28], and the 2.0 nm cutoff distance was for van der Waals forces calculations. At first, these initial structures were minimized for 100,000 conjugate gradient steps and then performed with 1-ns NVT and 1-ns NPT cMD simulations. Additionally, the final structures from our cMD simulations were performed in 20-ns GaMD equilibration and 300-ns GaMD production simulations [29]. GaMD simulation trajectories were recorded every 0.2 ps for analyzing the C.M. distances and backbone RMSDs in our PMF profiles. Four times 300-ns GaMD production simulations were selected to analyz the C.M. distance and backbone RMSD of the cyclic peptides and UEV domain protein with the AmberTools cpptraj software [30]. For 2D PMF profiles, the PyReweighting toolkit was applied to reweight the four time GaMD each cyclic peptide-UEV domain protein systems [31]. The backbone RMSDs of UEV binding domain (residue numbers: 61-66, 90-100 and 139-145) and the distances (C.M) between the centers of the two domains (UEV domain protein residue numbers: 61-66, 90-100 and 139-145; the three cyclic peptides) ( Figure 1D) were used as reaction coordinates. From our 2D PMF profiles, the cyclic peptide-VEU domain protein complex structures with the lowest PMF values (PMF is equal to zero) were applied with 10-ns cMD simulations for SIE free energy calculations and binding mode (hydrogen bonds and nonbonding interactions) analysis by using AmberTools cpptraj software [30]. Detailed information concerning GaMD and SIE are listed in the supporting information.

Gaussian Accelerated Molecular Dynamics Simulation (GaMD)
GaMD is a biasing method for enhanced sampling of biomolecular conformations by adding a harmonic boost potential to smooth the biomolecular potential energy surface [29]. Where V is a bio molecular system potential energy, E is a systemic referenced energy, and ∆V is a harmonic boost potential energy. If V > E, ∆V = 0. If V < E, a ∆V is applied in the biosystem as follows: For GaMD simulations, the biased system potential energy (V*) is shown as where K is a harmonic force constant If V1 < V2, V1* < V2*. By replacing V* with Equation (2), the relationship is expressed as If V1 > V2. By replacing V* with Equation (2), the relationship is expressed as Combing Equations (3) and (4) is expressed as where Vmin and Vmax are the minimum and maximum potential energies. From Equation (5), we can derive where K constant is defined as K0 is the magnitude of the applied boost potential. Standard deviation (SD) of ∆V must be sufficiently small to ensure accurate reweighting [31].
where the Vave and σV are the average and SD of the potential energies, respectively, and σ∆V is the SD of ∆V with σ0 as a user-specified upper limit for accurate reweighting of potential energies. The SDs of the total potential and dihedral potential boosts are 10 kcal/mol each GaMD simulations. If E = Vmax, we can derive Equation (5) and the result is expressed as According to Equations (6) and (7), K0 can be defined as If E = Vmin + 1/k, we can derive Equation (8) and the result is expressed as The GaMD boost potential (∆V) is given as follows: For our GaMD simulations, the magnitude K0 = 1.0, Vmax and Vmin were obtained through normal MD (cMD) simulations.

SIE Free Energy Calculations
The binding free energies of the UEV domain protein with the three cyclic peptide inhibitors were calculated for 10 ns cMD simulations. The SIE [32] function is expressed as where E c and E vdw are the intermolecular electrostatic and van der Waals interaction energies, respectively. ∆G R bind is the reaction field energy between the bound states and unbond state, and ∆G R bind is solved by the Poisson equation [33,34]. ∆MSA is the change in the molecular surface area upon binding [35]. Additionally, the other parameters are as follows: α = 0.1048, Din = 2.25, ρ = 1.1, γ = 0.0129 kcal/(mol Å2), and C = −2.89 kcal/mol [36].

Cyclic Peptide-VEU Domaim Protein Docking Simulations
These cyclic peptides were docked into the UEV domain protein structure. The initial cyclic peptide-UEV domain proteins complex structures were selected after sorting the weighted score values. Our cyclic peptide docking results were shown in Table 2 and Figure 2. From Table 2 and Figure 2, Our results showed that these cyclic peptides could bind to the motif of the UEV domain protein, and these initial complex structures were reasonable for further cMD and GaMD simulations.

Prediction of Cyclic Peptide Binding Conformations through 2D PMF Profiles
GaMD simulations has been prove to refine the cyclic peptide-protein receptor complex structures. The C.M. distances and backbone RMSDs of the four times 300 ns-GaMD simulation trajectories were analyzed using AmberTools cpptraj software. The GaMD simulations were reweighted using the PyReweighting toolkitto calculate the 2D PMF profiles (Figure 2 and Table 3) Because the structure with the lowest PMF values might be reasonable and possible, the chosen complex structures were selected for SIE free energy calculations, and confirmed with experimental free energies. For the CP1 simulation system, the lowest PMFs were at the UEV binding domain backbone RMSD of 1.0 and the C.M. distance of 5.0 Å in the GaMD simulations, respectively. For the CP2 simulation system, lowest PMFs were at the UEV binding domain backbone RMSD of 2.0 and the C.M. distance of 9.0 Å. For the CP3 simulation system, the lowest PMFs were at the UEV binding domain backbone RMSD of 2.0 and the C.M. distance of 7.0 Å. When the cyclic peptide docking simulations were compared with GaMD refined structures, significant conformational changes and the reduced PMF values were observed in UEV domain protein during cyclic peptide binding ( Figure 3 and Table 3). These GaMD refined structures were selected to use in 10-ns cMD simulations for SIE free energy calculations and binding mode analysis. From the 2D PMF profiles and Table 3, GaMD simulations significantly refined the three cyclic peptide-UEV domain docking conformations.

SIE Free Energy Calculations
The predicting and experimental binding free energies (∆G) of the three cyclic peptides with UEV domain protein are presented in Tables 4 and 5. For the three cyclic peptides, the predicted ∆G values are −6.05, −5.31 and −6.52 kcal/mole, respectively. For the predicting binding free energies erroneous analysis, the error of our predicting binding free energies (∆Gdiff) are −0.51, 0.314 and 0.020 kcal/mol, respectively. Therefore, the GaMD refined complex structures could present a satisfactory correlation between the predicting and experimental binding free energies of all three cyclic peptides. Moreover, GaMD refined complex structures could provide more accurate and reasonable predictions.

Analysis of the Binding Modes of UEV Domain Protein with the Three Cyclic Peptide Inhibitors
Crystallographic structures revealed the binding modes of the UEV domain proteins with the native linear peptide, and this binding mode is shown in Figure 4. In the UEV domain protein-native peptide complex, the native peptide residues Glu6, Pro7, Thr8 and Pro10 can form hydrogen bonds with the Asn69 and Ser143 residues of UEV domain protein, whereas the native peptide residues Pro5, Glu6, Pro7, Thr8, Ala9, Pro10, Pro11 and Glu13 have nonbonding interactions with the Asp34, Thr58, Tyr63, Arg64, Tyr68, Ile70, Thr92, Met95, Val141 and Phe142 residues of UEV domain protein.
For analyzing the binding modes of the three cyclic peptides-UEV domain proteins, the 10-ns cMD simulations (SIE free energies calculation trajectories) were applied in the analysis. For finding out the hot-spot residues, the binding mode analysis of our 10-ns cMD simulations, which are more than a fraction of 30%, are shown in Table 6. The binding mode analysis of the UEV domain protein with the three cyclic peptides through ClusPro Dock simulations are shown in Figure 5, Figure 6, Figure 7. In the UEV domain protein-CP1 simulations, the CP1 residues Gly144, Trp145, and Tyr147 can form hydrogen bonds with the Arg61, Tyr60, and Ser140 residues of UEV domain protein, whereas the CP1 residues Trp145, Tyr147, and Trp148 have nonbonding interactions with the Asn63, Pro142, Pro139, Pro136, and Ile67 residues of UEV domain protein. In the UEV domain protein-CP2 simulations, the CP2 residue Gly144 can form hydrogen bonds with the Tyr65 residues of UEV domain protein, whereas the CP2 residues Trp145, Ser143, and Trp148 have nonbonding interactions with the Tyr60, Phe139, Pro142, and Gly62 residues of UEV domain protein. In the UEV domain protein-CP3 simulations, the CP3 residues Tyr147 and Trp148 can form hydrogen bonds with the residues Tyr65, Ser140, and Val139 residues of UEV domain protein, whereas the CP3 residues Trp145, Tyr147, Ile146, and Trp148 have nonbonding interactions with the Arg61, Gly62, Phe139, Val58, Asn63, Ile67, Pro136, and Pro137 residues of UEV domain protein. Moreover, these binding modes analysis could reveal the essential residues of the three cyclic peptides in the cyclic peptides-UEV domain proteins interactions.

Discussion
Our validation of the Tsg101 UEV protein and the linear peptide (amino sequence: PEPTAPPEE and PDB ID: 3obu) is reasonable and shown in the supporting information. We are confident about our predictions. Cyclic peptides are a yet underexploited candidates for drug discovery. Compared to their linear counterparts, cyclic peptides have enzymatic resistance and lower 3D conformational freedom, which allow them to bind to protein receptors with higher binding affinity and specificity [37]. For speeding up the cyclic peptides drug discovery, many peptide-protein computational docking approaches are developed to address this issue. Unfortunately, most peptide docking approaches remain a great challenge for most computational docking programs because of the peptide flexibility. Even though these docking methods are fast, these methods still lack thermodynamic refinements and cannot provide accurate peptide binding conformational predictions. These peptide docking methods can't perform docking simulations with cyclic peptides and can only perform related simulations with linear peptides [16,38,39]. To overcome this issue, we used an approach combined global peptide docking and GaMD methods to investigate the cyclic peptides-protein receptors interactions. This approach has been confirmed to succeed in predicting the linear peptide-protein binding conformations [13]. Compared to their linear counterparts, cyclic peptides indeed reduce 3D conformational freedom [40]. Therefore, this approach is suitable for studying the interactions between the three cyclic peptides and UEV domain proteins.
GaMD refinement simulations revealed the 2D PMF profiles (Figure 2 and Table 3), and we determined the most possible and reasonable complex structures with the lowest PMF values. Next, the 10-ns normal MD simulations were applied for SIE free energy calculations and the cyclic peptides-VEU domain protein binding mode analysis. Our predicting binding free energies have been confirmed to succeed in the three cyclic peptides binding conformational predictions. Therefore, we believe that our predictions could provide a reasonable and accurate insight into the interactions of UEV domain proteins with the three cyclic peptides. Since the mot of peptides can bind to the largest pocket on protein receptor surfaces, the complex aspect of the peptide binding still remains challenging. When our predictions were compared with the UEV domain protein experimental information (Table 4,  Table 5, Table 6 and Figure 4, Figure 5, Figure 6, Figure 7), the three cyclic peptide binding modes were found to be very different from the native peptide-UEV domain protein binding information. The cyclic CP1 and CP3 peptides can prevent the binding of UEV domain protein to P6 motif, thereby stopping HIV and Ebola virus budding [41]. The cyclic CP2 peptides indeed have weaker binding ability with UEV domain proteins.
Since the single point residue mutations of peptides might have different 3D conformations and affect the binding affinities with protein receptors [42,43]. The two cyclic peptides (CP2 and CP3) are the single residue mutation of the CP1 peptide. Therefore, we analyzed the binding modes of the cyclic peptides with UEV domain protein and found out the hot-spot residues each the cyclic peptide. Regarding the CP1-UEV domain protein interactions, the residues Gly144, Trp145, and Tyr147 of CP1 can form hydrogen bonds with UEV domain proteins, and residues Trp145, Tyr147, and Trp148 of CP1 can interact with UEV domain proteins via the nonbonding interactions. Regarding the CP2-UEV domain protein interactions, the Gly144 residue of CP2 can form hydrogen bonds with UEV domain proteins, and Trp145, Ser143, and Trp148 residues of CP2 can interact with UEV domain proteins via the nonbonding interactions. Regarding the CP3-UEV domain proteins interactions, our results show that the Tyr147 and Trp148 residue of CP3 can form hydrogen bonds with UEV domain proteins, and Trp145, Tyr147, Ile146, and Trp148 residues of CP3 can interact with UEV domain proteins via the nonbonding interactions. Compared with the three cyclic peptides, the mutation of Tyr147 can reduce UEV domain protein binding abilities. To sum up, the residues Trp145, Tyr147, and Trp148 of the native cyclic peptide (CP1) indeed play essential roles in the cyclic peptides-UEV domain proteins interactions.

Conclusions
For speeding up the cyclic peptides drug discovery, an approach, which combins global docking, normal MD simulations, GaMD simulations, 2D free energy profiles (2D PMFs), and SIE free energy techniques, was applied in studying the interactions between the three cyclic peptides and UEV domain proteins. From the GaMD refinement simulations, our 2D PMF profiles found out the most possible cyclic peptide-UEV domain protein complex structures. Then our SIE binding free energy predictions were close to the experimental binding free energies, and these results confirmed that our cyclic peptide binding conformational predictions seemed reasonable and possible. Next, we analyzed the binding modes of the cyclic peptides with UEV domain protein from 10-ns cMD simulations each cyclic peptide. The hot spot residues Trp145, Tyr147, and Trp148 of the three cyclic peptides can influence UEV domain protein binding affinities. Our study not only can increase the accuracy of cyclic peptide-protein 3D conformational prediction, but this study also can help future cyclic peptide drug discovery. Advancements in computational methods and computing power are expected to further aid in addressing the challenges in cyclic peptide drug design.
Supplementary Materials: The following are available online at http://www.mdpi.com/2073-4360/12/10/2235/s1. We followed our strategy to simulate the refinement of the Tsg101 UEV protein and the linear peptide (amino sequence: PEPTAPPEE). Then our refinement was compared with the x-ray structure (PDB ID: 3obu). Our results are shown in Figures S1−3. Figure S1: The linear peptide (amino sequence: PEPTAPPEE) is build by VEGAZZ software, Figure