Virtual screening and molecular dynamics simulation for identification of natural antiviral agents targeting SARS-CoV-2 NSP10.

New variations of SARS-CoV-2 continue to emerge in the global pandemic, which may be resistant to at least some vaccines in COVID-19, indicating that drug and vaccine development must be continuously strengthened. NSP10 plays an essential role in SARS-CoV-2 viral life cycle. It stimulates the enzymatic activities of NSP14-ExoN and NSP16-O-MTase by the formation of NSP10/NSP14 and NSP10/NSP16 complexes. Inhibiting NSP10 can block the binding of NSP10 to NSP14 and NSP16. This study has identified potential natural NSP10 inhibitors from ZINC database. The protein druggable pocket was identified for screening candidates. Molecular docking of the selected compounds was performed and MM-GBSA binding energy was calculated. After ADMET assessment, 4 hits were obtained for favorable druggability. The analysis of site interactions suggested that the hits all had excellent binding. Molecular dynamics studies revealed that selected natural compounds stably bind to NSP10. These compounds were identified as potential leads against NSP10 for the development of strategies to combat SARS-CoV-2 replication and could serve as the basis for further studies.

Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the pathogen of coronavirus disease 2019 (COVID-19), which has caused and still causes great health and economic burden on the world [1]. In addition, the SARS-CoV-2 variants (i.e., Delta and Omicron) have been proposed significant challenges to human health worldwide [2]. As of May 28, 2022, this deadly virus has claimed more than 6 million lives (https://www.worldometers.info/coronavirus/). The progress in vaccine and drug development is encouraging in terms of preventing the SARS-CoV-2 pandemic [3]. The continuous mutation of the virus increases the risk of SARS-CoV-2 epidemics, which urgently needs to be addressed through the development of innovative and effective strategies for preventive and therapeutic drugs.

SARS-CoV-2 is a single-stranded positive-sense RNA virus, belongs to the genera Betacoronavirus from the order Nidovirales [4]. The replication of SARS-CoV-2 is strongly dependent on the non-structural proteins (NSPs), which are encoded by two overlapping open reading frames (ORFs), ORF1a and ORF1b [5]. These two ORFs are translated into two polyproteins (pp1a and pp1ab), and pp1a is autocatalytically processed to NSP1-11 and pp1ab is processed to NSP1-16 by viral proteases 3CLpro and PLpro [6]. The various NSPs have their own structures and functions, and they are responsible for viral replication and transcription, suppression of host innate immunity, and modulating other cellular functions [7].

Among these NSPs, NSP14 is a bifunctional enzyme, having RNA cap guanine N7-methyltransferase (N7-MTase) and 3′–5′ exoribonuclease (ExoN) activities, and is necessary for viral replication and transcription [8]. During coronavirus replication, the proofreading of NSP14-ExoN can prevent lethal mutagenesis in RNA-dependent RNA polymerases (RdRps) replication [9]. Although the NSP14–N7-MTase is independent of NSP10, the activity of NSP14-ExoN is highly regulated by NSP10 [10,11]. In order to escape from the degradation by host nucleases, SARS-CoV-2 has modified its genome at the 5′-untranslated region by adding a cap [12]. NSP16 is an S-adenosyl-L-methionine (SAM)-dependent 2′-O-methyltransferase (2′-O-MTase), and plays central role in RNA cap methylation and directly related to immune escape [13]. It catalyzes the ribose 2′-OH methylation of viral RNA cap structures, which prevents virus detection by cell innate immunity mechanisms [14]. In a word, the methylation of RNA cap at ribose 2′-O is catalyzed by NSP16-2′-O-MTase, and the proofreading role for prevention of lethal mutagenesis is catalyzed by NSP14-ExoN. It has been proved that the activities of both NSP14-ExoN and NSP16-2′-O-MTase are strongly dependent on stimulation of NSP10 [15]. So NSP10 plays an important role in minus-strand RNA synthesis of related coronavirus.

As noted, NSP10 stimulates the enzymatic activities of NSP14-ExoN and NSP16-2′-O-MTase through protein-protein interactions [16]. The NSP10 can increase the NSP14-ExoN activity by more than 35 times. The NSP16 turns the cap-0 into a cap-1 by methylation with the assistance of NSP10. NSP10 is a bifunctional protein essential for SARS-CoV-2 replication [17]. It has previously been shown that breaking the interaction of NSP10 with NSP14 and/or NSP16 can block the viral replication cycle. Furthermore, these viral proteins are also involved in the dysregulation of antiviral IFN signaling [18]. It has been reported that SARS-CoV replication can be successful inhibited by TP29, a peptide inhibitor targeting NSP10 [19]. NSP10 plays a critical role not only in viral replication but also in the regulation of the innate immune response, and it is a potential target for COVID-19 treatment [11].

The present study aims to identify natural candidates as inhibitors of NSP10 using virtual screening approach in silico. The natural product library was downloaded from ZINC database and preliminary screened, the selected compounds were docked to NSP10, and the candidate compounds with promising pharmacological properties were obtained by ADMET analysis. Finally, molecular dynamics simulations were performed to assess the stability of complexes of NSP10 with the hits. Results reveal that all four selected compounds exhibit well binding affinity to NSP10, and are potential SARS-CoV-2 antagonists.

Materials and methods

Preparation of protein and compound database

The crystal structure of SARS-CoV-2 NSP10 (PDB ID: 6ZPE) [20] was downloaded from RCSB protein data bank (https://www.rcsb.org/) and imported into the Protein Preparation Wizard in Maestro (Schrödinger, LLC, New York, NY, 2020) [21]. It was added hydrogens, as well as adjusted the protonation states of histidine residues, removed waters based on their hydrogen bonding pattern and performed restrained minimization using OPLS3e force field [22].

The 224,205 natural compounds were downloaded from the ZINC database (https://zinc.docking.org/), their structures were optimized using OPLS3e force field with LigPrep [23]. The ionization states and the tautomers were both generated under the pH 7.0 ± 2.0 using Epik [24], with a maximum of 32 isomers obtained for each natural compound.

Binding pocket prediction

NSP10 potential binding pocket was predicted using DoGSiteScorer (DGSS) from Protein Plus server (https://proteins.plus/) and the druggability was detected by Difference of Gaussian (DoG) filter. The pocket with detected score greater than 0.4 is considered druggable [25].

NSP10 potential binding pocket was also validated by Protein Allosteric Sites Server (PASSer) (https://passer.smu.edu/). Based on FPocket algorithm and ensemble learning method, the allosteric pockets were also identified and detected [26]. These allosteric pockets were calculated as average probability of eXtreme gradient boosting (XGBoost) and graph convolutional neural network (GCNN) models. The top ranked NSP10 pocket was considered as potential binding site.

Screening the compound database

Based on the determined binding site, virtual screening of natural compounds from the ZINC database was performed using Glide module (Schrödinger, LLC, New York, NY, 2020) [27]. The prepared database was prefiltered using LigFilter module to remove molecules that were not in accordance with the drug-likeness properties. Firstly, Lipinski's Rule of Five (RO5) [28] was used to filter out compounds that had unsuitable properties and reactive functional groups. Then, the obtained compounds were initially screened using the predicted binding pocket, and sequenced by Glide Score (GScore). Subsequently, standard precision (SP) and extra precision (XP) docking were employed for precise screening [29]. Docking flexibility and Epik-based penalty counts in docking mode were set as default.

Binding free energy calculation

The binding free energies (ΔGbind) of XP docking were calculated by the Prime module of Schrödinger [30]. The structures were minimized for binding energy calculation through molecular mechanics-generalized Born surface area (MM-GBSA) approach [31]. The variable-dielectirc generalized Born (VSGB) solvation model [32] was used for polar solvation and the OPLS3e force field was used for energy calculation.

The binding free energies (ΔGbind) were calculated using following equation.here, ΔEMM is molecular mechanics energy, including electrostatic energy (ΔEelec) and van der Waals energy (ΔEvdW), ΔGpolar is polar solvation, ΔGnonpolar is non-polar solvation, and TΔS represents the conformational entropy.

ADMET prediction

The physicochemical and ADMET properties of the selected compounds, such as molecular weight, topological polar surface area (TPSA), hERG inhibition, AMES toxicity and gastrointestinal absorption (GIA), were predicted and analyzed. SMILES notations of the hits were uploaded to the ADMETlab 2.0 server (https://admetmesh.scbdd.com/) to estimate the toxicological properties. The pharmacokinetic parameters, such as absorption, distribution, metabolism, and excretion were evaluated through the SwissADME server (http://www.swissadme.ch/).

Molecular dynamics simulation

Molecular dynamics (MD) simulations were performed using GROMACS 5.1.4 with GROMOS96 53a6 force field [33]. The topology file of NSP10 was generated by GROMACS, and the ligand topology files were created by PRODRG web server (http://davapc1.bioch.dundee.ac.uk/cgi-bin/prodrg/) [34]. Each system was solvated in a dodecahedral box with SPCE model, and neutralized by adding sodium or chloride ions, the energy was minimized using steepest descent method. Then system was gradually heated up to 310 K with NVT ensemble and equilibrated at 1 atm and NPT ensemble for 1 ns. And MD simulation was carried out using a time step of 2 fs for 100 ns. Eventually, root mean square deviation (RMSD), root mean square fluctuation (RMSF) and radius of gyration (Rg) were calculated based on MD trajectories.

Interfacial residues studies

The key residues in binding interfaces of NSP10/NSP14 and NSP10/NSP16 complexes were analyzed using EMBL-EBI's online software Interfaces and Assemblies (http://www.ebi.ac.uk/msd-srv/prot_int/pistart.html). The conformational changes of NSP10 alone compared to that in NSP10/NSP14 and NSP10/NSP16 were analyzed as well as after binding inhibitor, especially the changes near the binding interfaces.

Results and discussion

Druggable pocket identification

As shown in Fig. 1 , six binding sites were predicted by DGSS, they distributed quite evenly around the protein surface, and ranked by Drug Scores. From the pocket parameters (Table S1), it can be seen that P0 has the largest volume, surface area, the highest Drug Score, and is a well enclosed hydrophobic site, is the most druggable. Its structural elements include helix α2, α3, α4, a long C-terminal loop and α3/β3 intervening loop etc.. This binding pocket was mainly consisted of Leu75, Tyr76, Pro84, Asn85, Phe89, Asp91, Leu92, Asn114, Thr115, Val116 (Fig. S1).

Fig. 1

Predicted druggable pockets created by DoGSiteScorer. (a) The pockets shown in mesh. (b) The pocket surface represented in translucent model.

The potential binding pockets were further identified using PASSer tool and nine allosteric sites were obtained. Based on the prediction probabilities, the top three binding pockets were analyzed (Fig. S2), and their interacting residues were analyzed (Table S2). Site 1 is composed of Leu75, Tyr76, Pro84, Asn85, Phe89, Asp91, Leu92, Thr111, Leu112, Asn114, Val116, and has the highest probability. It completely overlaps the potential druggable pocket P0 provided by DGSS. In addition, Sites 2, 3, 6, 7 and 8 have similarity to pockets P4, P2, P1, P5 and P3, respectively. P0 (Site 1) is the highest ranked binding site predicted by both methods and is discussed in this study.

Virtual screening and binding energy calculation

First, the natural product database was filtered by Lipinski's rules and removed the compounds that have reactive groups, 172,009 compounds were left. The predicted binding site, P0, was used for HTVS, the top 30% (51,603) were exported for SP docking, and the top 10% (5,160) were taken to XP docking. After these, 844 compounds were identified after calculation of Epik's state penalty. Eventually, based on the XP docking results, the defined dock cutoff score is considered to be −8.0 kcal/mol [35]. So 34 natural molecules with XP GScore ≤ −8.0 kcal/mol were selected as potential inhibitors (Fig. S3), and their scores were between −10.5 and −8.0 kcal/mol and glide energies ranging from −37.6 to −76.3 kcal/mol (Table 1 and Table S3), and they were for further analyzation by MM-GBSA.

Table 1

The docking scores and MM-GBSA of top 15 molecules.

No.	ZINC ID	XP GScore (kcal/mol)	glide energy (kcal/mol)	MM-GBSA ΔGbind (kcal/mol)
1	ZINC000085489178	−8.854	−51.025	−103.05
2	ZINC000085626263	−8.743	−64.415	−97.65
3	ZINC000085488238	−10.500	−45.146	−95.82
4	ZINC000085488163	−9.250	−46.016	−93.09
5	ZINC000085488189	−10.354	−53.364	−92.04
6	ZINC000095486258	−8.549	−41.364	−91.32
7	ZINC000085488272	−9.246	−45.361	−89.98
8	ZINC000085488307	−9.657	−40.624	−89.31
9	ZINC000085626242	−8.834	−63.638	−86.42
10	ZINC000014694403	−8.117	−51.955	−86.02
11	ZINC000085597458	−8.096	−61.919	−85.78
12	ZINC000095918985	−8.701	−57.053	−77.13
13	ZINC000085488288	−8.421	−37.648	−76.41
14	ZINC000085490847	−8.113	−55.576	−72.18
15	ZINC000248252675	−8.279	−52.026	−71.93

From MM-GBSA binding free energies (ΔGbind, Table 1), it can be seen that ZINC000085489178 (1) has the lowest binding energy (−103.05 kcal/mol), followed by ZINC000085626263 (2, −97.65 kcal/mol), and ZINC000085488238 (3, −95.82 kcal/mol), ZINC000085488163 (4, −93.09 kcal/mol), ZINC000085488189 (5, −92.04 kcal/mol), ZINC000095486258 (6, −91.32 kcal/mol), and there are 11 compounds with ΔGbind lower than −80.00 kcal/mol and have high affinity to NSP10.

ADMET properties

From the toxicities predicted by ADMETlab of above 11 compounds, including AMES Toxicity, hERG Blocker, Oral Rat Acute Toxicity (LD50), Carcinogencity (Table 2 ), it is shown that all of the compounds showed no or less toxicity, except 2, 6 and ZINC000085626242 (9). All of them do not inhibit human ether-a-go-go related gene (hERG), that is, there is no risk of cardiotoxicity caused by hERG inhibition, and they all have high LD50. None of them displays positive result for carcinogenicity, and they are all safety.

Table 2

Toxicity parameters of the selected compounds.

No.	AMES Toxicity	hERG Blocker	Oral Rat Acute Toxicity (LD50)	Carcinogencity
1	(- - -)	(−)	(- -)	(−)
2	(+)	(- - -)	(- - -)	(- - -)
3	(- - -)	(- - -)	(- - -)	(- - -)
4	(- - -)	(- -)	(- - -)	(- - -)
5	(- - -)	(- -)	(- - -)	(- - -)
6	(+)	(- - -)	(- - -)	(- - -)
7	(- - -)	(- -)	(- - -)	(- - -)
8	(- - -)	(- -)	(- - -)	(- - -)
9	(+)	(- - -)	(- - -)	(- - -)
10	(−)	(- - -)	(- - -)	(- -)
11	(−)	(−)	(−)	(- -)

Note: The prediction probability values are expressed as six symbols: 0–0.1(---), 0.1–0.3 (--), 0.3–0.5 (−), 0.5–0.7 (+), 0.7–0.9 (++), and 0.9–1.0 (+++).

The predicted ADME properties (Table 3 ) demonstrate that all of the hits have acceptable lipid solubility (0 < logP <5), TPSA values of 6 and 9 are larger than 140 Å2, indicating that they had bad molecular permeability. Among them, 8 compounds display high human gastrointestinal absorption, and the exceptions are 6, 9 and ZINC000085597458 (11). Except for 2, the other compounds exhibit substrates of P-glycoprotein, meaning that they can penetrate cell membranes via the ATP-binding cassette (ABC) translocator. Similarly, none of them exhibits central nervous system (CNS) and blood-brain barrier (BBB) permeability, suggesting that they have less CNS side effects. Therefore, due to their inappropriate properties, compounds 2, 6, 9 and 11 were excluded from further discussion. The others are expected to have good druggability. Based on ADMET properties and their binding free energies, compounds 1, 3, 4, 5 were picked out for the next study.

Table 3

ADME properties of the selected compounds.

No.	MWa	HBAb	HBDc	Fraction Csp3	TPSAd	cLogP	CNSe	BBBf	P-gp substrateg	GIAh
1	470.69	5	4	0.62	95.58	4.85	−2	No	Yes	High
2	458.54	7	5	0.50	127.45	3.26	−2	No	No	High
3	412.61	6	4	0.96	93.03	2.59	−2	No	Yes	High
4	442.63	7	6	0.96	122.05	1.87	−2	No	Yes	High
5	470.69	7	6	0.96	122.05	2.55	−2	No	Yes	High
6	418.39	9	7	0.38	167.91	0.03	−2	No	No	Low
7	444.65	7	6	0.96	122.05	2.43	−2	No	Yes	High
8	372.54	6	5	0.95	101.82	2.08	−2	No	Yes	High
9	420.45	8	6	0.41	147.68	1.55	−2	No	Yes	Low
10	348.39	6	5	0.37	110.38	1.92	−2	No	Yes	High
11	482.57	7	5	0.39	127.45	3.45	−2	No	Yes	Low

Molecular weight.

Number of hydrogen bond acceptors.

Number of hydrogen bond donors.

Topological polar surface area in Å2.

CNS permeant, −2 indicates low, central nervous system penetration.

BBB permeant, No, low blood-brain barrier penetration.

P-gp substrate, Yes, substrate, No, non-substrate.

Gastrointestinal absorption, High, good absorption, Low, low absorption.

NSP10-ligand interaction

The 2D and 3D diagrams of molecular interactions between NSP10 with 1, 3, 4, 5 are illustrated in Fig. 2 . 1 forms three hydrogen bonds with NSP10, the amino group (NH2) as a donor with Asp91, the phenolic hydroxyl as a donor with Lys113 and as an acceptor with Thr115. Furthermore, the phenyl ring forms a π-cation interaction with Lys28, and imino, carboxyl, and amino involve in non-covalent interactions by forming salt bridges with residues Asp22, Lys28 and Asp91.

Fig. 2

2D and 3D interaction diagrams of the binding poses of the compounds within the NSP10 binding pocket. (a, b) 1-NSP10, (c, d) 3-NSP10, (e, f) 4-NSP10, (g, h) 5-NSP10. In the 3D diagrams, the compound shown in ball-and-stick representation (yellow carbons), NSP10 backbone presented as cartoon (light blue), and interacting amino acids presented as sticks (gray). Dashed black lines displayed hydrogen bonds, and π-cation bonds are shown as pink dashes. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

3 forms six hydrogen bonds with NSP10, 3-hydroxyl group as a donor with Asp22, NH as a donor with Asp22 and Tyr76, and the octahydro-quinolizinyl nitrogen as a donor with Asp91, the octahydro-quinolizinyl hydroxyl as a donor with Asp91 and as an acceptor with Leu92. Asp22, Lys28 and Asp91 involve in the salt bridge formation. This reveals that 3 has high binding affinity with NSP10.

4 forms five hydrogen bonds with NSP10, NH of 4 as a donor with Asp22, the oxygen of the carbonyl group as an acceptor with Lys28, 3-hydroxyl group as an acceptor with the hydrophobic Tyr76, the piperidyl hydroxyl as an acceptor with the hydrophobic Val116, and aliphatic hydroxyl group as a donor exhibits H-bond interactions with Leu112. The protonated nitrogen of secondary amine and the terminal oxygen of fatty acid in 4 establish salt bridges with residues Asp22 and Lys28, respectively. Moreover, Asp91 forms salt bridge with the protonated nitrogen atom of the piperidine ring of 4.

5 forms six hydrogen bonds with NSP10, 3-hydroxyl group of 5 as a donor and NH as a donor with Asp22, aliphatic hydroxyl group which reaching deep into the binding cavity as an acceptor with the hydrophobic Leu75 and as a donor with Thr111, and the piperidyl nitrogen as a donor with Asp91, the piperidyl hydroxyl as an acceptor with the hydrophobic Val116. The three salt bridges form by Asp22, Lys28 and Asp91 facilitate binding with the hydrophobic binding cavity.

Notably, the acidic amino acid residues, i.e., Asp22 and Asp91 can be observed to form hydrogen bonds or salt bridges in all four complexes. Other amino acid residues, such as Lsy28, Tyr76 and Val116, usually forming interactions with the compounds, are also potential key residues. Results of molecular docking and interaction analysis indicate that hydrogen bonding, hydrophobic and salt-bridge interactions are determinant forces in the formation of complexes. These interactions stabilized the binding of the molecule to NSP10, meaning that the four compounds have a high potential for NSP10 inhibition. The top four hits were selected for MD simulation to further verify the stability of the docking.

The result reveals that RMSDs (for α-C of protein backbone, Fig. 3 a) rise continuously at the beginning of the simulation, and gradually become equilibrium for all complexes. For 1-NSP10, the RMSD is in equilibrium after 10 ns, then remains balanced, and reaches to the average around 0.357 nm until the end. The RMSD of 3-NSP10 is balanced at 12 ns and then fluctuates around 0.332 nm. It is observed that 1-NSP10 and 3-NSP10 exhibit small RMSD fluctuations, which indicates both systems are highly stable. In case of 4-NSP10, the RMSD increases during the first 5 ns and then maintains at around 0.300 nm. At 28 ns, it continues to increase to 0.480 nm. Although the RMSD fluctuates slightly till 86 ns, it stabilizes at 0.490 nm later. The RMSD fluctuation is due to the flexible loop (79–92) and the loop between α1 and α2 (19–22). For 5-NSP10, the RMSD increases until 47 ns, then reaches equilibrium with fluctuation around 0.465 nm. The RMSD trend for 5-NSP10 is consistent with that of 4-NSP10. Compared to 4-NSP10, 5-NSP10 possesses higher protein stability. All RMSDs are with little fluctuations after reaching balanced, except 4-NSP10. The RMSDs for the Cα of the four systems are less than 0.6 nm, suggesting that the systems are very stable.

Fig. 3

MD simulations results of the complexes of compound 1 (yellow), 3 (orange), 4 (cyan) and 5 (brown) with NSP10. RMSD (a), RMSF (b) and Rg (c) for the Cα atoms of NSP10 complexes. (For interpretation of the references to colour in this figure legend, the reader is referred to the Web version of this article.)

The RMSF shows similar tendency (Fig. 3b). The residues 30–36 and 83–92 in 1-NSP10 have high flexibility. The fluctuation of 3-NSP10 appears at residues 19, 45–54 and 83–92. In the case of 4-NSP10, fluctuation is significantly high around residues 19–28, 32–35, 79–92 and 116–123. Similarly, RMSF of 5-NSP10 exhibits large fluctuation at residues 16–28, 30–38, 79–92 and 99–107. Furthermore, the residues in the loop (79–92) close to the ligand exhibit greater fluctuation and flexibility. In particular, 4-NSP10 flexible loop exists with fluctuation exceeding 0.5 nm. Other residues with higher peaks suggest their higher flexibility and participation in ligand binding. The results show that the binding to NSP10 dramatically increases flexibility of NSP10/NSP14 or NSP10/NSP16 binding interface. RMSF shows that the terminal part of each complex is flexible. Maximum RMSFs are 0.355, 0.419, 0.512 and 0.404 nm for the protein backbone atoms in the four systems, respectively. The overall RMSF results reveal that these complexes are sufficiently stable.

The Rg values (Fig. 3c) decrease gradually in all systems, except 5-NSP10, showing that protein structure is compressed when binding to compounds 1, 3, or 4. The Rg of 1-NSP10 decrease from 1.470 nm to 1.355 nm. The Rg average of 3-NSP10 is 1.433 nm, reflecting the system has lower degree of fluctuation. Similarly, Rg of 4-NSP10 is in the range from 1.357 to 1.472 nm, the average is 1.419 nm, and that of 5-NSP10 complex ranges from 1.400 to 1.521 nm, with mean of 1.439 nm, and the last 15 ns, Rg of the protein is stable at around 1.500 nm. The results of Rg show that these systems are compact and converged well.

Conformational analysis upon inhibitor binding

It has reported that the contact area in NSP10/NSP14 complex is ∼2275 Å2, involving 15–19, 40–45, 68–80 and 93–96 residues [36]. The interactions are mainly intermolecular hydrogen bonds and hydrophobic interactions. And the important residues involved in the NSP10/NSP16 complex (buried surface area of ∼933 Å2), including Cys41, Lys43, Leu45, Ala71, Ser72, Arg78, His80, Gly94 and Tyr96 [16]. The NSP10/NSP14 and NSP10/NSP16 interfaces have extensive overlapping areas, and in NSP10/NSP14, the N-terminal helix α1 is essential for NSP10 binding to NSP14 [37]. The conformation of NSP10 is altered after molecules binding to it, and the changes disrupt the binding of NSP14 or NSP16.

The results show that after the inhibitor binding, Phe16 of NSP10 rotated about 10° and Phe19 deviated from the equilibrium position. Compared with 4 and 5, greater deflection of Phe16 aromatic ring can be observed in 1 and 3. This affects the van der Waals interaction between NSP10 and NSP14 with helix H4, suggesting important hydrophobic interactions are lost at the interface. Ala20 and Val21 also have greater deflection after the inhibitor binding. In addition, larger flips of Asn40 and Cys41 are found in intervening loops for 1 and 3 compared to 4 and 5. Leu45 and Pro59 observed small offsets at four hits. The aromatic ring of Tyr76 exhibits significant rotation in hits, especially 3 and 5. Whereas His80 shows larger shift in hits, especially 1 and 3. Moreover, 93–96 also all have remarkable deflections in hits, where Lys93 exhibits towards the opposite direction after binding. Residues including Phe68, Gly70, Ala71 and Arg78 are slightly shifted in interface region. They are involved in polar interactions binding NSP14 and NSP16 and thus may interfere with the fit between the two proteins. The conformational changes of residues in NSP10 binding four hits are shown in Figs. S4–7.

As previously mentioned, the binding of these different molecules may affect ExoN and 2′-O-MTase activity by inducing interfacial residues. The results show that the residues participate in complex formation involving helices α3, η1, strands β2, β3 and intervening loops, displaying large conformational changes. Thus, binding of allosteric agents results in changes for residues at the distal interface and may block NSP14 or NSP16 binding at the interaction interface.

Conclusions

Inhibition of NSP10 has emerged as a potential therapy for SARS-CoV-2. The search for small molecule drugs targeting NSP10 may provide an effective therapeutic strategy for antiviral therapy. In the current study, virtual screening from natural product database was used to identify potential small molecule inhibitors targeting SARS-CoV-2 NSP10 in silico. The docking results show that the selected compounds form hydrogen bonds with potential key amino acid residues in the binding pocket, including Asp91, Tyr76 and Asp22. The MM-GBSA assessment reveals 11 compounds with binding energies in excess of −80.00 kcal/mol. The higher binding free energy indicates that these NSP10-ligand complexes have high binding affinity. Based on docking scores, binding free energies and ADMET properties, four representative compounds 1, 3, 4 and 5 were obtained. They form complexes with excellent binding free energies of −103.05, −95.82, −93.09 and −92.04 kcal/mol. In addition, all four compounds conformed to the Lipinski's rules and had good druglike properties. Toxicity prediction results show low mutagenic and carcinogenic toxicity, indicating positive safety profiles for the selected compounds. Furthermore, 100 ns MD simulations were performed to validate the binding modes and stability of the inhibitors with NSP10. MD analysis demonstrate that all candidate compounds are stable in binding to the protein system with fine molecular dynamics equilibrium. Taken altogether, the study of potential natural inhibitors of SARS-CoV-2 NSP10 as a therapeutic target. The identified candidate compounds targeting NSP10 have potential as inhibitors against SARS-CoV-2 and might be used as a starting point for the development of novel drugs.

Declaration of competing interest

The authors declare no conflict of interest.