Comprehensive virtual screening of 4.8 k flavonoids reveals novel insights into allosteric inhibition of SARS-CoV-2 MPRO

Abstract

SARS-CoV-2 main protease is a common target for inhibition assays due to its high conservation among coronaviruses. Since flavonoids show antiviral activity, several in silico works have proposed them as potential SARS-CoV-2 main protease inhibitors. Nonetheless, there is reason to doubt certain results given the lack of consideration for flavonoid promiscuity or main protease plasticity, usage of short library sizes, absence of control molecules and/or the limitation of the methodology to a single target site. Here, we report a virtual screening study where dorsilurin E, euchrenone a11, sanggenol O and CHEMBL2171598 are proposed to inhibit main protease through different pathways. Remarkably, novel structural mechanisms were observed after sanggenol O and CHEMBL2171598 bound to experimentally proven allosteric sites. The former drastically affected the active site, while the latter triggered a hinge movement which has been previously reported for an inactive SARS-CoV main protease mutant. The use of a curated database of 4.8 k flavonoids, combining two well-known docking software (AutoDock Vina and AutoDock4.2), molecular dynamics and MMPBSA, guaranteed an adequate analysis and robust interpretation. These criteria can be considered for future screening campaigns against SARS-CoV-2 main protease.

Figure 1

SARS-CoV-2 M^PRO and its three Putative Binding Sites (PBS). (a) Protomers are in cartoon and surface representations, each in a different colour; black boxes indicate the position of PBSs: (A) Substrate Binding Site, (B) Dimer Site, (C) Cryptic Site. (b) Residues conforming the binding sites are depicted in the green protomer as sticks according to colour code: blue for dimerization and magenta for cryptic binding sites, respectively. Substrate binding site residues and labels are coloured according to 5 different subregions: red (S1’: His41, Gly143, Ser144 and Cys145^,), blue (S1: Ser1 (from the other protomer, not shown), Phe140, Leu141, Asn142, His163, Glu166 and His172), purple (S2: His41, Met49, Tyr54, Met165 and Asp187), orange (S3: Met165, Glu166 and Gln189), grey (S4: Leu167, Pro168, Phe185, Thr190, Ala191, Gln192 and Gln189^,) and cyan (S5: Pro168, Ala191). Residues with two colours belong to two subregions; in this case, text labels are coloured black. (c) Close-up view of SBS. Labels from catalytic residues are enclosed in a red rectangle. (d) Close-up view of DS. DS residues were defined according to the literature: Arg4, Met6, Ser10, Gly11, Glu14, Asn28, Ser139, Phe140, Ser147, Glu166, Glu290 and Arg298. (e) Close-up view of CS. CS residues were predicted using CryptoSite: Lys5, Met6, Pro108, Gly109, Arg131, Trp218, Phe219, Tyr239, Glu240, Leu271, Leu272, Leu287, Glu288, Asp289, Glu290, Arg298, Gln299 and Val303.

Figure 2

Structural alignment between positive controls and the original M^PRO crystal. Docked non-covalent inhibitors X77 (PDB ID: 6W63), X7V (PDB ID: 7KX5) and YD1 (PDB ID: 7LTJ) were aligned to their original M^PRO crystals. Docking was carried out following the exhaustive-ranking procedure using AutoDock-GPU. RMSD was calculated without considering hydrogen atoms. (a) Representative ligand X77 in the substrate binding site (black circle). All controls were also bound to this region. The M^PRO crystal is shown as in Fig. 1. Ligands are shown in sticks: orange for the control and green for the cross-docked inhibitors. Non-carbon atoms are coloured following the CPK colouring convention. Close-ups are provided to observe the docked pose in detail compared to the original position: (b) X77, RMSD = 1.352 Å, MBE = − 10.77 kcal/mol; (c) X7V, RMSD = 2.027 Å, MBE = -9.93 kcal/mol; (d) YD1, RMSD = 1.048 Å, MBE = − 8.00 kcal/mol. MBE, Mean Binding Energy. X77 was selected as the positive control for further analysis because it had the best docking energy.

Figure 3

RMSD plots of flavonoid-M^PRO selected ligands along molecular dynamics simulations. Positive and negative controls are also shown. As negative control only remained associated during the first 10 ns, MD was not extended. RMSD of the protein, complex and flavonoid ligand is shown in red, black and blue, respectively. Green vertical lines reflect snapshots taken for further analysis with MMPBSA methods. The grey box indicates residence time (RT, time in which the ligand is interacting with the protein). In plots without a grey box, the ligand never detaches from the protein. SBS: (a) Dorsilurin E, (b) euchrenone a11 (detaches at 500 ns), (c) X77 (positive control, detaches at 200 ns), (d) TDZD-8 (negative control, detaches at 10 ns). In the SBS, the second lowest RMSD and fluctuations within RT pertained to dorsilurin E (< 5 Å during the first 100 ns, < 16 Å during all the RT, with exception of some spikes), which suggest stable binding. DS: (e) sanggenol O. CS: (f) CHEMBL2171598.

Figure 4

Sanggenol O-induced conformational changes in region A. (a) Snapshots of sanggenol O throughout the MD are shown each 200 ns. Residues from region A are shown in green surface while the rest of the protein is shown in cyan. Ligand backbone is coloured blue and other atoms follow the CPK colouring convention. (b) Ligand-induced RMSD variations in SBS residues are shown in red. RMSD of the apo-monomer structure is provided for comparison (black).

Figure 5

CHEMBL2171598 bound to region B and promoted the reorientation of M^PRO domains I/II. (a) Superimposition of the initial M^PRO complex (cyan; frame 0) to the 600 ns snapshot (green; frame 59,999) over the domain III. CHEMBL2171598 is shown as sticks. Movements of key residues are shown with black arrows. (b–d) MD snapshots detailing domain reorganization and the hinge mechanism. Key residues Arg4 (yellow), Arg298 (green) and Gln299 (gray) are displayed as sticks along with their interactions (red dotted lines for π-cation, black for electrostatic). Residues enclosed in [brackets]* are engaging in a polar contact. C-terminal helices (residues 244–257 and 293–306) that interact with the ligand are coloured orange. Lower panels are rotated 90° with respect to upper panels to display interactions between domains I/II and III. Residues participating in those interactions are shown as cyan sticks. Gln299 was removed for clarity. (b) CHEMBL2171598 induced the rotation of Arg4. Two polar contacts are observed between domains I/II and domain III: Gly109-Asn203 and Thr111-Asp295. (c) Arg4’s rotation reoriented Arg4 and Arg298 side chains in syn conformation. In their new syn conformation, Arg298 and Arg4 guanidinium groups interrupt the interdomain polar contacts. (d) Arg4 displaced Arg298 while CHEMBL2171598 adopted its final position, named Region B. Arg4’s insertion to the cleft made impossible for interdomain polar contacts to form again. The hinge movement was triggered shortly after, at 465 ns.

Figure 6

RMSF and hydrogen bond plots of flavonoid-SARS-CoV-2 M^PRO selected ligands. In both plot types, ligands are shown according to previous colour code: green for substrate, blue for dimeric and magenta for cryptic sites. (a–d) Protein RMSF plots. RMSF of apo- dimeric or monomeric (protomer A) protein is shown in black. Red arrows show the most significant differences between apo- and holo- structures. (e–g) Total number of hydrogen bonds formed between SARS-CoV-2 M^PRO and flavonoid ligands. (a/e) Dorsilurin E restrained the flexibility of residue 49. This ligand also had a long-range effect by reducing the flexibility of domain III residues 200–275. Dorsilurin E held the lowest average number of hydrogen bonds, 0.8 (s = 0.99). (b/f) Euchrenone a11 restrained the flexibility of residue 49. Additionally, euchrenone a11 increased fluctuation in domain II (S1’ subpocket) residues 143–146. This ligand also held the second lowest average number of hydrogen bonds, 0.9 (s = 0.94). (c/g) Sanggenol O restrained the protein movement except for residues 186–198, which showed flexibilization. Sanggenol O shows remarkable fluctuations in its hydrogen bonds, holding an average of 1.6 hydrogen bonds through time (s = 1.30). (d/h) CHEMBL2171598 increased the protein flexibility throughout its sequence, with a marked increase in domain II (residues 100–182) and residues 210–307 of domain III. CHEMBL2171598 held heterogenous hydrogen bonds through MD with an average of 0.9 (s = 1.01).

Figure 7

Mean and standard deviation values of binding energies calculated by the MMPBSA method. Mean binding energies are calculated from a N = 150 sample. Whiskers indicate standard deviations. Each N = 150 set was sampled by using a particular snapshot of the MD. Three snapshots were retrieved for each ligand. Darker blue filling indicates that snapshots have been taken from further along the holo MD. PBS: (a) SBS. Dorsilurin E had the most favourable energies, that started from − 22.66 kcal/mol (20 ns, s = 0.85) and gradually increased to − 4.34 (600 ns, s = 0.77). Its mean overall binding energy was − 8.25 kcal/mol. Positive control X77 had the second-best mean binding energy over time, with a peak at − 12.24 kcal/mol (150 ns, s = 1.37). Euchrenone a11 binding energies started from − 12.09 kcal/mol (250 ns, s = 0.88), decreased to − 14 kcal/mol (350 ns, s = 0.93) and then, closer to the unbinding event, stayed at − 4.27 kcal/mol (450 ns, s = 0.95). (b) DS. CHEMBL2171573 showed the lowest mean binding energy at − 18.33 kcal/mol, which decreased over time with a peak at − 23.63 kcal/mol (600 ns, s = 0.89). CHEMBL2171584 also showed decreasing mean binding energies over time that peaked at − 18.30 kcal/mol (550 ns, s = 1.04). Kanzonol E had an initial favourable value (− 20.26 kcal/mol, 175 ns, s = 0.64) that increased with time (− 12.27 kcal/mol, 600 ns, s = 0.97). This ligand also showed the least standard deviation values. (c) CS. CHEMBL2171598 started with a mean binding energy (− 21.72 kcal/mol, 400 ns, s = 0.98) that increased at 550 ns (− 17.23 kcal/mol, s = 1.22), only to decrease again at 600 ns (− 18.57 kcal/mol, s = 1.25). Overall, it showed the lowest mean binding energy at − 19.18 kcal/mol.

Figure 8

Interaction types for selected ligands as percentage of total interactions. Interaction types are limited to hydrogen bonds (blue), hydrophobic interactions (red), π-Cation contacts (yellow), and π-π stacking contacts (green). Residues with under 2.5% of interactions are not shown. Asterisks (*) indicate residues as belonging to protomer B, the default being protomer A. Note that not all ligands share the same number of interactions, so equal percentages do not mean equal number of interactions. SBS: (a) dorsilurin E, (b) euchrenone a11. DS: (c) sanggenol O. CS: (d) CHEMBL2171598.

Figure 9

SASA values from SBS in each protomer along simulation time. SASA values correspond to protomer A and protomer B SBSs for (a) apo dimeric M^PRO or holo dimeric M^PRO with (b) euchrenone a11, (c) dorsilurin E or (d) X77 (positive control). Circles represent the real measure of SASA during simulation time for protomer A, while the same is true for triangles and protomer B. Lines are the smoothed values for each protomer (A, red; B, blue) that help in tendency recognition. Smoothing was done through locally estimated scatterplot smoothing (LOESS) using a span of 0.005. As was expected, X77 occupancy substantially increased overall SASA in protomer A SBS, thereby heightening the alternation pattern even without pronounced SBS closing in protomer B. Therefore, it is hypothesized that X77 would be exerting the same effect seen in the natural substrate and thus its main inhibition mechanism would not be the disruption of the alternation.