In Silico Design of New Dual Inhibitors of SARS-CoV-2 M PRO through Ligand-and Structure-Based Methods

SARS-CoV-2 is a highly transmissible β-coronavirus with a positive single-stranded RNA genome, related to SARS-CoV and MERS-CoV. In the context of the COVID-19 pandemic, drug repurposing provided rapid therapeutic options to complement vaccination. Among viral targets, the main protease (Mpro, 3CLpro) is pivotal for viral replication, cleaving viral polyproteins at 11 conserved sites to generate non-structural proteins. Mpro specifically cleaves after glutamine residues and no human proteases share this specificity, making it an attractive selective antiviral target. Structurally, Mpro functions as a homodimer; each monomer comprises three domains, and the substrate-binding pocket spans subsites S1′ (catalytic, Cys145/His41), S1, S2, and S3/S4. Besides the catalytic site, a potential allosteric site at the dimerization interface (between domains II and III) has been identified, involving the N-finger (Ser1–Gly11) and residues around Asn214 and the C-terminal region, suggesting that interfering with dimerization may inactivate the protease. This study aims to identify small molecules capable of dual-site inhibition—targeting both catalytic and allosteric/dimerization sites—using integrated in silico methods.

Prior studies established Mpro’s structure, substrate specificity, and active-site subsite architecture (S1′, S1, S2, S3/S4) and reported numerous active-site inhibitors, including covalent warheads and peptidomimetics (e.g., nirmatrelvir). Allosteric inhibition at the dimer interface has also been explored, with pelitinib and other ligands targeting the inter-domain site demonstrated by structural and biophysical studies. Computational campaigns (ligand- and structure-based virtual screening, docking, MD, QM/MM) have been widely applied to discover Mpro inhibitors. However, dual-site (catalytic and dimerization) inhibition remains less explored. This work builds on that literature by combining ligand-based screening, structure-based docking, multivariate analysis (PCA) relative to pelitinib, and MD simulations to prioritize dual-binding benzo[b]thiophene and benzo[b]furan derivatives.

A hierarchical, hybrid in silico workflow was applied: - Ligand-based screening (DRUDIT BPT): Constructed an Mpro catalytic site template using docking of known drugs; screened an in-house database (~10,000 heterocycles) with DRUDIT Biotarget Predictor Tool using standard parameters (N=500, Z=50, G=a). Ranked hits by Drudit Affinity Score (DAS, 0–1; higher indicates stronger predicted binding). A cutoff of DAS ≥ 0.8 selected 24 benzo[b]thiophene (series 1) and benzo[b]furan (series 2) derivatives (general scaffold: ethyl 3-benzoylamino-5-[(1H-imidazol-4-yl-methyl)-amino]-heteroaryl-2-carboxylates). - ADME filtering: Evaluated drug-likeness and medicinal chemistry alerts using SwissADME (Lipinski, Veber, Egan, Ghose; PAINS). Thirteen of 24 compounds exhibited zero rule violations; all had zero PAINS alerts; no compounds were excluded from further analysis. - Structure-based docking at catalytic site: Protein preparation for PDB 7VH8 (Mpro bound to PF-07321332, nirmatrelvir) with covalent bond to Cys145 first broken; prepared with Protein Preparation Wizard (protonation via Epik at pH 7.0±0.4; waters removed; H-bond network optimized; restrained minimization with OPLS-2005). Ligands prepared using LigPrep (all tautomers/stereoisomers; pH 7.0±0.4; Epik ionization) and minimized with OPLS-2005. Glide XP docking with grids centered on co-crystal ligands; post-docking minimization; up to two poses per conformer. Validation by re-docking reference ligands achieved RMSD 0.51 Å. Induced Fit Docking (IFD) was then performed to account for receptor flexibility (IFD score = 1.0 Glide Gscore + 0.05 Prime Energy; more negative values indicate better binding in reported tables). - Multivariate analysis (PCA) for allosteric potential: Merged molecular descriptor matrix (from DRUDIT) of the 24 compounds with those of pelitinib (non-covalent dimerization-site inhibitor). PCA (TIBCO Statistica) on the matrix (Supplementary Matrix S1) captured 50% variance in first two components. Compounds proximal to pelitinib in PC1–PC2 space were selected for allosteric docking; outliers identified via Grubb’s test were excluded. - Structure-based docking at allosteric/dimerization site: Used PDB 7AXM (Mpro complexed with pelitinib). Protein and ligand preparation as above; Glide XP docking and IFD at the allosteric site; grids centered on pelitinib. Amino acid interaction maps were analyzed within 3 Å for key residues (N-finger, Arg298, Cys300, Ser301, etc.). - Molecular dynamics (MD) simulations: MacroModel (Schrödinger) with OPLS-2005 for the two top complexes: 1d at catalytic site (7VH8), 1c at allosteric site (7AXM). Systems minimized (1000 steps) and simulated at ~300 K for 15 ns (1d) and 20 ns (1c). Protein and ligand binding energies over time were monitored; both systems showed stable plateaus (Figures S1–S2).

- DRUDIT ligand-based screening: Applying a DAS cutoff of 0.8 prioritized 24 compounds across two scaffolds (benzo[b]thiophene series 1 and benzo[b]furan series 2). Representative high DAS values (Table 1) included: 2b (0.940), 2a (0.922), 2c (0.922), 1b (0.918), 2g (0.916), 2f (0.910), 2l (0.910), 2i (0.904), 2h (0.900), 1d (0.900), 1c (0.898), among others. - ADME (SwissADME): 13/24 compounds had zero Lipinski/Veber/Egan/Ghose violations; all compounds had 0 PAINS alerts, supporting drug-likeness and medicinal chemistry friendliness. - Catalytic site docking/IFD (PDB 7VH8): Multiple compounds achieved IFD scores more favorable than nirmatrelvir (reference IFD score −673.142). Top performers (Table 3) included: 1d (−675.768; docking −8.979), 2l (−675.108; −12.040), 1j (−674.838; −7.595), 1f (−674.292; −9.781), 1i (−674.180; −9.222), 2i (−674.046; −11.050), 2h (−674.040; −10.969), 1l (−674.037; −7.673), 1a (−673.969; −8.573), 1k (−673.927; −10.744), 1c (−673.740; −8.008), 1b (−673.730; −8.314). Key interactions spanned S1′ (Cys145/His41 dyad), S1 (Phe140, Asn142, His163, His164, Glu166, His172), S2 (Met49, Tyr54, Met165, Pro168, Val186), S3/S4 (Gln189, Ala191, Gln192, Gly251). Several compounds formed more interactions than nirmatrelvir across subregions S1′, S1, S2, S3/S4 (Table 4). Overlaps for 1b, 1c, 1l, and 2l showed consistent poses, with the imidazole moiety engaging S1 and S3/S4 via H-bonds (Asn142, Gly143, His163, His164, Glu166) and the carboxyethyl moiety stabilized in S1 mimicking peptide interactions. - PCA and allosteric selection: PCA (50% variance in PC1–PC2) guided selection of compounds near pelitinib: 1b, 1c, 1g, 1h, 1i, 1l and 2a, 2b, 2c, 2g, 2h, 2i, 2l (Table 5 provides coordinates and distances to pelitinib). - Allosteric/dimerization site docking/IFD (PDB 7AXM): Several compounds exceeded pelitinib’s IFD score (pelitinib −691.09). Top IFD scores (Table 6) included: 1c (−693.48; docking −7.005), 1b (−692.66; −7.214), 2l (−692.66; −6.327), 1l (−692.59; −6.72), 2i (−692.24; −7.522), 1g (−692.05; −6.368), 1i (−691.57; −5.981), 2b (−691.36; −6.187), 2c (−691.03; −6.675). Predicted interactions at the dimer interface included residues Thr154, Pro252, Gln256, Val297, Arg298, Cys300, Val303, Thr304 and N-finger residues (Ser1, Gly2, Phe3, Arg4), suggesting interference with monomer–monomer interactions. - Dual-site candidates: Integrating catalytic and allosteric IFD results, compounds 1b, 1c, 1i, 1l and 2i, 2l emerged as promising dual-site inhibitors of Mpro with favorable ADME profiles. - MD simulations: Complexes 1d–Mpro (catalytic) and 1c–Mpro (allosteric) remained stable over 15 ns and 20 ns, respectively, with binding energies reaching plateaus, supporting robustness of docking poses.

The integrated ligand- and structure-based approach identified benzo[b]thiophene and benzo[b]furan derivatives as strong binders at Mpro’s catalytic site, often surpassing nirmatrelvir in IFD scoring, with extensive engagement of conserved residues across S1′, S1, S2, and S3/S4. Structural analyses suggest the imidazole substituent mimics the glutamine side chain of natural substrates, stabilizing interactions in S1/S3–S4, while the carboxyethyl group engages S1 akin to peptide bonds, and hydrophobic/aromatic fragments fit the flexible S2 pocket with possible π–π interactions (e.g., His41). PCA-guided selection and subsequent allosteric docking indicate that several analogues can also occupy the dimerization interface and engage key residues (N-finger, Arg298, Cys300, Ser301, etc.), potentially disrupting dimer formation and hence protease activation. The convergence of strong binding at both sites for select compounds (1b, 1c, 1i, 1l; 2i, 2l) supports the hypothesis that dual-site modulators could simultaneously inhibit proteolysis and dimerization, offering a potentially more robust antiviral mechanism less susceptible to resistance. The high conservation of interacting residues suggests possible efficacy across SARS-CoV-2 variants and related coronaviruses. MD simulations further support the stability of representative catalytic and allosteric complexes, adding confidence to the predicted binding modes.

A hybrid virtual screening workflow integrating DRUDIT ligand-based prioritization, SwissADME filtering, IFD docking at both catalytic (7VH8) and allosteric/dimerization (7AXM) sites, PCA against pelitinib descriptors, and MD simulations identified benzo[b]thiophene and benzo[b]furan derivatives as promising SARS-CoV-2 Mpro inhibitors. Notably, compounds 1b, 1c, 1i, 1l and 2i, 2l exhibited favorable ADME properties and strong predicted binding at both catalytic and allosteric sites, suggesting dual-site inhibitory potential. The study highlights a feasible strategy to discover dual-action antivirals against Mpro, a concept with potential breadth across variants due to conserved residue interactions. Future work should synthesize and experimentally validate these hits via biochemical assays (catalytic inhibition, dimerization disruption), structural biology (co-crystallography), and antiviral cell-based studies, followed by pharmacokinetic and toxicity profiling.

The study is entirely based on in silico methods (ligand-based screening, docking, PCA, and short-timescale MD) without experimental biochemical or cellular validation. Docking/IFD scoring and short MD simulations, while informative, may not fully capture entropic effects, solvent dynamics, or long-timescale conformational changes at the dimerization interface, which could affect generalizability to in vitro/in vivo outcomes.