Search
Close this search box.

In silico design of peptide inhibitors for Dengue virus to treat Dengue virus-associated infections – Scientific Reports

Peptide inhibitor modeling and docking

The sequence of peptide inhibitors taken from the literature was uploaded to PepFold3 software and modeled. The PepFold3 server provided a total of five models of the peptide inhibitor. The model was then validated by PROCHECK server. Figure S1 shows the Ramachandran plot for the developed model and Fig. S2 shows the ERRAT plot for the model. According to Ramachandran plot analysis 76% of the residues were found in the most favored region, 24% of the residues were found in the additional allowed region while 0% of the residues were found in the disallowed region. ERRAT plot indicates an overall quality of the model as 100. The 3D structure of the peptide inhibitor was then docked with the envelope protein of the dengue virus. A total of 10 poses were generated. All the poses were ranked based on the docking score. The docking score of the complex was predicted as − 11.90 kcal/mol. The pose that revealed good docking score and interaction was then used for further study. Figure 1 shows the peptide inhibitor docked with the enveloped protein. The wild-type peptide inhibitors formed a total of seven hydrogen bonds with the envelope protein.

Figure 1
figure 1

The binding interface of the peptide inhibitor docked with the enveloped protein.

Candidate peptide inhibitor generation

In this work, the dAffinity (changes in binding affinities) of the original peptide and the new peptide derivatives were calculated using the residue scanning function of the MOE software. The reference peptide’s non-essential residues that have no role in interacting with the E protein were found by the alanine scanning approach and substituted with different residues using a process called residue scan mutagenesis. When mutated to alanine, several residues (I21, L14, G16, V17, T19, Q27, V28, and F29) have positive dAffinity values. In order to improve the proposed peptides’ binding affinity for the E protein, these residues were substituted with other residues. To create a peptide library, only residues that did not interact with E protein were selected for substitution. Some single-residue substitutions were found to have significantly negative dAffinity values when compared to the reference peptide. A native-to-specific mutation that has raised the binding affinity for the receptor is indicated by a negative dAffinity score. By mutating a few residues, the suggested peptides’ affinity for the envelope protein was improved (Table 1). By replacing I21 with Q21, L14 with H14, and V28 with K28 the binding affinity of the peptide inhibitors was increased. The newly designed peptide inhibitors with single residue mutation that improved binding affinity are listed in Table 1. To determine how these alterations affected complex stability, MD simulation was performed on the top four peptide complexes.

Table 1 Designed peptide inhibitors with reference to wild-type peptide.

Stability and estimated half-life of the peptide inhibitors

The ProtParam server (http://web.expasy.org/protparam/) was used to calculate the molecular weight (MW) and stability of the designed peptide inhibitors39. Table 2 shows the molecular weight, estimated half-life, and instability index of the newly designed peptide inhibitors. A protein whose instability index is smaller than 40 is predicted as stable40. All the peptide inhibitors showed an instability index of less than 40 indicating that all the peptide inhibitors were stable.

Table 2 Instability index analysis, molecular weight and estimated half-life analysis.

RMSD analysis

Root Mean Square Deviation (RMSD) are usually used to determine the stability of the receptor-ligand complex in MD simulations. In the present study, we performed explicit MD simulations of the most potent peptide inhibitors with their respective target, which is the envelope protein. The structural drift and stability were measured by plotting the root-mean-square deviation (RMSD) of protein Cα atoms as a function of time. All systems remained significantly stable over the course of the production run (Fig. 2). The RMSD plot in the WT peptide complex was stable from the beginning with a slight increase and decrease in the RMSD value but after 61 ns the RMSD reached its maximum value but then remained stable until the end of the simulation. In the I21Q complex, the RMSD plot was observed to be stable from the beginning and then remained unchanged for 50 ns. Between 50 and 60 ns, the RMSD increases slightly but then immediately decreases to its initial value for the entire simulation runtime (Fig. 2A). The L14H mutant complex RMSD plot shows a slightly higher RMSD value than the wild-type peptide complex but remained stable until the end of the MD run (Fig. 2B). The V28K mutant complex shows a static pattern of deviation as shown in Fig. 2C and remains stable throughout the simulation time. The V28S mutant complex showed that the RMSD plot was stable from 0 to 38 ns and then increased very slowly and slightly and the RMSd value reached from 3 to 9 Å (Fig. 2D). The observation is further supported by RMSF analysis which highlights the decreased mobility of active site residues in the case of peptide–bound complexes.

Figure 2
figure 2

RMSd plots for the wild type (Black), I21Q (Red), L14H (Blue), V28K (green), and V28S (purple) peptide-protein complexes. The RMSd value is present on the Y-axis.

Residual flexibility analysis

To examine peptide binding-mediated impacts on the structural flexibility of Envelope protein, root-mean fluctuations (RMSF) were computed by using the coordinates of the Cα atoms. The RMSF graph showed that the residues of the reference complex exhibit the highest instability at positions 1–25 (9.0 Å), 90–110 (15.0 Å), 120–150(18.5 Å), 210–260 (14.0 Å) and 270–310 (17.0 Å) in comparison to the mutated peptide complexes. However, I21Q mutant shows less fluctuations in MD simulation. The L14H mutant complex shows the highest fluctuation at positions 17–25 (23.0 Å) followed by a very slight increase and decrease in the RMSF plot. In the third mutant complex V28K, the RMSF graph showed no major fluctuations indicating structural stability. The V28S mutant complex gave the same result as L14H, initially showing a high peak at positions 15–19 followed by a slight fluctuation in the plot. In addition, the root mean square fluctuation (RMSF) analysis was performed based on residues indicated that all the designed peptides exhibited lower fluctuation in the binding sites compared to the reference peptide (Fig. 3).

Figure 3
figure 3

Root mean square deviation (RMSF) plots for the wild type (Black), I21Q (Red), L14H (Blue), V28K (green), and V28S (purple) peptide-protein complexes.

Hydrogen bonding interaction analysis from simulation trajectory

Stable hydrogen bonds are crucial for molecular recognition. In order to assess the stability of hydrogen bonds, it is necessary to employ molecular dynamics simulations, as comparing crystal structures alone is insufficient. Hydrogen bonding contact profile were analyzed from the simulation trajectory of all four designed peptide complexes and wild type complex of envelope protein (Fig. 4). According to hydrogen bond analysis, the H-bond connections between the peptide inhibitors we designed and the envelope protein are stronger than those between the wild type and the inhibitor. According to our research, generated peptides show a strong affinity for the envelope protein and may find application in therapy.

Figure 4
figure 4

Hydrogen bond (Hb) analysis plots for the wild type (Black), I21Q (Red), L14H (Blue), V28K (green), and V28S (purple) peptide-protein complexes.

Radius of Gyration

To understand the structural compactness, the radius of gyration (Rg) was calculated. The wild-type Rg value was initially 35.2 Å, but at 30 ns the Rg value gradually changed to 36.0 Å and then stabilized at the end of the simulation. The I21Q mutant shows an Rg value similar to that of the wild type from the beginning up to 10 ns, but after 10 ns the Rg value drops from 35.2 Å to 34.4 Å and shows a lower average Rg value than the wild-type complex. The L14H mutant complex shows a lower Rg graph than the wild type, but after 75 ns the value of Rg increases from that of the wild-type complex. A third mutant complex, V28K, experienced a slight increase and decrease in Rg between 5 and 50 ns and then exhibited the same pattern of Rg as the wild-type complex. The V28S mutant complex shows a small Rg value of about 33.6 Å from the beginning to 40 ns, but a large deviation was observed after 40 ns and the plot reached 35.2 Å and remained constant until the end of the simulation. The Rg of each designed and wild type complex was displayed in Fig. 5.

Figure 5
figure 5

Radius of gyration (Rg) plots for the wild type (Black), I21Q (Red), L14H (Blue), V28K (green), and V28S (purple) peptide-protein complexes. The Rg value is present on the Y-axis.

Analysis of conformational dynamics of top four Peptides

To understand the conformational dynamics of the reference peptide complex and the designed peptide complex, principal component analysis was performed using the AMBER v22 program. The fluctuation of the envelope protein is represented by its principal modes, which also indicate its direction (eigenvector) and magnitude (eigenvalue). PCA (Principal Component Analysis) contributes to characterizing the dynamic behavior of protein movement trajectories during MD simulations. Projection of top two PCA modes (PC1 and PC2) from reference peptide and designed peptides were plotted against 100 ns MD simulated conformation (Fig. 6). In order to define a unique conformation of the simulated system, PC1 and PC2 projections on the X and Y-axes, respectively, are essential subspaces that overlap with structurally similar conformations. It was observed that the protein clusters in the I21Q system, which were separated by the substrate from light blue to dark brown, had greater distribution and more local motions. The PC1 and PC2 regions are covered by dots in the I21Q system, extending from blue dots to dark brown dots with − 200 and + 200 and − 200 and + 100, respectively. The system L14H covers − 300 and + 200 along PC1 and − 200 and + 200 along PC2. The V28K PCA graph represents a region that lies between − 190 and + 200 along PC1 and between − 90 and + 190 along PC2.For the V28S system, the PCA plot covers − 450 and + 150 along PC1 and − 80 and + 200 along PC2.In contrast, regions − 250 and + 250 along PC1 and − 200 and + 180 along PC2 were covered by the wild-type system.

Figure 6
figure 6

Principal component analysis (PCA) plots for the I21Q (A), L14H (B), and WT complex (C), and V28K (D) and V28S (E), peptide-protein complexes.

Free energy landscape

To further comprehend the variation of the stabilized conformation of the protein–ligand during the 100 ns simulation, the free energy morphology was mapped. Free energy topography was plotted using data from the three dimensions of RMSD, Rg (Radius of gyration), and Gibbs free energies, indicated by color codes from higher (red) to lower (blue) energies, and two types of free energy topography were plotted in 3D for a more intuitive presentation. In addition, the minimum energy conformations were extracted. Each complex has multiple conformations at the lowest energy, as demonstrated in the rightmost column of Fig. 7. In this analysis, the PCA (principal component analysis) is utilized to assess the free energy landscape, providing a more reliable presentation of the time and energy-dependent conformation space of proteins. Since, free energy landscape approach effectively differentiates between the kinetic and thermodynamic characteristics of proteins, rendering it a highly precise tool for comprehensive assessment of protein dynamics. The FEL plots for PC1 and PC2 exhibit variations in Gibbs free energy ranging from black to dark brown. This indicates a transition of protein energy from an unstable state with high energy to a more stable one with reduced energy. The conformational dynamics of PC1 and PC2 of the wild-type peptide and the top four design peptides are shown by FEL analysis. When compared to the wild-type reference peptide, the design peptide complexes, I21Q, L14H, V28K, and V28S, achieved low energy conformations and more stable behavior. Therefore, compared to the wild-type peptide, the binding of these design peptides is more favorable due to changes in conformation and dynamics, as shown by the results.

Figure 7
figure 7

Free energy landscape (FEL) plots for the I21Q (A), L14H (B), and WT complex (C), and V28K (D) and V28S (E), peptide-protein complexes.

MMGBSA analysis

The MMGBSA analysis was carried out for all the protein-peptide complexes and the result was compared with the wild-type complex. The binding energy calculation revealed that among all the systems the V28S peptide revealed strong binding toward the receptor. The delta G value of the V28S was good (− 51.7962) as compared to all other systems. The delta G value for the mutant system V28K was found as (− 50.9924) followed by L14H. As compared to the wild-type system the binding energy was good for the mutated systems. Table 3 represents the binding energy for all the systems.

Table 3 MMGBSA analysis of all the complexes.