Structure-based virtual screening and molecular dynamics simulation of SARS-CoV-2 Guanine-N7 methyltransferase (nsp14) for identifying antiviral inhibitors against COVID-19.

The recent pandemic caused by Severe Acute Respiratory Syndrome Coronavirus-2 (SARS-CoV-2) calls the whole world into a medical emergency. For tackling Coronavirus Disease 2019 (COVID-19), researchers from around the world are swiftly working on designing and identifying inhibitors against all possible viral key protein targets. One of the attractive drug targets is guanine-N7 methyltransferase which plays the main role in capping the 5'-ends of viral genomic RNA and sub genomic RNAs, to escape the host's innate immunity. We performed homology modeling and molecular dynamic (MD) simulation, in order to understand the molecular architecture of Guanosine-P3-Adenosine-5',5'-Triphosphate (G3A) binding with C-terminal N7-MTase domain of nsp14 from SARS-CoV-2. The residue Asn388 is highly conserved in present both in N7-MTase from SARS-CoV and SARS-CoV-2 and displays a unique function in G3A binding. For an in-depth understanding of these substrate specificities, we tried to screen and identify inhibitors from the Traditional Chinese Medicine (TCM) database. The combination of several computational approaches, including screening, MM/GBSA, MD simulations, and PCA calculations, provides the screened compounds that readily interact with the G3A binding site of homology modeled N7-MTase domain. Compounds from this screening will have strong potency towards inhibiting the substrate-binding and efficiently hinder the viral 5'-end RNA capping mechanism. We strongly believe the final compounds can become COVID-19 therapeutics, with huge international support.[Formula: see text]The focus of this study is to screen for antiviral inhibitors blocking guanine-N7 methyltransferase (N7-MTase), one of the key drug targets involved in the first methylation step of the SARS-CoV-2 RNA capping mechanism. Compounds binding the substrate-binding site can interfere with enzyme catalysis and impede 5'-end cap formation, which is crucial to mimic host RNA and evade host cellular immune responses. Therefore, our study proposes the top hit compounds from the Traditional Chinese Medicine (TCM) database using a combination of several computational approaches.Communicated by Ramaswamy H. Sarma.

Introduction

Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2), a recent pandemic outbreak, which was first reported in Wuhan (Hubei province), China, causing novel coronavirus pneumonia (COVID-19), resulted in more than 298,000 deaths, as on 14th May 2020. Even though COVID-19 presents relatively mild symptoms compared to the recent outbreaks like Ebola or Middle East Respiratory Syndrome (MERS), they are highly transmissible and can lead to severe pneumonia and death (Boopathi et al., 2020; DiMaio et al., 2020). The SARS-CoV-2 is an enveloped virus belongs to beta coronavirus family with a large and complex positive-sense, single-stranded RNA genome (∼30kb) similar to the other six known human coronaviruses - CoVs (SARS-CoV, MERS-CoV, HKU1, NL63, OC43, and 229E) (Fung & Liu, 2019). Genome analysis revealed SARS-CoV-2 is highly similar to bat CoV (SARSr-CoV-RaTG13), which has evolved and pre-adapted in intermediate animal species to bind human ACE2 receptor with increased affinity compared to other known CoVs (Andersen et al., 2020). SARS-CoV-2 also enters the host cell through receptor binding with viral spike protein and membrane fusion, thereby resulting in the release of genomic RNA, which further replicates and produces several subgenomic RNAs (Wang et al., 2020). The genomic and sub-generic RNAs are capped, thereby limiting the degradation by cellular 5′-3′ exonuclease (Ferron et al., 2012). Lai and Stohlman in 1981 discovered the RNA cap structure produced by a rat CoV (murine hepatitis virus) and later well studied in viruses from family Reoviridae and Poxviridae (Lai & Stohlman, 1981).

The SARS-CoV-2 genomic RNA is polycistronic with two large ORFs (1a and 1b) which encodes several non-structural proteins (nsp) forming the viral replicase-transcriptase enzyme complex (RTC) and several other ORFs which encode structural and accessory proteins, by cap-dependent translation and ribosomal frameshifts (Abdelli et al., 2020; Nakagawa et al., 2016). The viral RNAs can be targeted by the host immune system, therefore they have evolved a mimicking canonical capping mechanism and cap the 5′-ends of the RNA molecules similar to its eukaryotic host for escaping innate cellular immune response and efficiently replicate (Koonin & Moss, 2010). The essential enzymes for RNA methylation in SARS-CoV-2 are encoded by ORF1b, namely nsp14 which is bifunctional, includes exoribonuclease (ExoN) (Khan et al., 2020; Minskaia et al., 2006) in the N-terminus that plays an essential role in proofreading and C-terminus Guanine-N7 methyltransferase (N7-MTase) (Chen et al., 2009) along with nsp16 or (Nucleoside-2′-O) methyltransferase (2′O-MTase) which function with its co-factor nsp10 activator protein are essential for methylation of 5′-end RNA cap. In coronaviruses, the initial cap core structure (cap-0) is formed at the 5′-end of RNA by two capping enzymes RNA 5′-triphosphatase (RTPase/nsp13) hydrolyzes nascent RNA to pp-RNA, which is followed by GTP hydrolysis and transfer the product GMP to pp-RNA forming Gppp-RNA catalyzed by RNA guanylyltransferase (GTase). Later this minimal cap structure is methylated at the N7-position of the capped Guanine by N7-MTase and then finally 2′-O position of first nucleotide’s ribose by 2′O-MTase resulting in a cap-1 and methylation of the subsequent nucleotide form cap-2 structure. During the methylation reaction mechanism, the methyl group is transferred from donor S-adenosyl-L-methionine (SAM/AdoMet) and added to the RNA substrate by forming S-adenosyl-L-homocysteine (SAH/AdoHcy) as a byproduct (Chen & Guo, 2016; Decroly et al., 2011). Therefore, a 5’cap consist of a 7-methylguanosine (m7G) linked to a 5′-5′ triphosphate bridge followed by a methylated nucleotide 2’hydroxyl group, represented as m7GpppNm2′O, where N represents the first nucleotide (Figure 1).

Figure 1.

Schematic pathway representation of coronaviral (CoV) RNA capping mechanism of showing the target step for proposed antiviral inhibitors against the SARS-CoV-2 N7-MTase resulting in viral RNA degradation by the host immune response and chemical structure of a viral RNA cap-2 formed at the 5’-end of genomic and sub-genomic RNAs.

Some RNA viruses (e.g. flaviviruses) has evolved single enzyme to catalyze both methylation steps (Hercik et al., 2017). Methylated 5′-end cap structures confer stability and allow the ribosome to recognize the viral mRNAs by eukaryotic translation initiation factor 4E (eIF4E) for translation and its crucial in viral host genome replication within a host cell cytoplasm (Ma et al., 2015). The N7-MTase domain structure exhibits a noncanonical MTase fold with a rare β -sheet insertion and a peripheral zinc finger. Even though the overall structure of the viral N7-MTase domain has high similarity to its cellular counterpart, there are several structural differences that should be considered to develop a viral target specific antiviral drug (Elmezayen et al., 2020; Ferron et al., 2012; Ma et al., 2015). In this study, we have modeled SARS-CoV-2 N7-MTase protein structure in order to decipher the molecular architecture of the substrate-binding site and to perform virtual screening and docking analysis to identify an attractive and novel lead compound from TCM@Taiwan database, which can assist researchers in designing a potent antiviral drug against COVID-19 by blocking the SARS-CoV-2’s RNA capping mechanism required for evading the host defense.

Materials and methods

Homology modeling

The experimental structure of SARS-CoV-2’s Guanine-N7 methyltransferase enzyme is unavailable in the Protein Data Bank (PDB), therefore the amino acid sequence of nsp14 is retrieved from the NCBI Reference Sequence: YP_009725309.1 (https://www.ncbi.nlm.nih.gov/) for the purpose of protein homology modeling. Suitable template structures are searched within Protein Data Bank using BLASTp with default parameters (Rayalu et al., 2012; Selvaraj, Sivakamavalli, Vaseeharan, et al., 2014). Template deciding factors are structure with maximum identity with a high score, lower e-value, and functionally similar proteins (Fazil et al., 2012; Sivakamavalli et al., 2014). The sequence alignment between the template with the target protein is performed by using the PRALINE tool for better visualization and to compare/analyze the high and low conserved amino acids between template and target sequences (Bawono & Heringa, 2014; Selvaraj, Sivakamavalli, Baskaralingam, et al., 2014). For 3D protein structure generation, the academic version of MODELLER 9.2 (https://salilab.org/modeller/) is used along with the obtained sequence alignment (Webb & Sali, 2016, 2017). Ten models are generated and from that low objective function model is subject for the stereochemical check for standardizing the normal bond length, dihedrals, and non-bonded atom-atom distances (Gupta et al., 2020; Sarma et al., 2020). The model protein is validated for initial protein quality checks using SAVES server(Pontius et al., 1996; Selvaraj, Sivakamavalli, Vaseeharan, et al., 2014).

Protein and ligand preparation

The final model structure is proceeded to the preparation for direct usage in the molecular modeling environment, due to the misled bond orders, topologies, or formal atomic charges (Kumar et al., 2018; Lionta et al., 2014; Vijayalakshmi et al., 2013). For obtaining the correct structure, the protein preparation wizard in Schrödinger is used to apply for correct bond orders, missing atoms, and side-chain refinement along with correctness of partial atomic charges (Selvaraj, Bharathi Priya, et al., 2014; Tripathi et al., 2012). Amino acid flipping is incorporated for side-chains angles of amino acids (Asn, Gln, and His) and it may influence the hydrogen bond formation and generate tautomer/ionized states (Sastry et al., 2013). Finally, the whole structure is optimized with the intramolecular hydrogen bonds and by applying the OPLS-3e FF, the whole structure is minimized until the RMSD (root-mean-square deviation) of all atoms reaches 0.30 Å (Selvaraj & Singh, 2014). For ligand preparation, the substrate of the Guanine-N7 methyltransferase namely G3A (Guanosine-P3-Adenosine-5′,5′-Triphosphate) is downloaded from the PubChem (ID: 135450590) and for screening approaches, 22122 molecules from TCM Database@Taiwan is used (Ma et al., 2016; Wang et al., 2018). The pH values are maintained with 7.0 ± 2.0 for attaining the ligand ionization state and stereoisomers are generated for each ligand structure (Aanouz et al., 2020; Onufriev & Alexov, 2013). Both the substrate and ligand molecules are subjected to the LigPrep minimization OPLS3e force field with the output up to ten conformations per ligand will be generated (Harder et al., 2016; Selvaraj, Singh, et al., 2014b).

Substrate docking and substrate based grid generation

The substrate molecule G3A is matched with template protein and the binding substrate region is marked. The G3A is allowed to bind in the same region the modeled protein and prepared with similar protocol discussed in the protein preparation (Das et al., 2020; Zhou et al., 2009). The protein complexed with G3A region is taken for substrate-based grid generation by manual picking of the substrate and customize grid box around it. By this, the focused region of docking is set to 2 Å around the substrate-binding region. The process of substrate matching is performed using the Maestro (Schrödinger) and the grid generation using the Glide (Alogheli et al., 2017; Pant et al., 2020).

High throughput virtual screening (HTVS)

In order to identify the appropriate compounds that perfectly fits the substrate-binding site of the G3A, the Glide incorporated virtual screening protocol is used (Grover et al., 2012; Kawatkar et al., 2009; Muralidharan et al., 2020). The prepared TCM compounds are considered for virtual screening. ADME screening is not considered because the compounds have a long history of being used for medicinal purpose (Gyebi et al., 2020; Wang et al., 2019). The whole database is initially screened with High Throughput Virtual Screening (HTVS), for identifying appropriate compounds that fit in the substrate-binding site (Lobo-Galo et al., 2020; Patidar et al., 2016). The succeeded compounds in the HTVS will be many in number and from that top 10% be forwarded to SP (Standard Precision) docking and top 10% compounds will from SP docking will again forward to XP (Extra Precision) docking (Aksoydan et al., 2018; Verma et al., 2018). The top twenty compounds final pose is ranked based on Prime MM/GBSA (molecular mechanics/Poisson–Boltzmann (Generalized Born) surface area) ranking and Glide score and from that top five compounds are subjected to MD simulations (Rapp et al., 2011; Selvaraj, Singh, et al., 2014a).

Molecular dynamics simulation

The molecular dynamics (MD) is performed for the G3A complex nsp14 containing Guanine-N7 methyltransferase domain and for the top five ligand-bound complexes through the virtual screening using the GROningen MAchine for Chemical Simulations (GROMACS 5.1.4, http://gromacs.org) (Chavez Thielemann et al., 2019; Rakhshani et al., 2019). Both the substrate-bound complex and ligand-bound complexes are simulated for the 50 ns of timescale for understanding the dynamic behavior and to understand the interaction pattern (Selvaraj et al., 2018; Shafreen et al., 2013; Swegat et al., 2003). The top five ligand-bound structures are prepared using the OPLS-AA force field solvated by the TIP3P water model within a periodic boundary box of distance 1.0 nm, fixed in between the protein and cubic box (Bandaru et al., 2017; Nguyen et al., 2014; Umesh et al., 2020). The charge topology of ligand molecules is externally added through the PRODRG server (http://prodrg1.dyndns.org/) (Schuttelkopf & van Aalten, 2004; van Aalten et al., 1996). For neutralizing the whole system, the accurate concentration of (Na+/Cl−) ions is added based on the rebalancing charges. Initial energy minimization is appending with the prepared complex systems for 1,000 steps of steepest descent algorithm via a tolerance of 10 kJ/mol/nm to avoid the steric clashes (Kumar et al., 2020; Sliwoski et al., 2014). Thermostat coupling is set with a reference temperature of 300 K using Berendsen thermostat and pressure coupling with 1.0 bar reference pressure using Parrinello-Rahman along with periodic boundary conditions with cut-offs for Lennard-Jones and Coulomb interactions. Particle-Mesh Ewald method is used for calculating the long-range interactions for biomolecular systems (Chinnasamy et al., 2019; Selvaraj et al., 2015; Wohlert & Edholm, 2004). After the initial minimization step, the whole systems are again well equilibrated 1,000 ps at 300 K and 1 bar pressure in NVT and NPT ensembles (Childers & Daggett, 2018; Chinnasamy et al., 2020). Final MD simulation step is processed for all the protein-ligand complex with respect to the timescale of 50 ns. For MD simulation analysis, the results of RMSD, RMSF, bonding interactions, and PCA are performed using the VMD and GROMACS tools (Martinez, 2015).

Results and discussion

Model searching and sequence alignment

The experimental structure of the Guanine-N7 methyltransferase enzyme is unavailable in the PDB and BLASTp results suggest the crystal structure of the SARS-CoV nsp14-nsp10 complex (PDB ID: 5c8s) as the template for structure modeling. The model template structure used is B-chain from the protein complex by showing a 95.06% identity and 98% similarity. The template structure is having higher closeness with the target sequence in both phylogeny and functional behavior. Both the template and target are functionally holding the role of forming the cap-1 structure at the 5′-end, which assists in translation and evading host defense. Thus, the template, Guanine-N7 methyltransferase from SARS-CoV is taken for the construction of the homology model. Sequence information of both the proteins is taken for the sequence alignment with PRALINE tool for understanding the conserved and unconserved regions. As the template structure shows a 95.06% identity and 98% similarity in the BLAST, the PRALINE tool predicts the high conserved regions between both proteins as shown in Figure 2. The query target sequence with the showing the dominance of red color indicating highly conserved between both sequences, and there is 13 unconserved amino acid in the ExoN domain and 10 unconserved amino acids in the N7-MTase domain, which showing the differences. We have also noted that the template protein is having the missing residues from the 454th to 464th position.

Figure 2.

Target-template sequence alignment of bifunctional nsp14 containing N-terminal exoribonuclease (ExoN) and C-terminal guanine-N7 methyltransferase domains (shown side-by-side), and lower panel show the domain organization and boundaries of nsp14 protein.

Homology modeling and structure validation

From the outputs of sequence alignment, the 3D structure model of nsp14 protein is created using the homology modeling method, as the similarity between the template and target sequence is extremely high. MODELLER predicted ten models and from that, the low Z-score protein conformation is considered. We also noticed, the gap in the template structure is filled with loop regions in the final modeled protein. The detailed secondary structure analysis is performed with PDBsum, which shows that the template structure from SARS-CoV is having the topology arrangement of 6 beta-sheets, 8 beta-hairpins, 3 beta-bulges, 23 beta-strands, 13 helices, 10 helix-helix interacs, 60 beta-turns, and the modeled protein is having the topology arrangement of 6 beta-sheets, 8 beta-hairpins, 2 beta-bulges, 24 beta-strands, 14 helices, 6 helix-helix interacs, 63 beta-turns, 2 gamma-turns as shown in Figure S1(a, b). The secondary structure of both template and target structure shows similar topology, and thus the 3D structure of proteins is viewed and provided in Figure 3. Figure 3(a) shows the overall structure of nsp14 with the C-terminal Guanine-N7 methyltransferase domain structure of SARS-CoV used as a template for modeling and Figure 3(b) shows the modeled structure of SARS-CoV-2’s nsp14 structure with C-terminal Guanine-N7 methyltransferase domain. Both the structures look similar to each other and we noticed some structural changes between them. From Figure 3(a, b), we noticed four regions. The first region in the template showing an alpha helix, and in the same region, the model protein shows the loop region and we observed the order to disorder transition in the SARS-CoV-2.

Figure 3.

Overall structures of nsp14 protein displaying N-terminal exoribonuclease (ExoN) and C-terminal guanine-N7 methyltransferase (N7-MTase) domains: (a) template nsp14 structure of SARS-CoV (PDB ID: 5c8s) used for modeling. (b) modeled target nsp14 structure from SARS-CoV-2. (c) Structural superimposition of nsp14 from SARS-CoV with the homology modeled nsp14 from SARS-CoV-2.

The second region in the template shows the alpha helix, and for the model protein, the same region shows the loops, and again the order to disorder transition is noted. Similarly, the third and fourth regions are also showing the order to the disorder transitions. For a clear view of similar protein topology, the template (blue) and modeled protein (magenta) are structurally aligned and shown in Figure 3(c). The RMSD values between the template and the model structure show 0.377 Å and this RMSD values indicate that there is a high structural similarity between the proteins. The geometric structural validation of the modeled protein is shown through the Ramachandran 2D plot, which shows the statistical distribution of the combinations φ-ψ dihedral angles of the protein backbone. This can be used for the structure validation based on the distribution of Phi/Psi values plotted from the protein structure. The actual role of the Ramachandran plot is to visualize the energetically allowed and disallowed regions of the dihedral angles. If the modeled protein shows a higher allocation percentage in disallowed regions, the protein model stands back as the poor quality. This is due to the placement of backbone dihedral angles ϕ and ψ in the wrong orientation or statistical distribution indicate the problems with the structure. The Ramachandran plot for the model structure of SARS-CoV-2’s nsp14 with the C-terminal Guanine-N7 methyltransferase domain shows the good quality of the model is generated through homology modeling. Figure S2 Ramachandran plot shows the absence of disallowed regions in the modeled protein and the statistical distribution of the combinations φ-ψ dihedral angles of the protein backbone is perfect.

Substrate binding site of N7-MTase from SARS-CoV vs SARS-CoV-2

The G3A is the main substrate for the activation of the Guanine-N7 methyltransferase enzyme. The reaction mechanism encoded by G3A of the N7-MTase enzyme is completely relied on this substrate binding for possessing the cap structure at 5′-ends of viral genomic RNA and sub-genomic RNAs. Thus, it is mandatory to study the binding interactions of G3A with the N7-MTase for both SARS-CoV and SARS-CoV-2. As the template with a high sequence similarity-based model is performed, we expected a similar binding mode of G3A substrate binding in N7-MTase for SARS-CoV-2. The binding interaction of G3A with N7-MTase for SARS-CoV shows the interactions with Cys309 and Asn388 as shown in Figure 4(a). But the binding interactions of G3A with N7-MTase for SARS-CoV-2 show the interactions with Asn306, Arg310, Trp385, and Asn388 as shown in Figure 4(b). In between both the SARS-CoV and SARS-CoV-2 N7-MTase enzymes, the G3A substrate binds with the common Asn388 residue. The Cys309 in the SARS-CoV N7-MTase enzyme holds the charge of polar hydrophobic neutral and it has been replaced with Asn306, Arg310, Trp385 residues in SARS-CoV-2 N7-MTase. Here, Asn306 is the polar hydrophilic neutral amino acid, Arg310 is the polar hydrophilic positively charged amino acid, Trp385 is aromatic hydrophobic neutral amino acid. Based on analyzing these interactions, we found Asn388 is core important residue for the substrate binding and, the SARS-CoV-2 N7-MTase may have the additional support of Asn306, Arg310, and Trp385 residues in the reaction mechanism.

Figure 4.

Substrate-binding site of the N7-MTase domain of nsp14 highlighting the bound Guanosine-P3-Adenosine-5’,5’-Triphosphate (G3A) (a) template (PDB ID: 5c8s) from SARS-CoV and (b) modeled structure from SARS-CoV-2, the interacting amino acid residues are shown as sticks.

Virtual screening of the traditional Chinese medicine (TCM) database

To disturb the G3A binding towards the Asn388, the whole TCM database is searched for a suitable compound using the HTVS method. The funnel-based screening method yields the best compounds that suit the substrate-binding pocket of SARS-CoV-2 N7-MTase. The XP docking and Prime MM/GBSA provides the TCM 57025 (5a), TCM 3495 (5b), TCM 20111 (5c), TCM 31007 (5d) and TCM 5376 (5e) from the TCM database as the best compounds that can bind with the substrate-binding site of SARS-CoV-2 N7-MTase as shown in Figure 5(a–e). The scoring values of Glide XP docking and prime MM/GBSA are provided in Table 1, show all the screened compounds with the docking score > −8.5 kcal/mol and binding energy values > −45 kcal/mol.

Figure 5.

Molecular interaction of the top hits of TCM compounds resulted from virtual screening: TCM ligands interacting with the substrate-binding site (a) TCM 57025, (b) TCM 3495, (c) TCM 20111, (d) TCM 31007, and (e) TCM 5376 showing labelled interacting residues obtained from XP docking.

Table 1.

TCM database screened hit compounds scoring in extra precision (XP) docking and prime MM/GBSA.

Compound name	Docking score	Glide emodel	Glide ecoul	Glide evdw	Glide energy	ΔGbind
TCM 57025	−11.486	−78.55	−26.461	−33.703	−60.164	−55.640
TCM 3495	−9.889	−68.837	−13.451	−42.042	−55.492	−46.982
TCM 31007	−8.867	−85.896	−14.34	−48.268	−62.608	−48.879
TCM 20111	−8.775	−99.491	−20.796	−51.131	−71.927	−51.621
TCM 5376	−8.529	−108.033	−15.241	−53.811	−69.052	−43.296

unit - (kcal/mol).

This shows that the compounds are exceptional in binding towards the substrate-binding pocket and can readily inhibit the reaction mechanism that possess the cap structure at 5′-ends of viral genomic RNA and sub genomic RNAs. The insights of interacting amino acids are analyzed and tabulated in Table 2. The table values show that except the compound TCM 5376, all the other screened compounds are having direct hydrogen-bonding interactions with Arg310 amino acid. The compound TCM 57025 and TCM 5376 are having interactions with Asn388, which plays a major role in holding the G3A substrate. Apart from these amino acids, Asn386 is showing a high interacting pattern with new screened compounds and we believe that Asn386 may hold the strong lead in the inhibition of SARS-CoV-2 N7-MTase. For the interactions with the compound TCM 57025, the aromatic ring of Phe426 showing the Π-Π interactions with the 1,2-dihydroxybenzene. The imidazole functional group present in the His424 also shows the Π-Π interaction with another 1,2-dihydroxybenzene. In the same 1,2-dihydroxybenzene, one Π-cationic interaction is formed with Lys423. Overall, the compound TCM 57025 shows six hydrogen bonds, two Π-Π interaction, and one Π-cationic interaction along with the strong scoring of −11.486 kcal/mol docking score and −55.64 kcal/mol binding energy. The aromatic ring of Phe426 is also showing the Π-Π interaction with the compounds namely TCM 31007 and TCM 20111.

Table 2.

Screened hit compounds interaction with substrate-binding site of SARS-CoV-2 N7-MTase.

Compound name	No. of H-bond	Interacting amino acids	Π-Π interaction
TCM 57025	06	Glu295, Arg310, Lys336, Asn386, Asn388, Leu526	Yes
TCM 3495	06	Cys309, Arg310, Gln313, Leu526	No
TCM 31007	06	Trp292, Arg310, Gln313, Lys336, Asn386, Leu526	Yes
TCM 20111	06	Arg310, Lys336, Lys423, His424, Leu526	Yes
TCM 5376	04	Asn306, Asn386, Asn388	No

Molecular dynamics (MD) simulations

The MD simulations are processed with substrate and ligand bounded SARS-CoV-2 N7-MTase for the timescale of 50 ns for understanding the stability and dynamic behavior in an aqueous environment.

For both substrate and the lead molecules, the MD simulation is performed in a similar protocol and the values of RMSD are plotted in Figure 6. The SARS-CoV-2 N7-MTase complexed with substrate G3A shows a stable conformation throughout the simulations, which is shown in black color. This substrate complex is dynamically stable and found to place in the position between ∼−0.45 nm to ∼−0.55 nm and throughout the simulation, we did not see any sudden surge or sliding from this position. This clearly says that, the substrate-bound complex is dynamically stable due to the tight amino acid contributions from Asn306, Arg310, Trp385, and Asn388. Along with the substrate binding, the new compounds from the screening also analyzed for the RMSD values and compared in Figure 6. The RMSD plot shows that all the new compounds are stable throughout the simulations and RMSD values did not show any sudden surge or sliding throughout the simulations. From these compounds, notably the compound TCM 3495 and TCM 20111 shows the RMSD values ∼>0.94 nm, which is higher than the substrate complex and other ligand complexes. The average mean value of substrate binding is ∼0.54 nm, while the average RMSD of the compounds TCM 57025, TCM 3495, TCM 5376, TCM 20111 and TCM 31007 is ∼0.68 nm, ∼0.88 nm, ∼0.58 nm, ∼0.76 nm, and ∼0.56 nm respectively. Analyzing the residue wise fluctuations through RMSF may suggest the important core residues for the RMSD value fluctuations. Thus, the RMSF values for the substrate complex and ligand complexes are plotted in Supplementary Figure S3. The RMSF plot shows notable deviations have occurred in the N-terminal and C-terminal regions for all the complexes, but the other residues are not fluctuating. This indicates the reason for stable RMSD plots for all the substrate and ligand-bound complexes. Understanding the binding phenomenon of substrate and ligand complex is analyzed for its hydrogen-bonding interactions. The hydrogen-bonding interactions are plotted in Figure 7, which shows the substrate and ligand complex are showing the strong bonding interactions throughout the 50 ns of the MD simulations. The substrate G3A bound complex show 0–3 hydrogen bonds up to 3 ns, and after the 3 ns the hydrogen bond stabilizes the complex and showing from ∼4–6 hydrogen bonds. The substrate molecule bound with SARS-CoV-2 N7-MTase showing the average hydrogen bond of ∼3.6 hydrogen bonds throughout the MD simulations. But for the average hydrogen bond tendency of the compounds TCM 57025, TCM 3495, TCM 5376, TCM 20111, and TCM 31007 is ∼3.40, 3.35, 4.00, 3.97 and 3.34, respectively. On analyzing the binding of substrate and ligand molecules, the SARS-CoV-2 N7-MTase enzyme conformations are also necessary to be analyzed. For that, the principal component analysis (PCA) is performed from the trajectories by calculating the atomic motions. The determination values are calculated from the overall motion of the protein-substrate and protein-ligand complex by its eigenvector and eigenvalues by applying principal components and covariance matrix for examining the necessary phase space. Trajectories are subjected to spread over on the eigenvectors PC1 and PC2 (well-defined clusters) as represented in Figure S4(a–f). The PCA results indicated that the substrate complex is compact in space and resembling that, except the TCM 20111 compound, all the other ligand complex shows compactness, which implicates higher stability seen over the protein-ligand interactions. Overall, the substrate and ligand compounds show potential binding over the 50 ns MD simulations by showing narrow energy potential and stable conformations.

Figure 6.

The RMSD graph for the entire timescale (50 ns) of MD simulation shown for substrate (G3A) and the best TCM compound complex of nsp14 from SARS-CoV-2 for exploring the conformational landscapes, lower panel display the zoomed view of the highlighted region.

Figure 7.

Overall hydrogen-bonding interaction plot of ligand-bound complexes with SARS-CoV-2 substrate-binding site of N7-MTase and its substrate G3A.

Conclusion

SARS-CoV-2 has evolved with the ability of capping their RNAs through similar canonical RNA capping pathway found in higher eukaryotes, for successful replication inside the host cell. The N7-MTase enzyme is found in the C-terminal part of nsp14 protein, which also includes the exoribonuclease domain in its N-terminus. This RNA capping mechanism is playing a vital role in the viral RNA escaping from the immune cells and that failure of RNA capping leads to viral RNA degradation, eventually hindering the replication cycle. Thus, the N7-MTase is an attractive drug target to inhibit the SARS-CoV-2 thereby curing COVID-19. In this present work, we have solved the 3D model structure of complete nsp14 from SARS-CoV-2, which includes both the guanine-N7 methyltransferase and the exoribonuclease domain. Substrate binding interactions are analyzed for N7-MTase and as predicted Asn306, Arg310, Trp385, and Asn388 are important residues to inhibit the functional reaction mechanism initiated by N7-MTase. Based on the molecular docking and simulation studies, we found that TCM 57025, TCM 3495, TCM 5376, TCM 20111, and TCM 31007 are the compounds from the TCM database, which can occupy and interact nicely to the substrate-binding site of N7-MTase. These new compounds are potent and found to interact in the common binding pattern, especially in showing the interactions with the Arg310 and Asn388. These compounds on further in vitro and in vivo experiments can become strong candidates for clinical trials. Overall, we believe that these five compounds from the TCM database have a strong potential as the anti-viral phytochemicals that may inhibit SARS-CoV-2 N7-MTase. In addition, the structure and mechanistic insights provided in this study may provide a deep understanding for developing anti-SARS-CoV-2 inhibitors.