3.1 Oral toxicity predictions of ligands
The 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 commercial drugs such as acyclovir, tecovirimat, and zanamivir was predicted to be in the fifth class with the recommended dosage values of 5000 mg/kg, 2028 mg/kg, and 5000 mg/kg respectively. Drugs like brincidofovir and sofosbuvir exhibited sixth-class toxicity with a dosage value of 5010 mg/kg and 12000 mg/kg respectively. The oral toxicity of cidofovir, efavirenz, maraviroc, remdesivir, and telaprevir was found to be in the fourth class, and the recommended dosage values for these drugs are 1681 mg/kg, 600 mg/kg, 1000 mg/kg, 1000 mg/kg, and 400 mg/kg respectively. Concerning natural compounds, the oral toxicity of artemisinin B and bergenin were found to be in the sixth class with LD50 values of 9000 and 10000 mg/kg respectively. While baicalein, cycloolivil, isoquercetin, isorhamnetin, and punicalagin belonged to the fifth-class level of toxicity with LD50 of 5000 mg/kg. The oral toxicity of chebulagic acid and swertiachoside B was found to be in the fourth class, while glucoevatromonoside exhibited a higher level of toxicity. Table 2 represents the predicted oral toxicity of commercial and natural compounds from the ProTox-II server.
3.2 Active site residues of MPXV E8
The probable ligand-binding sites of the MPXV E8 protein were predicted using the PrankWeb tool (https://prankweb.cz/). Among the different ligand-binding pockets, we selected pocket 2 with a pocket score of 3.12 and a probability score of 0.108. The number of active amino acids surrounding the binding pocket of MPXV E8 was found to be 11 with the recorded conservation score of 1.342. The predicted amino acids include both polar (Asn18, Arg20, Thr168, Tyr169, Asp228, Glu230, Tyr322) and non-polar (Phe56, Pro58, Leu170, Val181) residues corresponding to ligand binding.
3.3 Molecular docking analysis
To 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 software [31]. To find a potential lead candidate, a total of 20 compounds (10 each of natural and 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 -5.1 to -9.1 kcal/mol. Considering the docking analysis of commercialized drugs, the lowest binding energy was observed between the MPXV E8-maraviroc complex with the energy value of -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 (Figure 2).
By comparing the docking results, the most stable and favorable binding interactions were observed between punicalagin and the target envelope protein with the lowest binding energy of -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 20 docked compounds have been represented in Table 2.
3.4 Protein-ligand interaction analysis
After 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 (Table). 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Å. The higher the number of H-bonds formations, the greater the stability of the protein-ligand complex during docking. In addition, a salt bridge was formed between Arg20 and the interacting ligand atom (2404 and 2405) with a distance of 4.35 Å. The interaction profiles of the selected docked conformers are shown in Table 3.
3.5 MD simulations analysis
To understand the structure and interaction dynamics of the envelope protein MPXV E8 with the proposed drug punicalagin, we performed 100 ns MD simulations using GROMACS 2023.2 software [34], and it was compared with the control drug maraviroc by performing another simulation. 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.
3.5.1 RMS deviations
Root-mean-square deviation (RMSD) is used to measure deviations in the backbone atoms of proteins from the initial conformation till the end of the simulations[39]. In this study, the structural deviations in the monkeypox virus E8 protein in complex with both the ligands were determined by calculating the RMSD over 100 ns MD simulations and then plotted together for comparisons. The RMSD plot of the MPXV E8-punicalagin complex (Figure 3) attained equilibration around ~10 ns with ~0.7 nm deviations (red) on the y-axis, while the MPXV E8-maraviroc complex was found to have several significant intermittent deviations between 10-20 ns, 40-45 ns, and 75-100 ns. The average RMSD of the MPXV E8-punicalagin complex was calculated to be 0.62 nm, and for the MPXV E8-maraviroc complex, it was found to be 0.78 nm (Table 4). 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.
3.5.2 RMS fluctuations
Root-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 fluctuations [40]. 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 (Table 4). Figure 4a and 4b 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 1. Thus RMSF results suggest that most of the residues of MPXV E8 acquire higher stability and less fluctuation with punicalagin binding.
3.5.3 Radius of gyration
The 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, etc [41]. 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 (Figure 5). 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.73 nm respectively (Table 4). 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.
3.5.4 Intermolecular hydrogen bonding analysis between protein-ligand complexes
The 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 6 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 a total of 0-2 H-bonds with an average of 0.04 bonds. While there are 0-8 H-bonds formation in the MPXV E8-punicalagin complex. This complex was stabilized by the continuous formation of 2-4 bonds at an average of 3.1 bonds (Table 4), which is higher when compared to the number of bonds formed between MPXV E8-maraviroc complexes. This result was in good accord with the H-bonding interactions predicted between the ligand and the active site residues of MPXV E8 through docking.
3.5.5 Binding poses of protein-ligand complexes at different time scales of simulations
Different poses of the protein-ligand complex were taken using visual molecular dynamics (VMD) software to observe the structural changes during the simulation. The snapshot was taken at every 20 ns time interval until 100 ns. It can be observed in Figure 7A 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 (Figure 7B), 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 8 shows the initial and final conformational changes in the protein-ligand complexes during the dynamic simulations. Other intermediate timestep interaction plots are provided in Supplementary File 2. 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 (Figure 8A-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 0ns), 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) (Figure 8C-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.
3.5.6 Analysis of solvent accessible surface area in MPXV E8
Water 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 functions [42]. In this study, SASA was used to measure the binding effect of ligands on the conformational behavior of MPXV E8. Figure 9 represents the solvent accessibility of MPXV E8 upon binding of maraviroc and punicalagin over 100 ns MD simulations.
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 to 128.27 nm2 with an average of 122.38 nm2 (Table 4). 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 120.64 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.