3.1 Convergence of MD Simulations
Evaluation of RMSD and RG revealed that systems reached equilibrium in 10–20 ns with average RMSG and RG values that fluctuated between 1.5 ± 0.2 and 2.4 ± 0.2 Å, and 25.7 ± 0.3 and 26.1 ± 0.2 Å, respectively (Table 1, supplementary material). Values determined here were similar to those previously observed for the bound dimeric system [7]. Therefore, the first 30 ns of the 100 ns simulation were excluded from further analyses.
3.2 MD Simulations of Ligand-Dimeric SARS-CoV2 Mpro Binding
MD simulations showed that ligands remained bound to both subunits of the dimer, with the exception of systems assessing rutin and ledipasvir with dimeric SARS-CoV2 Mpro. In this case rutin diffused from both subunits of SARS-CoV2 Mpro, and ledispavir remained bound to one dimeric SARS-CoV2 Mpro catalytic site. Peptide-like inhibitor N3 was bound to subunit 1 via hydrophobic contacts involving T25, L27, H41, C44, T45, S46, M49 and Q189. N3 also established polar interactions with the side chain of Q189 (Fig. 1A). At subunit 2 inhibitor N3 was stabilized by contact with L27, H41, M49, N142, G143, S144, C145, M165, P168, D187 and Q189, and hydrogen bonds with backbone C145, and side chain S144 (Fig. 1B). Of these residues, T25, H41, M49, N142, G143, S144, C145, M165, P168, D187 and Q189 were observed in the co-crystallized complex (PDB code 6LU7).
Darunavir in subunit 1 was bound by L27, H41, M49, F140, N142, G143, S144, C145, H163, M165, D187 and Q189 via non-polar interactions, and formed hydrogen bonds with backbone G143 and D187 residues (Fig. 1C). On subunit 2, coordination of darunavir occurred through non-polar interactions with T25, H41, C44, S46, M49, N142, C145, M165, D187, R88 and Q189, and polar interactions with T25 and N142 side chains and backbone atoms of H41 and C44 (Fig. 1D).
Indinavir was stabilized at subunit 1 through non-polar interactions with H41, M49, M165, E166, P168, R188, Q189, A191 and Q192, and established polar interactions with the side chain of Q189 (Fig. 1E). On subunit 2, the ligand was coordinated by H41, M49, C145, M165, E166, L167, P168, A191 and Q192 through hydrophobic contacts, and by polar interactions with backbone and side chain atoms of E166 (Fig. 1F).
Saquinavir bound subunit 1 through non-polar interactions with T25, H41, S46, M49, F140, L141, N142, S144, C145, M165, E166, P168, D187 and polar interactions with the side chain of Q189 and backbone atoms of D187 (Fig. 2A). Saquinavir bound subunit 2 via non-polar interactions involving H41, V42, M49, F140, L141, N142, S144, C145, M165, E166, P168, D187 and Q189, and through polar interactions with side chain groups of N142 and Q189, and with backbone atoms of E166 (Fig. 2B).
Tripanivir bound subunit 1 through hydrophobic contacts with T25, T26, L27, H41, V42, C44, M49, N142, G143, S144, C145 and M165, and polar interactions with backbone atoms of G143 and T26 (Fig. 2C). On subunit 2, tripanivir was bound via polar interactions with T24, T25, L27, C44, T45, S46, M49, N142 and G143, and polar interactions with the side chain of S46 (Fig. 2D).
Interactions between diosmin and subunit 1 were stabilized through hydrophobic interactions with T25, L27, H41, M49, P52, F140, N142, G143, S144, C145, H163, M165, E166, R188 and Q189, polar interactions with backbone atoms of M49 and S144, and side chain atoms of H41 and E166 (Fig. 2E). On the subunit 2 diosmin bound H41, M49, M165, L167, P168, D187, R188, Q189 and Q192 via non-polar interactions, and backbone atoms of R188 through polar interactions (Fig. 2F).
Hesperidin bound subunit 1 through hydrophobic contacts with L27, H41, T45, C145, H164, M165, E166, L167, P168, D187, Q189, T190 and A191, and hydrogen bonds with the side chains of H41 and E166 (Fig. 3A). On subunit 2, hesperidin bound via hydrophobic interactions with L27, H41, M49, N142, C145, M165 and Q189, and formed polar interaction with the side chain of the latter residue (Fig. 3B).
Raltegravir bound subunit 1 through non-polar interactions with T25, L27, H41, T45, S46, M49 and N142 (Fig. 3C). On subunit 2, binding to raltegavir occurred through polar contacts with H41, M49, M165, D187, R188 and Q189, and hydrogen bonding with R188 (Fig. 3D). Velpatasvir was stabilized at subunit 1 through hydrophobic contacts with T25, L27, S46, M49, L50, G143, C145, Q189, T190 and A191, and by one hydrogen bond with the side chain of S46 (Fig. 3E). In contrast, at subunit 2 the ligand bound via non-polar interactions with T25, T26, L27, S46, M49, L50, N142, G143, P168, Q189, T190, A191, and polar interaction with backbone atoms of T26, Q189 and A191 (Fig. 3F).
Ledipasvir bound subunit 1 via hydrophobic contacts with L141, N142, Q189, T190, A191, and polar interactions with both the backbone and side chain of N142 (Fig. 4A).
Rosuvastatin was stabilized through hydrophobic interactions with T25, L27, H41, M49, N142, C145, H163, M165, E166, Q189 and Q192, and formed polar interaction with the side chain of the latter residue (Fig. 4B). At subunit 2 rosuvastatin was stabilized via hydrophobic interactions with T25, L27, H41, T45, N142, G143, C145, M165 and Q189 (Fig. 4C).
Bortezomib bound subunit 1 through hydrophobic contacts with L27, H41, M49, L141, N142, G143, S144, F140, C145, H163 and M165, and established polar interactions with the side chain of N142, and backbone atoms of G143 and L27 (Fig. 4D). At subunit 2, bortezomib bound L27, H41, C44, T45, A46, M49, M165, D187, R188, and Q189 via hydrophobic interactions, and formed hydrogen bonds with backbone atoms of R188, and the side chain of Q189 (Fig. 4E).
Binding was primarily stabilized via non-polar interactions involving L27, H41, M49 and M165 residues. T25, T26, H41, C44, S46, M49, N142, G143, S144, C145, E166, D187, R188, Q189, A191 and Q192 established polar interactions through their backbones or side chains with some of the compounds assessed (Fig. 2–4). H41 and C145 are part of the catalytic site, and H41, M49, G143, S144, M165, E166, D187, R188, Q189, A191 and Q192 form part of the substrate binding region [28,29], indicating that all ligands fit into their respective active sites. In addition, we observed that complexes involving each subunit were generally stabilized by an uneven number of residues. This suggested that the conformations of the two catalytic sites of the dimer were variable, which was a result not observed via crystallographic methods.
3.3 Free Energy of Binding
Affinity for complex formation was calculated using the MM/GBSA approach. We observed that all binding was thermodynamically favorable and occurred through the formation of non-polar interactions, van der Waals energy (ΔEvdw) and non-polar free energy of desolvation (ΔGnpol,sol). ΔGbind values for compounds bound to the first subunit of dimeric SARS-CoV2 Mpro demonstrated the following binding tendencies: saquinavir > tipranavir > darunavir > diosmin > rosuvastatin > indinavir > bortizomib > peptide-like inhibitor N3 > velpatasvir > hesperidin > raltegravir. ΔGbind values associated with binding to the second subunit were as follows: saquinavir > indinavir > hesperidin > darunavir > velpatasvir = ledipasvir > peptide-like inhibitor N3 > ratelgravir > diosmin > tipranavir = rosuvastin > bortezomib (Table 1). Differing affinity for each subunit is consistent with observed differences in the number of interactions involved in binding to each subunit, which were determined through structural analyses (Fig. 1-4). This highlights the importance of evaluating affinity to both subunits using end-point free energy methods. Binding tendency determined here is in agreement with a previous study in which enhanced inhibitory properties were associated with saquinavir relative to darunavir [30], but contrasting with the enhanced inhibitory activity experimentally determined for nelfinavir compared to saquinavir. This may suggest nelfinavir inhibits activity of more than one target [31]. Saquinavir and darunavir had increased affinity for dimeric SARS-CoV Mpro than that which was previously determined for praziquantel, perampanel, ritonavir and nelfinavir using similar methodology, whereas affinity of darunavir was similar to that which had previously been reported for lopinavir [7].
Based on this analysis, it is clear that both saquinavir and darunavir, both known to possess potent antiviral protease inhibitory activity [32,33], are attractive anti-COVID-19 clinical drug candidates. The two drugs had stronger affinity for both subunits of Mpro than that which was experimentally determined for peptide-like inhibitor N3, which displayed strong antiviral effects in the μM concentration range in SARS-CoV-2 virus-infected Vero cells [34]. Assessment of ΔGbind values for indinavir and tipranavir indicated that they may also possess moderate activities against COVID-19. A comparison of the 12 FDA-approved compounds evaluated, after excluding ledispavir and rutin that exhibited ligand diffusion in one or both subunits, revealed that the worst candidates to inhibit activity of dimeric SARS-CoV Mpro were raltegravir and bortezomib. These compounds had the greatest affinities of all inhibitors assessed via virtual screening, which employed docking studies involving monomeric SARS-CoV Mpro [10,11].
Comparative analysis of the affinity tendency observed for three of the 12 evaluated compounds previously reported by Farag et al (Darunavir > rosuvastatin > saquinavir) [8], Chen et al (ledipasvir > velpatasvir) [9], Adem et al (hesperidin > rutin > diosmin) [35], and Kumar et al (tipranavir > raltegravir) [10] revealed discrepancies when compared with our findings. These observations highlight the degree to which the computational strategy employed to identify potential inhibitors of SARS-CoV Mpro impacted affinity tendency, and underscored the utility of combining docking, MD simulations and end-point binding free energy methods.
3.4 Per-residue Free Energy Decomposition
Investigation of the residues contributing to ΔGbind in ligands-dimeric SARS-CoV2 complex revealed that complex stabilization could generally be attributed to 5 to 14 residues (Table 2-4). Complex formation between dimeric SARS-CoV2 Mpro and peptide-like inhibitor N3 involved T25, C44, T45 and S46 of subunit 1. The energetic contribution of N142, G143, S144, C145, M165, P168 and D187 residues were only observed when ligands bound subunit 2 (Table 2), while L27, H41, M49, and Q189 contributed to the stabilization of binding to both subunits. For darunavir, L27, F140, G143, S144 and H163 stabilized binding subunit 1, and T25, C44, S46 and R188 stabilized binding to subunit 2. H41, M49, N142, C145, M165, D187 and Q189 residues were involved in binding to both subunits.
For indinavir, R188 and Q189 contributed stabilized binding to subunit 1, whereas tC145 and L167 facilitated subunit 2 binding. H41, M49, M165, D166, P168, A191 and Q192 were involved in binding to both subunits. For saquinavir, T25 and S46 of subunit 1 stabilized binding, and V42 stabilized interactions involving subunit 2. H41, M49, F140, L141, N142, S144, C145, M165, D166, P168, D187 and Q189 contributed energetically to binding to both subunits. T26, H41, V42, S144, C145 and M165 residues of subunit 1, and T24, T45 and S46 of subunit 2, were involved in binding to tipranovir. T25, L27, C44, M49, N142 and G143 facilitated binding to both subunits of Mpro. T25, L27, P52, F140, N142, G143, S144, C145, H163 and E166 stabilized diosmin binding to subunit 1, whereas L167 and P168 promoted binding to subunit 2. H41, M49, M165, R188 and Q189 contribute at both subunits of dimeric SARS-CoV2 Mpro (Table 3).
T45, H164, E166, L167, P168, D187, T190 and A191 facilitated hesperidin-subunit 1 binding, and M49 and N142 exclusively facilitated binding to subunit 2. L27, H41, C145, M165 and Q189 of both subunits of dimeric SARS-CoV2 were involved in binding hesperidin (Table 3). T25, L27, T45, S46 and N142 of subunit 1 were determined to be the principal stabilizers of raltegravir binding, whereas that M165, D187, R188 and Q189 of subunit 2 promoted binding to the inhibitor. H41 and M49 were involved in binding both Mpro subunits (Table 3). C145 was the only residue of subunit 1 involved in binding velpatasvir, while T26, N142 and P168 of subunit 2 stabilized binding. While T25, L27, S46, M49, L50, G143, Q189, T190 and A191 of both subunits were involved in binding (Table 3). Regarding ledipasvir, ligand only remained bound to subunit 2 via the energetic contributions of L141, N142, Q189 and T190 (Table 4). With regard to rosuvastatin, M49, H163, E166 and Q192 of subunit 1, while and T45 and G143 of subunit 2 stabilized interactions with ligands. L27, H41, N142, N145, M165 and Q189 of both subunits facilitated binding. For bortezomib is was observed that F140, L141, N142, G143, S144, C145 and H163 stabilized ligand binding to subunit 1. C44, T45, S46, D187, R188 and Q189 stabilized ligand binding to subunit 2. However, L27, H41, M49 and M165 stabilized ligand binding to both subunits.
An analysis of residues that enhanced ligand binding affinity revealed that, in general, L27, T25, T26, H41, M49, V42, T45, L50, S46, F140, N142, G143, S144, C145, H163, M165, D166, P168, D187, R188, Q189, T190, A191 and Q192 significantly contributed to ΔGbind values of some complexes, but only L27, H41, M49 and M165 were involved in protein-ligand complex formation consistently, This highlights the importance of residues that comprise the catalytic (H41) and substrate binding (M49 and M165) sites [28,29], and reveals the importance of the participation of other residues involved in ligand stabilization, which were also identified previously [7].
3.5 Principal Component Analysis (PCA)
PCA analysis allowed researchers to quantitatively approximate the degree of mobility change that occurred upon ligand complexation. Therefore, the trace of the diagonalized covariance matrix of backbone atoms was calculated for bound SARS-CoV2 Mpro systems. This analysis showed that darunavir, hesperidin and raltegravir binding was linked to a reduced degree of conformational change of dimeric SARS-CoV2 Mpro relative to that which occurred when the apo state [7] of the protein was bound. Peptide-like inhibitor N3, indinavir, saquinavir, tipranavir, velpatasvir, and bortezomib binding to dimeric SARS-CoV2 Mpro produced a small increase in conformational mobility. However, diosmin and ledipasvir binding did not produce conformational change upon ligand binding. Importantly, binding of rosuvastatin increased the conformational mobility of dimeric SARS-CoV2 Mpro, which suggested binding of most compounds was not likely to significantly impact affinity for rosuvastatin (Table 1), while decreases in the predicted affinity for rosuvastatin could be expected as a result of favorable entropic components that affect the degree of conformational mobility upon complex formation.