Computational investigations of potential inhibitors of monkeypox virus envelope protein E8 through molecular docking and molecular dynamics simulations

Active site residues of MPXV E8The probable ligand-binding sites of the MPXV E8 protein were predicted using the PrankWeb (https://prankweb.cz/) and CASTp server (http://sts.bioe.uic.edu/castp/index.html?201l) respectively. Multiple ligand-binding sites were predicted in the target envelope protein. Pocket 1 comprise of 18 amino acids with a pocket score of 20.63, a probability score of 0.845 and an average conservation score of 1.715. These residues included Thr39, Lys41, Leu42, Arg44, Ser64, Ser65, His67, Tyr69, Val94, Tyr104, Lys108, Ile116, Thr173, Ile174, Asn175, Ser177, Asn218, and Arg220. In pocket 2, we observed 12 active residues with an average conservation value of 1.342, a probability score of 0.108, and a pocket score of 3.12. These residues are Asn18, Arg20, Leu21, Phe56, Pro58, Asp111, Thr168, Tyr169, Leu170, Val181, Asp228, and Glu230, which were found to be actively involved in binding with the proposed drug punicalagin. While residues such as Phe47, Phe56, Pro58, Tyr61, Val62, Leu63, Asp112, Ile115, and Val181 were identified in pocket 3 with an overall pocket score of 1.20, 0.012 (probability score), and 1.818 (average conservation score). The majority of these residues were significant in binding with the control drug maraviroc (Fig. 3A, B). The different binding pose analysis of ligands with the target envelope protein have been represented in the Supplementary File 1 (Figure 1). Among the three binding pockets predicted; we opted for the binding sites with the lowest energy value computed through blind docking. The predicted energy value for each binding site of MPXV E8 has been shown in the Supplementary Table 1.Figure 3Active binding pockets of MPXV E8 protein predicted from (A) PrankWeb and (B) CASTp.Molecular docking analysisTo discover a new potential candidate for treating MPXV, we performed a structure-based virtual screening by molecular docking using AutoDock Vina implicated in the PyRx software20. To find a potential lead candidate, a total of 70 compounds (60 natural and 10 commercial) were virtually screened using the docking approach. During docking, the different conformations of ligands were generated and the binding energy profile of each conformer was found to be in the range of − 2.3 to − 9.1 kcal/mol. Considering the docking analysis, the stable and favorable interactions were observed in maraviroc and punicalagin respectively. The binding energy value of the commercial drug maraviroc was found to be − 7.4 kcal/mol. This binding is attributed to the formation of two H-bonds and other interatomic interactions formed between the protein–ligand complexes (Fig. 4). Similarly, the binding energy value of the MPXV E8-punicalagin complex was calculated to − 9.1 kcal/mol. This binding interaction is favored by the formation of four H-bonds and other non-covalent interactions that existed between the ligand and the protein’s active site residues during docking. The predicted binding energy values of all the 70 docked compounds have been represented in the Supplementary Table 2.Figure 4Binding mode analysis of ligands with the protein’s active site residues during docking (A) MPXV E8-Maraviroc complex (i) A 2D plot (ii) stick model and (iii) surface representations of the docked conformer. (B) Structural representations of the docked MPXV E8-punicalagin complex in (i) 2D plot (ii) stick model and (iii) surface.Protein–ligand interaction analysisAfter retrieving the best-docked complexes, we subjected them to a protein–ligand interaction profile analysis using the PLIP server https://plip-tool.biotec.tu-dresden.de/plip-web/plip/index and found that the control drug maraviroc forms hydrophobic interactions and two H-bond interactions during docking. No salt bridge formations were observed in the protein-maraviroc complex. The hydrophobic interactions were observed in residues Val62, Trp96, Lys98, Tyr104, and the interacting ligand atoms (2363, 2356, 2363, and 2357, 2376, and 2378) respectively ( Fig. 5A and Table 2). Only a few residues (Phe47, Leu63) were involved in H-bonding interactions with the ligand atoms within the cut-off distance of 3.5–5 Å. In contrast, the natural compound punicalagin forms two hydrophobic interactions and four H-bonds along with one salt bridge formation. The hydrophobic interactions take place between the ligand atoms (2393 and 2372) and the interacting residues such as Phe56 and Tyr232 over a distance of 3.5 Å. The H-bonds are formed between the ligand atoms and the active site residues including Leu170, Glu230, and Asp228 within the cut-off distance of 3 Å (Fig. 5B and Table 2). In addition, a salt bridge was formed between Arg20 and the interacting ligand atom (2404 and 2405) with a distance of 4.35 Å (As shown in the Table 2).Figure 5Protein–ligand interaction profile analysis using PLIP. (A). MPXV E8-Maraviroc interactions (B). The visual representation indicates hydrogen bonding interactions with a blue line, while hydrophobic interactions were illustrated in light grey. The ligand was depicted in a sphere form, with oxygen (O), nitrogen (N), and carbon (C) atoms represented by red, blue, and light orange colours, respectively.Table 2 Protein–ligand interaction profile analysis from PLIP.Oral toxicity predictions of ligandsThe oral toxicity of the selected commercial and natural compounds was evaluated based on the LD50 values using the ProTox-II server (https://tox-new.charite.de/protox_II/). The oral toxicity of the commercial drug maraviroc was found to be in the fourth class, and the recommended dosage value for this compound is 1681 mg/kg. Concerning natural compounds, the oral toxicity of punicalagin belonged to the fifth-class level of toxicity with LD50 of 5000 mg/kg.MD simulations analysisTo understand the structure and interaction dynamics of the envelope protein MPXV E8 with the proposed drug punicalagin, we performed 100 ns MD simulations in three different states using GROMACS 2023.2 software25, and it was compared with the control drug maraviroc by performing another three independent simulations. MD analysis involved the calculation of deviations in the backbone of the protein–ligand complex from its initial conformation which is root-mean-square deviations (RMSD), root-mean-square fluctuations (RMSF), radius of gyration (ROG), intermolecular hydrogen bonding analysis, binding pose analysis, and solvent accessibility.RMS deviationsThe root-mean-square deviation was used to measure deviations in the backbone atoms of the target protein upon binding to the commercial drug maraviroc and the phytocompound punicalagin35. RMS deviations were calculated by superimposing the MD trajectories with respect to the initial structures as a function of time. Figure 6 demonstrates deviations in the backbone dynamics of each trajectory of the protein–ligand complexes corresponding to 100 ns simulations. From the figure, it can be observed that in the MPXV E8-maraviroc complex, the RMSD values of the docked complex fluctuate significantly with larger amplitudes over time. The complex exhibits intermittent deviations between 10–20 ns, 40–45 ns, and 75–100 ns, while it remains virtually stable for only a short period. While in the punicalagin bound state, the backbone RMSD of each trajectory increases steadily at the beginning of the simulations. After 28 ns, all the trajectories attain equilibrium and retain stable conformations till the simulation ends. The average RMS deviations of the protein–ligand complexes were calculated to be 0.78 ± 0.15 nm (MPXV E8-maraviroc) and 0.56 ± 0.07 nm (MPXV E8-punicalagin) respectively. Thus from the RMSD results, it can be suggested that our proposed drug punicalagin forms stable structures with the MPXV envelope protein E8 by lowering the RMSD when compared to the control drug maraviroc. Since each trajectory of the protein–ligand complexes exhibited structural convergence during the simulations. Therefore, only a single trajectory of the complexes was selected for further MD analysis. The average RMSD, RMSF, Rg, hydrogen bonding, SASA and PCA analysis for each replica of MPXV E8-maraviroc and MPXV E8-punicalagin complexes have been represented in the Supplementary File 1 (Supplementary Figures 2–7 and Supplementary Table 3).Figure 6Calculations of backbone RMS deviations for each trajectory of MPXV E8-Maraviroc and MPXV E8-punicalagin complexes from the initial position till 100 ns MD simulations. Colors: black, violet and cyan represent Traj1, Traj2 and Traj3 of MPXV E8-maraviroc complex; red, orange and green signify Traj1, Traj2 and Traj3 of MPXV E8-punicalagin complex.RMS fluctuationsRoot-mean-square fluctuation (RMSF) is a parameter used to determine the flexibility of each amino acid residue of the protein during simulations. RMSF describes changes in the structural conformations of the protein from the initial position till the end of the simulations based on the residue fluctuations36. To understand how structural deviations vary protein flexibility upon ligand binding, the RMSF was calculated from the 100 ns MD trajectory. Higher RMSF represents higher flexibility of the target protein during simulations. The RMS fluctuations of the envelope protein were found to be lower (i.e. 0.13 nm) when bound with punicalagin. The RMSF value increases with an average of 0.15 nm upon binding of maraviroc with the target protein during simulations (Supplementary Table 3). Figure 7a, b shows the location of the fluctuating residues in both the complex trajectories. From the figure, it can be observed that the MPXV E8-maraviroc complex (shaded in black) has a higher occurrence of fluctuating residues than the MPXV E8-punicalagin complex (shaded as red). The common residues found to be fluctuating in both complexes include Lys13, Glu30, Tyr85, His126, and Glu209. Besides, residues such as Leu57, His176, and Asn207 exhibited a higher rate of fluctuations in the MPXV E8-maraviroc complex. In the MPXV E8-punicalagin complex, the unique pattern of fluctuation was observed at the residue position Asp151. Overall, the highest RMSF peak was observed for the terminal end residues Asn207 (MPXV E8-maraviroc complex) and Glu209 (MPXV E8-punicalagin complex) with RMSF scores of 0.55 nm and 0.45 nm respectively. The residue on each peak was identified using the XMGRACE tool, and the detail of each identified residue number is provided in Supplementary File 2. Thus RMSF results suggest that most of the residues of MPXV E8 acquire higher stability and less fluctuation with punicalagin binding.Figure 7RMS fluctuations of the MPXV E8 residues with maraviroc and punicalagin corresponding to 100 ns MD simulations. (a). Higher fluctuations of residues in MPXV E8-maraviroc complex (black), and lesser residue fluctuations observed over the length of the protein in the punicalagin-bound state (red). (b) is the location of the peak residues in the MPXV E8 protein bound with maraviroc (black) and punicalagin (red) respectively.Radius of gyrationThe radius of gyration (ROG) is a parameter used to measure the equilibrium conformations of the total system, and this helps in determining protein compactness during simulations. The size or compactness of the proteins depends on their amino acid sequence composition; and it varies during protein–protein interactions, protein-nucleic acid interactions, protein–ligand binding, etc37. Thus to determine how changes in the residue flexibility affect protein compactness, we have calculated the ROG for both MPXV E8-maraviroc and MPXV E8-punicalagin complex trajectories (Fig. 8). We observed that when maraviroc binds to the target receptor protein, there is an increase in the Rg trajectory from the initial position till the end of the simulations. In this condition, the protein becomes less compact and unstable over different simulation timeframes. While in the MPXV E8-punicalagin trajectory, there is a gradual decrease in the Rg over different time scales of the simulation. This indicates that the protein structure is well-equilibrated, compact, and more stable after punicalagin binding. The lowest Rg value in the MPXV E8-maraviroc trajectory was found to be 1.72 nm and the maximum Rg value is 1.79 nm. The lowest and the highest Rg values of the MPXV E8-punicalagin complex were calculated to be 1.70 nm and 1.76 nm respectively. The average Rg values of MPXV E8-maraviroc and MPXV E8-punicalagin complexes were estimated to be 1.76 nm and 1.74 nm respectively (Supplementaery Table 3). Further to check the statistical significance of the Rg results, the p-values were calculated using the Kolmogorov–Smirnov test, and it was found to be < 0.05, suggesting that the difference between the two Rg groups is statistically significant.Figure 8The radius of gyration measuring compactness of the envelope protein during 100 ns MD simulations. Higher gyrations were observed from the initial conformations till the end of the simulations in the MPXV E8-maraviroc complex (black). The MPXV E8-punicalagin complex maintains more compact and stable conformations by decreasing gyrations over time (Red).Intermolecular hydrogen bonding analysis between protein–ligand complexesThe binding affinity and stability of ligands with the target receptor protein were analyzed by measuring the number of intermolecular H-bonds formed during the dynamic simulations. Figure 9 demonstrates the intermolecular H-bonding analysis of maraviroc and punicalagin with the target MPXV E8 over 100 ns MD simulations. The intermolecular H-bonding was examined using a cut-off value of 0.35 nm. The MPXV E8-maraviroc complex exhibited the formation of 0–4 H-bonds at an average of 0.21 bonds. The average H-bonding pairs formed between the commercial drug and the target protein was found to be 0.35 bonds. While on the other hand, there are 0–14 H-bonds formed between MPXV E8 and the natural compound punicalagin over different time scales of the simulation. This complex was stabilized by the continuous formation of H-bonds at an average of 6.78 bonds and the H-bonding pairs of 11.38 (Supplementary Table 3), which is higher when compared to the number of H-bonds formed between MPXV E8 and the commercial drug. This finding strongly supported the H-bonding interactions predicted between the ligand and the active site residues of MPXV E8 through docking.Figure 9Intermolecular H-bonding interactions of the protein–ligand complexes corresponding to 100 ns MD simulations. Lesser number of H-bonds in MPXV E8-maraviroc complex binding (Black). A higher number of H-bonds stabilize the MPXV E8-punicalagin complex (red). H-bonding pairs within 0.35 nm.Binding poses of protein–ligand complexes at different time scales of simulationsDifferent poses of the protein–ligand complex were taken using visual molecular dynamics (VMD) software to observe the structural changes during the simulation. The snapshot were taken at every 20 ns time intervals of 100 ns MD simulations (i.e. at 0 ns, 20 ns, 40 ns, 60 ns, 80 ns and 100 ns). It can be observed in Fig. 10A that the movement of maraviroc (green) in the binding pocket of the envelope protein MPXV E8 is very high compared to punicalagin (yellow). Punicalagin sits (topmost region) in the pocket with very little movement from its initial location (Fig. 10B), suggesting a relatively smaller RMSD. Further, the interacting amino acid residues were also analyzed at the different ns. LigPlot + v.2.2 software was used to generate the 2D plot of the interactions. Figure 11 shows the initial and final conformational changes in the protein–ligand complexes during the dynamic simulations. The interaction analysis shows that the MPXV E8-maraviroc complex has a very low number of interacting residues compared to the MPV E8-punicalagin complex (Fig. 11A, B). In the figure, it can be observed that at 0 ns, maraviroc was interacting with only one amino acid (Phe47) and, after 100 ns there were no closely interacting amino acid residues. In the case of the proposed drug punicalagin, the number of closely interacting amino acid residues was found to be very high. Initially (at 0 ns), punicalagin was found to interact with two amino acids (Leu170 and Phe56), while after the end of the MD simulation (at 100 ns), there was an increase in the number of interacting residues which were found to be seven (Arg20, Phe56, Asp111, Thr168, Leu170, Asp179, and Glu230) (Fig. 11C, D). This analysis confirms that punicalagin can exhibit better and more stable interactions with the target envelope protein when compared to the control drug maraviroc.Figure 10Binding pose analysis of protein–ligand complexes corresponding to 100 ns MD simulations (A). Structural changes in the MPXV E8-Maraviroc complex at 0 ns, 20 ns, 40 ns, 60 ns, 80 ns and 100 ns (B). Conformational transition upon punicalagin binding at 0 ns, 20 ns, 40 ns, 60 ns, 80 ns and 100 ns.Figure 11Interaction profile analysis of protein–ligand complexes at initial and end position of the simulations (A,B) Interactions of maraviroc with MPXV E8 at 0 and 100 ns (C,D) Interaction of punicalagin with MPXV E8 at 0 and 100 ns respectively.Analysis of solvent accessible surface area in MPXV E8Water plays an essential role in determining the structure, stability, dynamics, and function of proteins. Solvent accessible surface area (SASA) is a parameter that is used to measure the accessibility of protein residues to surrounding water molecules. Any changes in the accessibility might affect protein structure, dynamics, and functions38. In this study, SASA was used to measure the binding effect of ligands on the conformational behavior of MPXV E8. Figure 12 represents the solvent accessibility of MPXV E8 upon binding of maraviroc and punicalagin over 100 ns MD simulations.Figure 12SASA plot analysis corresponding to 100 ns MD simulations of MPXV E8-maraviroc and MPXV E8-punicalagin complex trajectory. Higher accessibility of MPXV E8 when treated with maraviroc (black). The accessibility decreases over a time scale of the simulation in the MPXV E8-punicalagin complex (red).By comparing the SASA plots, we found that both maraviroc and punicalagin exhibit similar trajectory patterns after binding with the target MPXV E8 protein. However, the solvent accessibility of MPXV E8 increases upon the binding of maraviroc. The SASA value of MPXV E8-maraviroc was found to be in the range between 117.97 and 128.27 nm2 with an average of 121.27 nm2 (Supplementary Table 3). On the other hand, the accessibility value of MPXV E8 gradually decreases upon binding of punicalagin over different simulation time frames. In the MPXV E8-punicalagin complex, the SASA value gradually increases from 112.79 to 132.55 nm2 at an average of 119.41 nm2. SASA results suggest punicalagin starts outside the binding pocket, gradually moving in, interacting with more amino acids, leading to tight binding and reduced solvent access area.Principal component analysisPCs are the eigenvectors of a covariance matrix, and the displacement of atoms along each eigenvector measure the concerted motion of proteins in conformational space. Because these motions help in determining protein structure and functions32,39. The PC analysis was computed for both maraviroc and punicalag in bound MPXV E8 protein corresponding to 100 ns MD simulations trajectory (Fig. 13). In the case of MPXV E8-maraviroc complex, the projection of the first two PCs (i.e. eigenvector 1 and eigenvector 2) exhibited widely spaced and extensive motions with an increase in the diagonalized covariance matrix value of 11.2506 nm2. While in the punicalagin bound state, protein motions were found to be well-clustered and more compact with a drop in the covariance matrix value of 4.24419 nm2.Figure 132D projections of the first two principal eigenvectors (PC1 and PC2) of MPXV E8-Maraviroc (black) and MPXV E8-Punicalagin (red) at 300 K corresponding to 100 ns MD simulations.MM-PBSA binding free energy calculationsTo characterize the strength of interaction of best docked ligands with their targets, the binding free energies of the docked complexes were computed using the MM-PBSA approach 35 considering a total of 1000 frames from the last 10 ns of MD simulation trajectories. The binding energies were decomposed into their energy components such as van der Waals energy, electrostatic energy, polar solvation energy and SASA non-polar solvation to gain insights into protein–ligand binding interactions. The binding free energy of MPXV E8-maraviroc complex was calculated to − 112.85 kJ/mol with − 181.49 kJ/mol (van der Waals energy), − 23.92 kJ/mol (electrostatic energy), − 119.65 (polar solvation), and − 27.09 kJ/mol (SASA/non-polar solvation). In the MPXV E8-punicalagin complex, the van der Waals energy contributed − 241.87 kJ/mol, − 54.39 kJ/mol (electrostatic interactions), 159.04 kJ/mol (polar solvation energy), and − 33.82 kJ/mol (non-polar solvation). The binding energy decomposition of the E8-punicalagin complex was estimated to − 170.68 kJ/mol. Table 3 summarizes the binding free energy decomposition values of the protein–ligand complexes corresponding to 100 ns MD simulations.
Table 3 MM-PBSA binding energy decomposition of the protein–ligand complexes computed from 100 ns MD simulations.Residue wise contribution analysis of the docked protein–ligand complexesThe binding energies of the docked complexes were decomposed to explore residue wise contribution towards the binding energy (Fig. 14). In the MPXV E8-maraviroc complex, the residues dominantly contributing to the binding energy include Phe47 (− 10.01 kJ/mol), Lys98 (− 12.42 kJ/mol) and Tyr104 (− 8.26 kJ/mol) while in the MPXV E8-punicalagin complex, residues such as Arg20 (− 18.01 kJ/mol), Phe56 (− 12.07 kJ/mol), Leu170 (− 11.10 kJ/mol), Glu230 (− 10.16 kJ/mol) and Tyr232 (− 13.01 kJ/mol) contributed significantly towards binding energy. The average decompositions of the ligand-binding residues were computed to − 8.79 kJ/mol (E8-maraviroc) and − 12.09 kJ/mol (E8-punicalagin) respectively. This firmly backed up the docking analysis suggesting the greater binding affinity and stability of the punicalagin with the target protein.Figure 14The binding free energy decomposition per residue of the protein–ligand complexes during 100 ns MD simulations. Residues with the lowest binding free energy profiles are highlighted.