3.2. Molecular dynamics (MD) simulation 3.2.1. Thermodynamic stability and flexibility analysis The 150 ns production simulations carried out for nine systems (complex of remdesivir, EGCG, TF1, TF2a, TF2b, TF3, hesperidin, myricetin, and quercetagetin with the SARS-CoV-2 RdRp) were stable on the basis of the potential energy and total energy (data not shown) of those complexes. Subsequently, the root-mean-square deviations (RMSDs) of backbone atoms relative to their respective initial positions were calculated for each complex and are shown in Figure 3(A). It is evident from Figure 3(A) that all the nine studied systems drifted from their initial positions during the first 50 ns, and after that, they reached equilibrium. The average RMSD values were estimated to be 2.30 Å, 2.45 Å, 1.87 Å, 2.28 Å, 1.68 Å, 2.47 Å, 1.90 Å, 2.03 Å and 1.88 Å for RdRp/remdesivir, RdRp/EGCG, RdRp/TF3, RdRp/TF2b, RdRp/TF2a, RdRp/myricetin, RdRp/quercetagetin, RdRp/hesperidin, and RdRp/TF1 complexes, respectively (Table 3). The least deviation was observed for RdRp/TF2a, while RdRp/myricetin displayed the highest deviation. We also investigated structural variations in the binding site, including all amino acids that fall within a radius of 5 Å from the inhibitor, and the same trend was observed (Figure S2A in the Supplementary Information). Overall, this suggests the convergence of our simulations. Figure 3. (A) Time evolution of root-mean-square deviations (RMSDs) of backbone atoms and (B) the root-mean-square fluctuations (RMSFs) of Cα atoms of nine complexes relative to their respective energy minimized structure. Table 3. The average backbone RMSD, radius of gyration (RoG), and solvent accessible surface area (SASA) for all nine complexes. The data are reported as average ± standard error of the mean (SEM). System RMSD (Å) RoG (Å) SASA (Å2) RdRp/Remdesivir 2.30 ± 0.03 29.96 ± 0.02 34973.20 ± 91.56 RdRp/EGCG 2.45 ± 0.05 29.52 ± 0.06 35026.03 ± 63.52 RdRp/TF3 1.87 ± 0.02 29.60 ± 0.01 34080.16 ± 53.41 RdRp/TF2b 2.28 ± 0.01 29.86 ± 0.02 35462.92 ± 50.40 RdRp/TF2a 1.68 ± 0.02 29.75 ± 0.02 34312.55 ± 112.89 RdRp/Myricetin 2.47 ± 0.03 29.88 ± 0.01 35395.35 ± 104.67 RdRp/Quercetagetin 1.90 ± 0.03 29.84 ± 0.01 34618.65 ± 51.25 RdRp/Hesperidin 2.03 ± 0.04 29.74 ± 0.02 34554.08 ± 47.25 RdRp/TF1 1.88 ± 0.03 29.86 ± 0.01 34420.69 ± 56.90 Next, we investigated the structural stability of remdesivir and eight polyphenols by estimating the temporal RMSDs of heavy atoms relative to their respective initial coordinates (see Figure S2B in the Supplementary Information). It is evident from Figure S2B that EGCG, myricetin, quercetagetin, TF1, TF3 displayed a rigid behavior in the bound form with an average RMSD of < 1.0 Å. However, remdesivir, TF2a, TF2b and hesperidin showed higher fluctuations as compared to the abovementioned polyphenols, and an average RMSD of > 2.0 Å was noted. To identify the regions which are flexible, the root-mean-square fluctuations (RMSFs) of Cα atoms of each residue are calculated and shown in Figure 3(B). From this analysis, we can get a better insight into what extent the binding of remdesivir and natural polyphenols affects the residual flexibility of RdRp (mainly Nsp12). Figure 3(B) indicates that RdRp/remdesivir and RdRp-polyphenols shared a similar RMSF pattern. Notable dynamic fluctuations were located in the non-active site domain, including both N-terminals and C-terminal. Regions around Asn150, Asp260, Arg305, Asn360, and Phe440 are found to be more flexible compared to the other area for all complexes. The binding pocket residues, such as Asp452, Lys545, Lys551, Tyr455, Arg553, Ala554, Arg555, Thr556, Asp618, Tyr619, pro620, Lys621, Cys622, Asp623 Arg624, Asn691, Asp760, Asp761, Phe793, Met794, Ser795, Lys798, Trp800, Glu811, Phe812, and Ser814 exhibited considerably low fluctuations for all the RdRp-inhibitor complexes. In the case of RdRp/TF3 and RdRp/TF2a, the binding site residues displayed lesser fluctuations compared to the other RdRp-polyphenol complexes. This suggests that TF3 and TF2a are likely to be bound to RdRp more strongly than the other polyphenols. Since the radius of gyration (RoG) helps us to understand the protein structural compactness, RoG of each complex was monitored and represented in Figure S3A in the Supplementary Information. The average values of RoG are 29.96 Å, 29.52 Å, 29.60 Å, 29.86 Å, 29.75 Å, 29.88 Å, 29.84 Å, 29.74 Å and 29.86 Å for RdRp complexed with remdesivir, EGCG, TF3, TF2b, TF2a, myricetin, quercetagetin, hesperidin and TF1 respectively (Table 3). This suggests that the structural compactness remained unchanged during simulations. Finally, the solvent-accessible surface area (SASA) was also explored, and the time evolution of SASA for four RdRp-polyphenol complexes are shown in Figure S3B in the Supplementary Information. The average values of SASA are reported in Table 3. Binding of an inhibitor to the substrate changes SASA and sometimes could greatly affect the protein structure. Here, a relatively higher SASA value was obtained for RdRp/TF2b (35462.9 Å2) compared to the other RdRp/inhibitor complexes. On the other hand, the lowest SASA value was noted for RdRp/TF3 (34080.2 Å2). Thus, it can be suggested that the binding of TF3 could potentially reduce protein expansion. 3.2.2. Binding free energy analysis We predicted the binding free energy of all nine complexes by utilizing the MM-PBSA scheme, and four polyphenols, namely EGCG, TF3, TF2b, and TF2a, displayed a higher estimated affinity compared to remdesivir as depicted in Figure 4. Various components of the binding free energy of EGCG, TF3, TF2b, and TF2a are reported in Table 4. The remaining four polyphenols which showed lower estimated affinity compare to remdesivir are shown in Table S2 in Supplementary Information. It can be noted from Figure 4 that the intermolecular van der Waals (ΔEvdW) and electrostatic (ΔEelec) terms are favorable for the ligand binding, whereas the desolvation of polar groups (ΔGpol) opposes the complex formation. Non-polar solvation free energy (ΔGnp) is favorable to the binding for all cases. A similar trend was observed in our earlier study (Sk, Roy, Jonniya, et al., 2020). Figure 4. Energy components (kcal/mol) for the binding of remdesivir and four polyphenols to RdRp receptor. ΔEvdW, van der Waals interaction; ΔEele, electrostatic interaction in the gas phase; ΔGpol, polar solvation energy; ΔGnp, non-polar solvation energy, and ΔGbind, estimated binding affinity. Table 4. Energetic components of the binding free energy of RdRp and natural polyphenols along with remdesivir complexes in kcal/mol. Data are represented as average ± SEM. Components Remdesivir EGCG TF3 TF2b TF2a ΔEvdW −31.85 ± 0.15 −25.11 ± 0.18 −37.82 ± 0.21 −30.66 ± 0.23 −22.55 ± 0.19 ΔEelec −98.40 ± 0.70 −69.38 ± 0.73 −123.63 ± 0.88 −47.18 ± 0.64 −95.28 ± 1.27 ΔGpol 109.97 ± 0.57 71.62 ± 0.48 124.47 ± 0.58 55.01 ± 0.49 94.94 ± 1.10 ΔGnp −4.29 ± 0.01 −4.15 ± 0.01 −5.29 ± 0.01 −3.91 ± 0.01 −4.28 ± 0.02 aΔGsolv 105.68 ± 0.57 67.47 ± 0.48 119.18 ± 0.58 51.1 ± 0.49 90.66 ± 1.10 bΔGpol + elec 11.57 ± 0.90 2.24 ± 0.87 0.84 ± 1.05 7.83 ± 0.80 −0.34 ± 1.68 cΔEMM −130.25 ± 0.71 −94.49 ± 0.75 −161.45 ± 0.90 −77.84 ± 0.68 −117.83 ± 1.28 ΔGbindSim −24.57 ± 0.91 −27.02 ± 0.89 −42.27 ± 1.07 −26.74 ± 0.83 −27.17 ± 1.69 a ΔGsolv = ΔGnp + ΔGpol, b ΔGpol + elec = ΔEelec + ΔGpol, c ΔEMM = ΔEvdW + ΔEelec. It is evident from Table 4 that for all complexes, ΔEvdW varies between −22.55 kcal/mol and −37.82 kcal/mol while ΔEelec ranges from −47.18 kcal/mol to −123.63 kcal/mol. Furthermore, in the cases of RdRp/remdesivir, RdRp/EGCG, RdRp/TF3, and RdRp/TF2b, ΔEele is over-compensated by the desolvation energy (ΔGpol), indicating that the sum of these two components, ΔGpol + elec, is unfavorable to the binding and varies between 0.84 kcal/mol and 11.57 kcal/mol (see Table 4) and similar results are found for RdRp/myricetin, RdRp/quercetagetin, RdRp/hesperidin, and RdRp/TF1 (see Table S2 in Supplementary Information). In contrast, in the case of RdRp/TF2a, ΔGpol + elec, is favorable to the complexation (ΔGpol + elec = −0.34 kcal/mol). Overall, this suggests that the complex formation is mainly driven by the van der Waals interactions between polyphenols as well as remdesivir and RdRp. Therefore, hydrophobic residues in the binding pocket played a crucial role in the complexation process. The estimated binding free energy (ΔGbind) of remdesivir, EGCG, TF3, TF2b, and TF2a were −24.57, −27.02, −42.27, −26.74 and −27.17 kcal/mol, respectively (Table 4) and myricetin, quercetagetin, hesperidin and TF1 show lower binding affinity compared to that of remdesivir (Table S2 in the Supplementary Information). This suggests that polyphenol TF3 binds most strongly to RdRp, followed by TF2a and EGCG. The potency of the five inhibitors decreases in the following order: TF3 > TF2a > EGCG > TF2b > remdesivir. TF3 binds most strongly to RdRp because both ΔEvdW and ΔEelec are more favorable to the binding compared to the other inhibitors. Similarly, TF2a binds more strongly to RdRp compared to EGCG or TF2b because ΔGpol + elec is favorable for TF2a (ΔGpol + elec = −0.34 kcal/mol) while it is found to be unfavorable for EGCG (ΔGpol + elec = 2.24 kcal/mol) and TF2b (ΔGpol + elec = 7.83 kcal/mol). 3.2.3. Essential residues for polyphenols binding Further, to gain a deeper insight into the best four RdRp/polyphenols and remdesivir interaction pattern, the total binding free energy was decomposed into polyphenols-residue pair based on the MM-GBSA scheme. The approach of per-residue based contributions is useful to determine the binding mechanisms of an inhibitor at an atomistic level, and it also reveals the individual residue contributions. The different energy contributions from the backbone and side-chain of each residue are shown in Figure 5 and listed in Table 5. Figure 5. Decomposition of the binding free energy into contributions from individual residues for RdRp complexed with remdesivir, EGCG, TF3, TF2b and TF2a. Table 5. Per-residue based decomposition of binding free energy for the complex of remdesivir, EGCG, TF3, TF2a and TF2b with the SARS-CoV-2 RdRp. Residue TvdW Telec Tpol Tnp Tside_chain Tbackbone Ttotal RdRp/Remdesivir Asp761 1.05 −21.73 18.81 −0.14 −1.82 −0.19 −2.01 Lys798 −2.74 −4.28 5.74 −0.47 −1.76 0.01 −1.75 Pro620 −1.42 −1.15 1.16 −0.18 −1.26 −0.33 −1.59 Asp760 −0.57 −4.27 4.24 −0.12 −0.20 −0.52 −0.72 Arg553 −2.07 −3.71 5.61 −0.51 −0.76 0.08 −0.68 RdRp/EGCG Asp452 1.65 −16.48 9.85 −0.08 −5.15 0.09 −5.06 Arg553 −3.83 −4.84 6.59 −0.59 −2.68 0.01 −2.67 Pro620 −1.64 −0.41 0.51 −0.33 −1.62 −0.25 −1.87 Asp618 0.49 −8.14 6.66 −0.09 −1.13 −0.05 −1.08 Lys621 −2.26 −5.12 6.92 −0.50 −0.87 −0.09 −0.96 RdRp/TF3 Asp761 2.52 −30.13 22.11 −0.25 −5.54 −0.21 −5.75 Arg836 −0.83 −14.95 12.29 −0.37 −3.85 −0.01 −3.86 Arg555 −5.80 −4.00 7.29 −0.70 −2.98 −0.23 −3.21 Thr556 −0.20 −3.96 2.44 −0.13 −0.21 −1.64 −1.85 Ile548 −1.03 −0.08 0.12 −0.13 −0.65 −0.47 −1.12 Ser814 −1.67 0.09 0.78 −0.08 −0.48 −0.40 −0.88 Val557 −0.65 −0.16 0.17 −0.24 −0.65 −0.23 −0.88 RdRp/TF2b His816 −3.00 −0.81 1.69 −0.32 −1.63 −0.81 −2.44 Asp833 −1.29 −0.23 0.46 −0.18 −0.21 −1.03 −1.24 Tyr877 −0.92 −2.06 2.06 −0.27 −0.43 −0.76 −1.19 Glu811 0.22 −9.24 8.11 −0.13 −2.07 1.03 −1.04 His810 −1.88 −0.41 1.66 −0.34 −0.10 −0.85 −0.95 Tyr831 −1.51 −0.36 1.19 −0.13 −0.21 −0.60 −0.81 Asn815 −0.60 −0.27 0.09 −0.01 −0.25 −0.54 −0.79 RdRp/TF2a Asp618 2.63 −20.71 13.40 −0.17 −4.84 −0.01 −4.85 Arg553 −3.80 −6.52 7.55 −0.69 −2.68 −0.78 −3.46 Lys551 −2.23 −0.80 2.23 −0.45 −0.84 −0.41 −1.25 Arg555 −1.33 −0.62 1.05 −0.27 −1.12 −0.05 −1.17 Glu167 0.42 −6.82 5.92 −0.12 −0.61 0.01 −0.60 As shown in Figure 5, it was observed that residues favoring the binding of the polyphenols with RdRp include Asp452, Arg553, Arg555, Val557, Asp618, Pro620, Lys621, Asp623, Arg624, Asp760, Asp761, and Glu811, Asp833, and Arg836. Most of these residues are located in the binding site of RdRp and can form direct contacts with polyphenols and remdesivir. Figure 5 shows that amino acids Pro620, Asp761 and Lys798 for RdRp/remdesivir; Asp452, Arg553, Pro620 and Lys621 for RdRp/EGCG; Ile548, Arg555, Thr556, Asp761 and Arg836 for RdRp/TF3; Glu811, His816, Asp833 and Tyr877 for RdRp/TF2b; Lys551, Arg553, Arg555 and Asp618 for RdRp/TF2a contributed more favorably towards the binding by contributing more than −1.0 kcal/mol in size. To complement the energetic analysis, we performed MD trajectory-based hydrogen bond (h-bond) analysis for all five complexes, and the h-bonds with occupancy are listed in Table 6. The h-bonds were determined by setting the acceptor-donor distance of ≤ 3.5 Å, and the angle cut off ≥ 1200. Important h-bonds between RdRp-inhibitors are shown in Figure 6. In the case of RdRp/remdesivir, key residues involved in the hydrogen bonding are Asp761, Asp760, and Ser759, respectively. Asp760 is found to form two h-bonds with remdesivir (Asp760@OD2 - Lig@O7, Asp760@OD2 - Lig@O6) with an occupancy of more than 15% (see Table 6 and Figure 6). In the case of RdRp/EGCG, both Asp618 and Asp760 form two h-bonds with the ligand with an occupancy in the range of 16.09 to 30.17%. On the other hand, Asp761 form an h-bond with TF3 (Asp761@OD1 - Lig@O11) with an occupancy of 69.84%, while Arg836 forms two h-bonds with the ligand (Arg836@NH2 - Lig@O14, Arg836@NE - Lig@O14) with an occupancy of 52.66%, and 48.70%, respectively. Glu811, Thr556 and Asp761 also formed h-bonds with the ligand during our simulations with an occupancy varying in the range of 44% to 58% (see Table 6). In the case of RdRp/TF2b, Glu811 is found to form two strong h–bonds with the ligand (Glu811@OE1 – Lig@O7, Glu811@OE2 – Lig@O7) with an occupancy of 22.45% and 18.89%, respectively. On the other hand, it can be observed from Table 6 that Pro832 and Tyr877 form strong h-bonds (Pro832@O -Lig@O8 and Lig@O10 -Tyr877@OH) with increased occupancy (> 24%). Finally, in the case of RdRp/TF2a, Asp618 is found to form two strong h-bonds with the TF2a (Asp618@OD1 – Lig@O10, Asp618@OD1 – Lig@O11) with an occupancy of 38.68% and 38.38%, respectively. Asp760 also forms a h-bond (Asp760@O – Lig@O11) with an occupancy of 20.83%. Figure 6. Five main hydrogen bond interactions between ligands and RdRp. Table 6. Main hydrogen bond interactions formed by RdRp with remdesivir and polyphenols along with the corresponding average distance and percentage of occupancy determined using the trajectories of production simulations. Acceptor Donor Avg. Distance (Å) Occupancy (%) RdRp/Remdesivir Asp760@OD2 Lig@O7 2.66 19.46 Asp761@OD1 Lig@O7 2.65 17.70 Asp761@OD2 Lig@O6 2.63 16.86 Asp760@OD2 Lig@O6 2.65 16.65 Lig@O6 Ser759@OG 2.80 11.63 Asp760@OD1 Lig@O7 2.66 10.59 RdRp/EGCG Asp618@OD1 Lig@O5 2.61 30.13 Asp618@OD1 Lig@O6 2.61 29.28 Asp618@OD2 Lig@O5 2.61 18.25 Asp618@OD2 Lig@O6 2.61 17.38 Asp760@OD1 Lig@O5 2.63 16.09 Tyr455@OH Lig@O11 2.83 10.97 RdRp/TF3 Asp761@OD1 Lig@O11 2.61 69.84 Glu811@O Lig@O10 2.76 58.43 Thr556@O Lig@O3 2.72 56.95 Lig@O14 Arg836@NH2 2.83 52.66 Lig@O14 Arg836@NE 2.86 48.70 Asp761@OD2 Lig@O20 2.62 44.51 RdRp/TF2b Pro832@O Lig@O8 2.77 26.31 Lig@O11 Tyr877@OH 2.75 24.39 Glu811@OE1 Lig@O7 2.65 22.45 Glu811@OE2 Lig@O7 2.65 18.89 Asp833@OD2 Lig@O8 2.65 12.01 Asn874@OD1 Lig@O11 2.68 8.58 RdRp/TF2a Asp618@OD1 Lig@O10 2.59 38.68 Asp618@OD1 Lig@O11 2.62 38.38 Asp760@O Lig@O11 2.70 20.83 Asp618@OD2 Lig@O11 2.62 16.59 Asp618@OD2 Lig@O10 2.58 16.49 Asp618@OD1 Lig@O15 2.67 16.37 Finally, we supplemented the above results by analyzing the final conformation of each production simulation with the help of 2D LigPlot+ software, and different h-bonds and hydrophobic interactions were shown in Figure 7. Hydrogen bonds are depicted in green dotted lines, while red semicircle residues are involved in hydrophobic interactions. For the RdRp/remdesivir complex, Figure 7(A) displayed nine hydrophobic interactions with Lys545, Ala547, Ser549, Arg553, Val557, Asp684, Ser759, Ser814, and Arg836. This large number of interactions account for the high stability and good binding affinity of remdesivir to RdRp. EGCG formed hydrophobic interactions with Lys551, Ala554, Arg553, Arg624, Pro620, (Figure 7(B)). In the case of TF3, eight hydrophobic interactions with His439, Ile548, Ser814, Phe812, Val557, Ser549, Tyr619 and Arg555 were formed as revealed by Figure 7(C). Figure 7(D) shows that seven hydrophobic interactions with Asp833, His816, Pro832, Gln815, His872, His810 and Ser434 were formed for RdRp/TF2b. Finally, Figure 7(E) shows that RdRp/TF2a formed hydrophobic interactions with Arg555, Ala554 and Lys551. Overall, TF3 has a higher binding affinity toward RdRp compared to the other polyphenols due to a larger number of stable hydrogen bonds and hydrophobic interactions. Figure 7. The RdRp-ligands interaction profile for (A) RdRp/remdesivir, (B) RdRp/EGCG, (C) RdRp/TF3, (D) RdRp/TF2b and (E) RdRp/TF2a. The polyphenols and remdesivir are shown in balls and sticks. Hydrogen bonds are depicted in green dotted lines, and red semicircles residues are involved in hydrophobic interactions.