The newly-emerged SARS-CoV-2 nucleotide gene (YP_009725309.1) was retrieved from the National Center for Biotechnology Information (NCBI) nucleotide database (https://www.ncbi.nlm.nih.gov/protein/YP_009725309.1). A homology model for the SARS-CoV-2 nsp14 was built using the Swiss Model web server (Waterhouse et al., 2018). The structure of the SARS-CoV nsp10/nsp14 complex (PDB ID: 5NFY) was the most sequelogous (99% sequence identity) to the SARS-CoV-2 nsp14 with a resolution of 3.38 Å. The quality of the homology model was analyzed using QMEAN score and Ramachandran plot generated by Swiss Model server itself (Waterhouse et al., 2018). It was further analyzed using the ProSA server (Wiederstein and Sippl, 2007). The active site pairwise structure alignment of CRN4 and homology model was carried out using TM align server and MATRAS protein 3D structure comparison tool (Kawabata, 2003; Zhang and Skolnick, 2005).

Molecular docking studies were performed using the AutoDock Vina software package. Results were analyzed using Autodock Tools and Chimera packages (Forli et al., 2016; Anil, 2019). Blind dock studies were carried out using previously identified three DEDDh/DEEDh subfamily nuclease inhibitors, pontacyl Violet 6R (PV6R), aurintricarboxylic acid (ATA) and 5,5′-dithiobis (2-nitrobenzoic acid) (DTNB) (Huang et al., 2016). The homology model of the SARS-CoV-2 nsp14 was used as the macromolecule.

MD simulations were performed using GROMACS software package. The RMSD of the PV6R molecule at the ExoN active site was calculated. Docking pose of the PV6R molecule was used to carry out the MD simulation. Ligands topology was prepared using CGenFF server and protein topology was prepared using Gromacs software itself. MD simulations were done using CHARMM36 all-atom force field. All other molecule preparations were done using Chimera software package unless otherwise stated.

Prepared complexes were immersed in a cubic box with a 1.5 nm distance between the protein surface and the boundary of the cubic box. The selected cubic box was filled with SPC/E water molecules to solvate the system and periodic boundary conditions were applied from all sides. To neutralize the protein and to maintain the net ionic strength of the system to a value closer to the physiological conditions Na⁺ and Cl⁻ ions were added replacing the solvent molecules. Energy minimization was done using the steepest descent algorithm with an energy minimization step of 0.01 kJ mol⁻¹ for 50,000 steps setting the maximum force value for 1,000 kJ mol⁻¹ nm⁻¹. NVT and NPT equilibrations were carried out for a total of 100 ps for each with 2 fs steps. Product MD run was performed for a total of 21.0 ns consist of 2 fs steps. Temperature and pressure coupling were performed using the modified Berendsen thermostat and Parrinello-Rahman method. Ewald particle mesh method was used to calculate long-range electrostatic interactions. Results were analyzed using the VMD software package.

The compound selected from molecular docking (lead compound) was virtually screened using the OpenEye Scientific vROCS application(ROCS 3.3.2.2: OpenEye Scientific Software, Santa Fe, NM. http://www.eyesopen.com.; Hawkins et al., 2007) Virtual screening was carried out using eMolecule and Kishida databases (ROCS 3.3.2.2: OpenEye Scientific Software, Santa Fe, NM. http://www.eyesopen.com). Best 20 molecules from each database were docked with homology model and average structure obtained from MD simulation using Autodock vina and the binding energies of selected compounds were calculated. Compounds that bind at the ExoN active site were selected. Control docks were carried out with endogenous ribonucleotides. All of the compounds were prepared to be optimized in their active forms in physiological conditions using Avagadro software (López, 2012). Binding of selected compounds at the ExoN active site of nsp14 was further examined using LigPlot application (Wallace et al., 1995). ADME parameters, pharmacokinetic properties, druglike nature and medicinal chemistry friendliness was evaluated using SwissADME server (Daina et al., 2017).

The lead molecule was used as the seed compound to generate novel compounds through generative shape-based neural network decoding using the LigDream web application (Skalic et al., 2019). The PV6R molecule was used as the seed compound to generate novel compounds from generative shape-based neural network decoding. Novel compounds generated by LigDream were arranged using RDock software based on the RDock score (Ruiz-Carmona et al., 2014). Five compounds with the best RDock score were directed to molecular docking using AutoDock vina for comparison purposes.

The SWISS MODEL server was used to build the homology model for the SARS-CoV-2 nsp14 using SARS-CoV nsp14/nsp10 complex (PDB ID: 5NFY) as the template. The template had a 99% sequence alignment with SARS-CoV-2 nsp14. The built homology model showed a QMEAN score of −3.18, GMQE score of 0.98, and MolProbity results were better compared to the template (Supplementary Table S1). ProSA server results showed a Z-score of -9.15, suggesting an accurate homology model. The Produced model was further evaluated using the TM align server and MATRAS pairwise 3D alignment web server shown in Supplementary Data 1 (Kawabata, 2003; Zhang and Skolnick, 2005). A 98% structure alignment (512 residues align with template out of 523 residues in the entire protein) of the homology model was obtained within a five Å range compared to the template and a root mean square deviation (RMSD) of 0.15 Å (Figure 1B). These results suggest a high-quality homology model, based on which our studies have been performed. The Ramchandran plots of the Homology model and the template are provided in Supplementary Figures S1, S2.

PV6R, ATA and DTNB have been previously identified as CRN-4 of and RNase T inhibitors. CRN-4 and RNase T are DEDDh exonucleases found in C. elegans and E. coli, respectively (Huang et al., 2016). Structure alignment of CRN-4 and nsp14 of SARS-CoV-2 using MATRAS pairwise 3D structure alignment and Chimera structure comparison shows close proximity (Figures 1C,D) of active site residues. Specifically, all DEDDh exonucleases share a conserved active site with a common set of residues (Huang et al., 2016; Ogando et al., 2019). In the nsp14 of SARS-CoV-2, residues Asp 90, Glu 92 (motif I), Glu 191 (motif II), Asp 273 (motif III) and His 268 (motif III) constitute the active site while the equivalent of E191 alternates between E and D in different nidovirus taxa. Nsp14 has Asp 90, Glu 92, Glu 191, His 268 and Asp 273 residues as the active site residues while CRN-4 having Asp15, Glu 17, Asp 115, His 179 and Asp 184 as active site residues (Huang et al., 2016; Ogando et al., 2019). This observation suggests that previously identified inhibitors for the DEDDh exonuclease may serve as effective inhibitors of nsp14.

Using the selected DEDDh exonucleases inhibitors, molecular docking revealed that PV6R binds to the ExoN binding site of nsp14 with a preferable calculated binding energy of -8.3 kcal/mol. We further analyzed the interactions between active site residues and PV6R. Figure 2A shows that the PV6R molecule is stabilized via two hydrogen bonds at the active site including a very short hydrogen bond between Ala 187 and fifth oxygen of PV6R along with hydrophobic interactions significantly contributing to the stabilization. Interestingly, most of these hydrophobic interactions are formed with five residues (Asp 90, Glu 92, Glu 191, His 268 and Asp 273) that are involved in the catalytic activity. Protonated state of PV6R shows a similar calculated binding energy (−8.3 kcal/mol) forming hydrogen bonds with Gly 93, His 268, Asp 273 and Asn 266N.

The stability of the PV6R and ExoN complex in an aqueous solution was evaluated using the parameter RMSD. MD simulation was performed for 21.0 ns using the docked structure of PV6R to ExoN. The average RMSD of the trajectories for bound protein backbone atoms showed relative stability (Figure 2B). During the simulation PV6R molecules shows a shift in binding mode with structural conformational change to achieve more stabilized interaction at the ExoN active site showing more stable calculated interaction energy in the later part of the simulation (Figures 2B–D). This flexibity is observed due to the presence of more polar groups in PV6R that can form stronger interaction with acidic residue rich active site of ExoN.

Next, virtual screening was performed using the vROCS application (OpenEye scientific software). This software package identifies the molecular features of a given molecule based on the structure, shape and electrostatics. To identify structurally and electrostatically similar molecules, the PV6R molecule was compared with the eMolecules database and Kishida database from the OpenEye software. Parameters were set to provide 20 best fits for each run. From obtained hits, molecules that bind at the active site with a calculated binding energy of less than -7.0 kcal/mol were selected and shown in Table 1. Two molecules from the eMolecules database (Molecule reference numbers: 22722115_4 and 6826969_6) and three molecules from the Kishida database (Molecule reference numbers: KS122-0741955_27, KS122-0742530_28, KS122-0742095_14) were selected for further studies as shown in Figure 3.

Figure 4 shows the molecular docking studies of these five selected compounds identified in Figure 3 (22722115_4, 6826969_6 from eMolecules database and KS122-0741955_27, KS122-0742530_28, KS122-0742095_14 from Kishida database). Each compound was further analyzed by the LigPlot application to identify the interactions that make these molecules stable at the binding site. Table 2 shows the interacting amino acids with the PV6R and the five selected compounds at the ExoN binding site. KS122-0741955_27, KS122-0742095_14 and KS122-0742530_28 exhibit a similar pattern of interactions at the active site forming one hydrogen bond with the Gly 93 residue and several hydrophobic interactions. All three molecules form hydrophobic interactions with Asp 90, Val 91, Glu 92, Pro 141, Phe 146, Phe 190, Val 191, His 268 and Asp 273 while KS122-0741955_27 and KS122-0742530_28 displays interactions with Asn 104 and KS122-0742530_28 form an extra hydrophobic interaction with Leu 149. Calculated Binding free energies of KS122-0741955_27, KS122-0742095_14 and KS122-0742530_28 were −8.2 kcal/mol, −7.9 kcal/mol and −8.7 kcal/mol, respectively.

From the calculated binding energies, KS122-0742530_28 represents the most favorable candidate due to the presence of extra hydrophobic interactions along with close H bond acceptor and donor localizations (2.95 nm). Calculated binding free energy of -8.7 kcal/mol suggests KS122-0742530_28 as a better ExoN inhibitor than PV6R. Molecules identified from eMolecules database showed higher binding free energies due to the presence of high electronegative groups and the formation of weak hydrogen bonds in 6826969_6 and the lack of strong intermolecular hydrogen bonding in 22722115_4.

Molecular docking studies between compounds identified with virtual screening and average structure of ExoN obtained from MD simulation showed further decrease for the KS122-0741955_27, KS122-0742530_28 and KS122-0742095_14 molecules with a calculated binding energy value of −8.3 kcal/mol, −9.1 kcal/mol and 8.0 kcal/mol respectively. Molecules from eMolecules database showed an increase in energy, as shown in Table 1. The drastic stability enhancement of KS122-0742530_28 can be due to its structural similarity to the PV6R molecule showing similar trend as did the PV6R molecule during the MD simulation. In general, the enhanced binding strength of molecules obtained from Kishida database shows their effectiveness and can be promising lead molecules for further experimental validation.

Next, a technique based on machine learning was used to generate novel molecules using PV6R as the seed compound. Three-dimensional shape and pharmacophoric features of the seed molecule (lead molecule) were used to produce lead-like compounds (Skalic et al., 2019). 92 identified novel compounds were further arranged using the RDock binding scores at the ExoN active site of the nsp14 (SMILES of 92 compounds have been provided in Supplementary Table S2). Five compounds (Supplementary Figure S3) with the best RDock scores were further docked using AutoDock vina. Their calculated binding energies at the ExoN site were tabulated in Table 1. Together with the virtual screening results, these identified molecules with preferable binding at the Exon active site serve as promising candidates for further investigations of their competitive inhibition of nsp14.

Physiochemical studies of the selected compounds (Table 3) showed that the calculated TPSA values were within acceptable limits. All selected compounds possess a large number of hydrogen bond acceptors and 6826969_6 has two hydrogen bond donors. The presence of moieties that enable the hydrogen bond formation increases the binding affinity at the active site of the ExoN. A higher polar surface area indicates the drug’s ability to permeate the cell membrane. Lipophilicity was calculated as an average value of iLOGP, XLOGP3, WLOGP, MLOGP and SILICOS-IT (Daina et al., 2017) indicated as log P, which can be implicated in blood-brain barrier penetration and permeability prediction. Further, log P along with the molecular weight (MW) can be used to predict the metabolism and excretion of xenobiotic compounds from the human body. All compounds possess better aqueous solubility compared to PV6R according to this analysis (Table 3).

High gastrointestinal (GI) absorption and impermeability of blood-brain barrier (BBB) make these promising candidates to be used as drugs against respiratory illness (Table 4). The pharmacokinetic parameters suggest that these compounds follow drug-likeness rules (Lipinski, Ghose, Veber, Egan and Muegge) confirming that these molecules are drug-like and have a strong possibility to be directed for further studies to be used as drugs (Ghose et al., 1999; Egan et al., 2000; Muegge et al., 2001; Veber et al., 2002; Lipinski et al., 2012). Two medicinal chemistry filters, pan assay interference compounds (PAINS) and Brenk, showed 0 alerts for four of the selected compounds, suggesting their lead likeness.

In summary, our molecular docking and molecular dynamics studies show that PV6R binds at the ExoN catalytic active site of SARS-CoV-2 with a strong calculated free binding energy of -8.3 kcal/mol. Virtual screening using PV6R as the lead molecule identified KS122-0742530_28 from the Kishida database with a better binding free energy. Both PV6R and KS122-0742530_28 confirmed to have drug-like physiochemical and pharmacokinetic properties. These molecules may serve as lead molecules in further experimental validation for the development of ExoN inhibitors against SARS-CoV-2, including profiling their 50% inhibitory concentration (IC50) values and virus-inhibitory effects.