SARS-CoV-2 nsp14 is a multidomain protein. The N-terminal domain functions as proofreading exoribonuclease, and the C-terminal is a methyltransferase. SARS-CoV-2 nsp14 shares 95.07% amino acid sequence identity (over complete protein length) with SARS-CoV nsp14 (Supplementary Figure 1). ExoN domain of SARS-CoV nsp14 resembles DEDD-type ExoNs (Ma et al., 2015). The DEDD superfamily members are defined by the presence of three canonical motifs—DXE (motif I), W(X)4EL (motif II), and DAIMTR (motif III) (Shannon et al., 2020). The presence of DEED instead of DEDD and an additional H makes the SARS-CoV ExoN a DEEDh-type ExoN (Ogando et al., 2019). In SARS-CoV-2, the catalytic residues—Asp90, Glu92, Glu191, His268, and Asp273, and the canonical motifs are conserved (Supplementary Figure 1). A 3-dimensional (3D) model of SARS-CoV-2 nsp14 was built using SARS-CoV nsp14 (PDB ID: 5C8S) as a template. A grid comprising the three conserved motifs was used for docking.

Ten thousand three hundred ninety-seven conformers generated from 2,240 approved small molecule drugs were screened using AutoDock VINA. Based on binding free energy, the top 20 binding poses were selected for further analysis (Figure 2 and Table 1). All 20 poses interact with catalytic residues. Dexamethasone metasulfobenzoate binds to the catalytic site of ExoN with the binding energy of −8.7 kcal/mol. Conivaptan, dutasteride, hesperidin, lumacaftor, and glycyrrhizic acid bind ExoN active site with the slightly higher energy of −8.6 kcal/mol. Interaction of ExoN domain with 12 unique drug molecules, corresponding to top 20 poses, was studied and is depicted in Table 2. Most of the analyzed poses interact with at least three of the five catalytic residues (Figures 3, 4).

Hesperidin is the only drug that interacts with all five catalytic residues. Hesperidin and glycyrrhizic acid have four and three ExoN catalytic residues, respectively, within hydrogen-bonding distance (Figure 4). The binding of these drugs to ExoN catalytic residues might potentially cause obstruction of substrate binding and catalysis.

Remdesivir, an investigational drug for the treatment of Ebola, was shown to inhibit SARS-CoV-2 RdRp, and inhibiting RNA synthesis (Gordon et al., 2020). As remdesivir is a drug without anti-ExoN activity, it was used here as a negative control. Remdesivir displayed a low binding affinity toward ExoN (−6.0 kcal/mol) than the top 20 poses (<−8.4 kcal/mol). The estimated binding energy of remdesivir with ExoN is -6.0 kcal/mol, higher than the top 20 poses from the approved drug category (Table 1).

Molecular dynamics simulation can provide atomistic insights on structural stability and the dynamic of protein–ligand interactions (Luthra et al., 2009; Prakash and Luthra, 2012; Wang et al., 2013; Panda et al., 2020). Based on molecular docking, interactions with active site ExoN residues and their antiviral properties, conivaptan, hesperidin, and glycyrrhizic acid were chosen for MD studies. As astemizole was shown to inhibit SARS-CoV-2 in in vitro assays, it was included in the MD studies (Riva et al., 2020).

The structural dynamics of glycyrrhizic acid, astemizole, conivaptan, and hesperidin in complex with ExoN displays maximum population density of stable conformation at ∼6.0, 6.5, 8, and 6 Å, respectively, relative to ExoN, which equilibrated at around 9.75 Å. Hence, drug molecules induced substantial rigidification in ExoN structure (Figure 5A). ExoN–glycyrrhizic acid exhibited the least structural fluctuations, suggesting the most stable protein–ligand complex. Although the complex of ExoN–conivaptan achieved a maximum population density of around 8 Å, the population density of conformational dynamics ranges from ∼4.0 to 9 Å. The ExoN–conivaptan complex shows a slightly smaller peak at ∼6.0 Å too. It suggests conivaptan might move between two conformations. The structure of ExoN and ExoN–glycyrrhizic acid, astemizole, conivaptan, and hesperidin had a maximum population density of radius of gyration (RoG) around 33, 33.5, 31.5, 32.2, and 32.2 Å, respectively (Figure 5B). During the simulation period of 200 ns, all five systems were stable around the solvent-accessible surface area (SASA) values of 2,700 to 2,900 Ų. RoG and SASA results suggest marginal or no structural compactness change of ExoN and ExoN–drug complexes (Figure 5C).

To understand the drifts in root mean square deviation (RMSD) plots (Figure 5 and Supplementary Figure 2A), the average distance of the four drug molecules from the center of the ExoN active site was measured. The time evolution distance plots show that the average distance of hesperidin and conivaptan remained consistent between 3.5 and 4.5 Å from the active site of ExoN (Supplementary Figure 3). Glycyrrhizic acid and astemizole move out from the binding pocket around ∼50 and ∼100 ns of simulation, respectively. The conformational adaptability of hesperidin and conivaptan during the simulation was explored by performing root mean square fluctuation (RMSF) analyses. The average RMSF peaks of all the amino acids of ExoN-hesperidin and ExoN–conivaptan complex are less than ExoN (Supplementary Figure 2 and Supplementary Methods). The RMSF values provide structural evidence of stable molecular interaction of hesperidin and conivaptan with ExoN. The average distances between metal ions (Mg²⁺) remain around 3.6 Å (Supplementary Figures 3, 4).

The efficacy of a drug molecule is largely dependent on molecular interactions at the active site (Schiebel et al., 2018; Mishra et al., 2021), and the network of H-bonds plays a crucial role in this interaction. Thus, H-bond interactions between ExoN and drug molecules were calculated (Chen et al., 2016). Distance cutoff of 3.5 Å, and angle cutoff of 135° were used for the calculation of H-bonds. Maximum occupancy of seven H-bonds between ExoN and hesperidin was observed, of which five to six H-bonds were observed consistently during the simulation (Supplementary Figure 6). The molecular interaction of conivaptan with ExoN shows the maximum possibility of two H-bonds. Out of that, only one H-bond remains consistent throughout the simulation time (0–200 ns).

To ascertain the molecular binding interaction of hesperidin and conivaptan with ExoN, a quantitative assessment of binding free energy (ΔG_binding) was carried out using MM-GBSA (Genheden and Ryde, 2015) on the conformational ensemble of protein–ligand complexes. Hesperidin shows more favorable binding free energy, ΔG_binding = −119.07 ± 0.69 kcal/mol, as compared to conivaptan (ΔG_binding = −85.86 ± 0.68 kcal/mol) (Table 3).

Free energy decomposition per residue at ExoN active site indicates energetically favorable molecular binding of hesperidin and conivaptan, largely contributed by the residues involved in van der Waals and electrostatic interactions (Figure 6). Hesperidin shows energetically favorable binding to catalytic residues Val91 and Hie268 and other active site residues Trp186, Ala187, Asn252, Leu253, Gln254, and Asn266. Interactions of these residues with hesperidin were observed in the molecular docking studies too (Figure 4A). Conivaptan interacts with catalytic residues Phe146, Phe190, and Glu191 and active site residues Val91, Gln145, Phe146, Trp186, Ala187, and Phe190. Like hesperidin, the interaction of the aforementioned residues with conivaptan was observed in molecular docking (Figure 4B).

The amino acid sequence of SARS-CoV-2 ExoN/nsp14 (P0DTD1) was used as a target sequence to build a 3D model using the SWISS-MODEL web server (Schwede et al., 2003). Nsp14 from SARS-CoV (PDB ID: 5C8S, chain B) had a high sequence identity to target sequence—95.07%. This 3D structure was used as a template for building the SARS-CoV-2 nsp14 model. The model quality was assessed using the QMEAN score. The QMEAN score for the generated model was −3.14. Scores greater than −4 indicate a good model.

Two-dImensional structures of approved small molecule drugs were downloaded from DrugBank in sdf format (2,454 structures) (Wishart et al., 2006). The structures were converted to 3D format using the OpenBabel –gen3d option (O’Boyle et al., 2011; Yoshikawa and Hutchison, 2019). Few of the drug molecule structures showed error at this stage. 3D structures of these molecules were downloaded from other sources (Supplementary Table 1). Multiple conformers for each 3D structure were generated using OpenBabel to increase the conformational space sampling.

The docking grid was defined to encompass the conserved motives of the ExoN domain (see section “SARS-CoV-2 ExoN Domain”). A total of 10,397 conformers from 2,240 compounds were subjected to virtual screening using the AutoDock VINA tool on SARS-CoV-2 nsp14 (Trott and Olson, 2010). Twenty poses were generated for each compound. Compound-nsp14 interactions were visualized using LigPlot⁺ v2.1 (Laskowski and Swindells, 2011).

SARS-CoV-2 ExoN crystal structure (PDB ID: 7MC6, length 291 amino acids) became available after the completion of our work. Pairwise structure comparison of SARS-CoV-2 ExoN crystal structure and homology model (this study) was done using the Dali web server (Holm, 2020). The comparison revealed the two structures to be highly similar (DALI Z score 39.7; RMSD 0.8 over 285 aligned residues).

The SARS-CoV-2 structure contains one Mg²⁺ ion in the catalytic site (Moeller et al., 2021). ExoN uses two metal ions to remove misincorporated nucleotides. In the structures of SARS-CoV-2 and SARS-CoV, ExoN coordinates of only one Mg²⁺ ion are observed. Lassa virus NP ExoN, an ExoN of DEDDh-family, contains two Mn²⁺ ions in its catalytic site. The two Mg²⁺ ions were modeled in the docked complexes based on the Mn2+ ion position in the Lassa NP ExoN-RNA complex (PDB ID: 4GV9). Simulations were performed using the Mg²⁺ ions containing docked structures.

All-atoms MD simulation was performed using the Amber16 with force field ff14SB (Case et al., 2005; Maier et al., 2015) for the metal ion–containing protein, ExoN, and the docked complexes with drug molecules, astemizole, conivaptan, hesperidin, and glycyrrhizic acid. Antechamber (Wang et al., 2006) is used to parameterize all selected ligands using the GAFFs force field (Wang et al., 2004). Divalent Mg⁺⁺ interactions are modeled using the 12-6-4 model compatible with the TIP3P water model. For the ligand preparation, topology, and parameter files were generated using the leap module of Amber (Maier et al., 2015). Keeping the protein at the center, a cubic box is prepared with 10 Å, padding the explicit TIP3P water molecules in all directions (Jorgensen et al., 1983), and the counter-ions (Na⁺Cl^–) added to neutralize the simulation box. Particle mesh Ewald approach (Essmann et al., 1995) was used for electrostatic interaction calculation and the SHAKE algorithm (Ryckaert et al., 1977) was used to constrain H-bonds. Energy minimization of prepared systems was performed in three stages, each of 10,000 steps of steepest descent (SD) and conjugate gradient (CG) to relax the system. Furthermore, each simulation system was gradually heated from 50 to 300 K in six steps, followed by 10,000 steps of SD and CG minimization, respectively. Under the NVT ensemble condition, each system is equilibrated for 1 ns. Finally, all five systems were submitted for the production run under NPT ensemble condition for 200 ns with a time step of 2 fs.

The obtained MD trajectories were analyzed for the structural stability of ExoN and binding with drug molecules through the RMSD, RMSF, RoG, SASA, and H-bond interactions between protein and ligands during the simulation. The Cα-backbone RMSD is calculated with reference to the starting structure of the protein. RMSF defines the average positional fluctuations of the protein residues from their initial position, which is important to determine the local dynamics of a protein. The binding free energy for each complex drug molecule was estimated using the molecular mechanics combined with the generalized Born MM-GBSA along with the weighted interactions active-site residues (Miller et al., 2012). The binding free energy components can be represented according to the equations;

where, ΔE_MM represents the enthalpic components, whereas ΔE_sol represents the polar and nonpolar electrostatic components from solvation. Here, the polar electrostatic component is calculated using the GB model, whereas the nonpolar electrostatic contribution is calculated by SASA. The last 50-ns simulation trajectory is used, which was sampled per 10-ps interval.