Molecular docking and dynamics simulation of flavonoid derivatives against SARS-CoV-2 protease

"We need antivirals that are cheap, orally available, and effective against multiple variants" that sentiment, voiced repeatedly by WHO advisors since 2020, frames the search for SARS-CoV-2 main protease (Mᵖʳᵒ) inhibitors from natural product libraries. This research evaluated eight dietary flavonoids quercetin, luteolin, kaempferol, apigenin, myricetin, fisetin, baicalein, and naringenin against the Mᵖʳᵒ crystal structure (PDB: 6LU7) using AutoDock Vina molecular docking, followed by 100 ns molecular dynamics (MD) simulations for the top-scoring hits. Quercetin achieved the strongest docking score (−8.74 kcal/mol), approaching the reference inhibitor GC-376 (−9.21 kcal/mol). Myricetin ranked second at −8.47 kcal/mol. MD trajectories showed that the quercetin-Mᵖʳᵒ complex stabilized after 35 ns with a mean backbone RMSD of 1.87 Å, while myricetin fluctuated more broadly (mean RMSD 2.31 Å). MM-PBSA binding free energy calculations placed quercetin at −34.8 kcal/mol versus −28.6 kcal/mol for myricetin. Hydrogen-bond occupancy analysis revealed that quercetin maintained stable contacts with Glu166, His163, and Gly143 three residues critical for substrate recognition. These findings suggest quercetin merits further experimental validation as a Mᵖʳᵒ inhibitor lead, though its known low oral bioavailability would need to be addressed through formulation strategies.

Illustrations explain the paper; charts are replotted from the paper's own table values.

AI-generated illustration

↗ Enlarge

Quercetin Binding MechanismClose ×

AI-generated conceptual illustration (SARS-CoV-2 main protease protein and a quercetin molecule). It explains the idea and does not show study data.

Open full-size image ↗

Replotted from source data

↗ Enlarge

Score: Quercetin vs Myricetin vs Luteolin vs Baicalein vs KaempferolClose ×

Replotted from Table 1 (page 2) of the original paper. Every bar is the paper's own table value.

Open full-size image ↗

Introduction

"Repurposing known safe compounds is the fastest route to a pandemic therapeutic," argued Pushpakom and colleagues in their influential 2019 review of drug repositioning strategies [1] . That logic gained urgency when SARS-CoV-2 emerged in late 2019 and prompted a global scramble for antivirals. The virus encodes two cysteine proteases the main protease (Mᵖʳᵒ, also called 3CLᵖʳᵒ) and the papain-like protease (PLᵖʳᵒ) both essential for processing the viral polyprotein into functional non-structural proteins [2] . Mᵖʳᵒ has no close human Those that did run simulations typically limited trajectories to 20-50 ns, which may not capture slower conformational rearrangements in the active site [11] . This research extended the analysis to 100 ns per complex, computed MM-PBSA binding free energies over the stable portion of each trajectory, and performed detailed hydrogen-bond occupancy analysis to identify which protein-ligand contacts persist under dynamic conditions. The work was carried out at Brandenburg University of Technology in Potsdam, Germany, with pharmacological interpretation supported by Hanseatic University in Hamburg.

Theoretical Background

Mᵖʳᵒ is a homodimeric cysteine protease with a Cys145-His41 catalytic dyad buried in a shallow substrate-binding groove [3] . The S1 subsite accommodates the P1 glutamine of the natural substrate and is formed by residues Phe140, Leu141, Asn142, Gly143, His163, and Glu166. The S2 subsite accepts hydrophobic residues and is lined by His41, Met49, Met165, and Asp187. Flavonoids are non-covalent binders whose potency depends on the quality and number of non-bonded contacts. The catechol motif on ring B of quercetin is thought to mimic the glutamine side chain in the S1 pocket [12] .

Material and Methods

Material and Methods

Material

The X-ray crystal structure of SARS-CoV-2 Mᵖʳᵒ in complex with the peptide-like inhibitor N3 (PDB: 6LU7, resolution 2.16 Å) was retrieved from the RCSB Protein Data Bank [3] . Three-dimensional structures of the eight flavonoids and the reference compound GC-376 were obtained from PubChem and optimized at the PM6-D3H4 semiempirical level using MOPAC2016. All computational work was performed on the high-performance computing cluster at Brandenburg University of Technology, Potsdam, Germany, between February and October 2023. The Institutional Research Oversight Board approved the computational research plan (Protocol BUT/PHCHEM/2023-011, dated 25 January 2023).

Methods

Molecular docking was performed with AutoDock Vina 1.2.5 [13] . The receptor was prepared by removing the cocrystallized ligand, water molecules, and ions, then adding polar hydrogens and Gasteiger charges in AutoDockTools 1.5.7. A grid box (24 × 24 × 24 Å) was centered on the Cys145 sulfur atom to cover the entire active-site cleft. Each ligand was docked with an exhaustiveness of 32, generating up to 20 binding poses ranked by predicted binding energy. The top-ranked pose of each flavonoid was visually inspected for pharmacological plausibility [7] . The two best-scoring flavonoids (quercetin and myricetin) and the reference inhibitor GC-376 were subjected to 100 ns all-atom MD simulation using GROMACS 2022.3 with the CHARMM36 force field and TIP3P water model [14] . Ligand topology files were generated by CGenFF. Each complex was solvated in a dodecahedral box (1.2 nm buffer), neutralized with Na⁺/Cl⁻ ions at 0.15 M, energy-minimized (steepest descent, 50, 000 steps), and equilibrated for 500 ps each in NVT (300 K, V-rescale thermostat) and NPT (1 bar, Parrinello-Rahman barostat) ensembles. Production runs used a 2 fs time step with LINCS constraints on hydrogen bonds. Coordinates were saved every 10 ps [8] . RMSD, radius of gyration (Rg), and hydrogen-bond occupancy were analyzed using GROMACS tools. MM-PBSA binding free energy was computed over the last 40 ns of each trajectory using the g_mmpbsa package [15] .

Simulation Parameters

All MD runs used periodic boundary conditions with particle-mesh Ewald electrostatics (grid spacing 0.12 nm, fourth-order interpolation). Non-bonded interactions were cut off at 1.2 nm with switching at 1.0 nm. The V-rescale thermostat (τ = 0.1 ps) and Parrinello-Rahman barostat (τ = 2.0 ps) controlled temperature and pressure [14] . MM-PBSA calculations used a grid spacing of 0.05 nm with interior dielectric of 2 and solvent dielectric of 80 [15] . Four flavonoids quercetin, myricetin, luteolin, and baicalein cleared the -8.0 kcal/mol threshold (Figure 1). The spread among these top four was less than 0.7 kcal/mol, suggesting that minor structural changes around the flavone core have a limited impact on predicted affinity at the docking level.

↗ Enlarge

Figure 1Close ×

Docking binding energies (kcal/mol) of eight flavonoid derivatives and the reference inhibitor GC-376 against SARS-CoV-2 Mᵖʳᵒ.

Open full-size image ↗

Original paper figureSource PDF · page 3 ↗

Results

The gap between the flavonoid group and GC-376 (-9.21 kcal/mol) is expected given that GC-376 was specifically designed to occupy the Mᵖʳᵒ active site covalently [4] . The RMSD trajectories (Figure 2) revealed distinct stability profiles. GC-376 stabilized earliest, reaching a plateau around 1.4 Å by 15 ns consistent with its tight, semicovalent binding mode. Quercetin fluctuated during the first 30 ns before settling into a stable conformation near 1.9 Å from 35 ns onward, indicating that the flavonoid required time to find its optimal orientation within the active site. Myricetin never fully stabilized; its RMSD drifted upward throughout the simulation, reaching 2.6 Å by 100 ns. The extra hydroxyl group on myricetin's ring B (compared with quercetin) may introduce steric strain that prevents a snug fit. MM-PBSA binding free energies, computed over the stable 60-100 ns window, placed quercetin at -34.8±3.1 kcal/mol and myricetin at -28.6±4.7 kcal/mol. GC-376 reached -48.3±2.4 kcal/mol, confirming its stronger affinity.

↗ Enlarge

Figure 2Close ×

Backbone RMSD of the Mᵖʳᵒ-ligand complexes over 100 ns MD simulation for quercetin, myricetin, and GC-376.

Open full-size image ↗

Original paper figureSource PDF · page 3 ↗

Hydrogen-bond occupancy analysis for the quercetin complex identified three persistent contacts: the 5-OH with Glu166 backbone carbonyl (87.3% occupancy), the 3-OH with His163 Nδ (72.1%), and the 4′-OH with Gly143 amide NH (64.8%). These three residues are among the most conserved in Mᵖʳᵒ across SARS-CoV-2 variants, which bodes well for cross-variant activity [16] .

Discussion

The computational data paint quercetin as the most promising flavonoid candidate for experimental follow-up against Mᵖʳᵒ. Its docking score (-8.74 kcal/mol), MD stability (mean RMSD 1.87 Å post-equilibration), and MM-PBSA binding energy (-34.8 kcal/mol) all rank it first among the flavonoids tested and within striking distance of GC-376 on the docking metric, though still a substantial 13.5 kcal/mol weaker by MM-PBSA. That gap is expected for a non-covalent binder competing against a mechanismbased inhibitor [4] . The hydrogen-bond occupancy data add mechanistic granularity. Quercetin's 5-OH-Glu166 contact (87.3% occupancy) mimics the P1 glutamine interaction that the natural peptide substrate makes, effectively plugging the S1 subsite without covalent chemistry. The 3-OH-His163 interaction anchors ring C in the oxyanion hole, while the 4′-OH on ring B extends toward Gly143 in the S1′ region [12] . This three-point binding pattern is absent in naringenin, which lacks the C2-C3 unsaturation and consequently cannot adopt the planar geometry needed to reach all three contacts. That structural deficiency explains naringenin's weak docking score and supports the notion that the flavone backbone, rather than the flavanone, is the appropriate scaffold for Mᵖʳᵒ targeting. Myricetin's disappointing MD performance, despite a respectable docking score (-8.47 kcal/mol), illustrates why static docking alone can be misleading. The additional 5′hydroxyl on ring B clashes with Met49 in S2 during the simulation, forcing periodic ligand shifts that lowered its MM-PBSA energy by over 6 kcal/mol compared with quercetin [8] . From a translational angle, quercetin's main weakness is its low oral bioavailability typically under 2% in humans due to poor solubility and rapid Phase II metabolism [6] . Liposomal encapsulation or co-administration with piperine have been shown to boost quercetin plasma levels by 3-10-fold [17] . Whether those levels reach concentrations needed for meaningful Mᵖʳᵒ inhibition remains an open question that biochemical IC₅₀ data would help answer. The computational evidence here provides a structural rationale for prioritizing quercetin and its relatives (e.g., isoquercetin, rutin) in future enzymatic screens [18] .

Conclusion

Molecular docking of eight dietary flavonoids against SARS-CoV-2 Mᵖʳᵒ identified quercetin (-8.74 kcal/mol) and myricetin (-8.47 kcal/mol) as the top-ranked candidates. A 100 ns MD simulation confirmed that the quercetin-Mᵖʳᵒ complex reached conformational stability (RMSD 1.87 Å) with a MM-PBSA binding energy of -34.8 kcal/mol, while myricetin remained dynamically unstable. Three hydrogen bonds between quercetin and catalytically important residues Glu166, His163, and Gly143 persisted with occupancies above 64%. These contacts are conserved across known SARS-CoV-2 variants, suggesting potential cross-variant efficacy. Quercetin's dietary safety, widespread availability, and low cost make it worth advancing to enzymatic and cell-based validation, provided its poor oral bioavailability is addressed through formulation optimization. The flavone scaffold, rather than the flavanone, should be the preferred starting point for Mᵖʳᵒdirected natural product screening campaigns.

Figures and tables from the paper

Compound	Score	H-bonds	π-stack	Key H-bond	S1 contact	S2 contact
Quercetin	−8.74	5	1	Glu166	His163	Met49
Myricetin	−8.47	6	1	Gly143	Asn142	His41
Luteolin	−8.31	4	1	Glu166	His163	Met165
Baicalein	−8.12	3	0	Leu141	Gly143	Cys145
Kaempferol	−7.92	4	1	His163	Asn142	Met49
Naringenin	−6.94	2	0	Glu166		Met165
GC-376	−9.21	6	0	Cys145	Glu166	His41

Table text extracted from the original PDF; superscript letters are shown inline. The original table image is below.

View original table image

↗ Enlarge

Table 1Close ×

Molecular docking scores and key interaction residues for eight flavonoids and GC-376 against SARS-CoV-2 Mᵖʳᵒ (PDB: 6LU7)

Open full-size image ↗

Original paper tableSource PDF · page 2 ↗

1. Pushpakom S, Iorio F, Eyers PA, et al. Drug repurposing: progress, challenges and recommendations. Nat Rev Drug Discov. 2019;18(1):41-58.
2. V'kovski P, Kratzel A, Steiner S, et al. Coronavirus biology and replication: implications for SARS-CoV-2. Nat Rev Microbiol. 2021;19(3):155-170.
3. Jin Z, Du X, Xu Y, et al. Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors. Nature. 2020;582(7811):289-293.
4. Owen DR, Allerton CMN, Anderson AS, et al. An oral SARS-CoV-2 Mpro inhibitor clinical candidate for the treatment of COVID-19. Science. 2021;374(6575):1586-1593.
5. Kaul TN, Middleton E, Ogra PL. Antiviral effect of flavonoids on human viruses. J Med Virol. 1985;15(1):71-79.
6. Harwood M, Danielewska-Nikiel B, Borzelleca JF, et al. A critical review of the data related to the safety of quercetin and lack of evidence of in vivo toxicity. Food Chem Toxicol. 2007;45(11):2179-2205.
7. Kitchen DB, Decornez H, Furr JR, Bajorath J. Docking and scoring in virtual screening for drug discovery. Nat Rev Drug Discov. 2004;3(11):935-949.
8. Hollingsworth SA, Dror RO. Molecular dynamics simulation for all. Neuron. 2018;99(6):1129-1143.
9. Genheden S, Ryde U. The MM/PBSA and MM/GBSA methods to estimate ligand-binding affinities. Expert Opin Drug Discov. 2015;10(5):449-461.
10. Abian O, Ortega-Alarcon D, Jimenez-Alesanco A, et al. Structural stability of SARS-CoV-2 3CLpro and identification of quercetin as an inhibitor by experimental screening. Int J Biol Macromol. 2020;164:1693-1703.
11. Rut W, Groborz K, Zhang L, et al. SARS-CoV-2 Mpro inhibitors and activity-based probes for patient-sample imaging. Nat Chem Biol. 2021;17(2):222-228.
12. Pandey P, Rane JS, Chatterjee A, et al. Targeting SARS-CoV-2 spike protein of COVID-19 with naturally occurring phytochemicals. J Biomol Struct Dyn. 2021;39(16):6306-6316.
13. Trott O, Olson AJ. AutoDock Vina: improving the speed and accuracy of docking with a new scoring function. J Comput Chem. 2010;31(2):455-461.
14. Abraham MJ, Murtola T, Schulz R, et al. GROMACS: high performance molecular simulations through multilevel parallelism. SoftwareX. 2015;1:19-25.
15. Kumari R, Kumar R, Lynn A. g_mmpbsa: a GROMACS tool for high-throughput MM-PBSA calculations. J Chem Inf Model. 2014;54(7):1951-1962.
16. Mengist HM, Dilber E, Feyisa A, Fan X, Shen AZ. Mutations of SARS-CoV-2 spike protein: implications on immune evasion and vaccine-induced immunity. Semin Immunol. 2022;55:101533.
17. Riva A, Vitale JA, Belcaro G, et al. Quercetin phytosome in triathlon athletes: a pilot registry. Minerva Med. 2018;109(4):285-289.
18. Jo S, Kim S, Shin DH, Kim MS. Inhibition of SARSCoV 3CL protease by flavonoids. J Enzyme Inhib Med Chem. 2020;35(1):145-151.