Computational investigation of potent inhibitors against SARS-CoV-2 2′-O-methyltransferase (nsp16): Structure-based pharmacophore modeling, molecular docking, molecular dynamics simulations and binding free energy calculations

The Coronavirus Disease 2019 (COVID-19), caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) infection, has created unprecedented public health and economic crises around the world. SARS-CoV-2 2′-O-methyltransferase (nsp16) adds a “cap” to viral RNA to maintain the stability of viral RNA, and inhibition of nsp16 activity may reduce viral proliferation, making this protein an attractive drug target. Here, we report the identification of several small molecule inhibitors of nsp16 by virtual screening. First, the nsp16-sinefungin complex (PDB ID: 6WKQ) was selected from the protein data bank. Asp6912, Cys6913, Asp6897 and Asp6928 were determined to be the key amino acids for sinefungin binding in the crystal structure of nsp16-sinefungin complex by molecular dynamics simulation. The complex structures in the stable binding trajectory of nsp16-sinefungin were than clustered through molecular dynamics RMSD analysis. Six clusters were generated, and six representative structures were selected to construct the pharmacophore based on the structure. These six pharmacophores were superimposed on the binding pocket to simplify and pick the common characteristics. The compounds obtained by the pharmacophore screening from Bionet and Chembiv databases were docked into the nsp16 active pocket. The candidate compounds were selected according to the molecular docking score and then screened by MM/GBSA. Finally, four candidate compounds were obtained. Four sets of 150ns molecular dynamics simulations were performed to determine whether candidate compounds could maintain stable interactions with key amino acids. The results of MD and MM/PBSA energy decomposition indicated that C1 and C2 could form a stable complex system with nsp16, and could form strong hydrogen bonds and salt bridges with the key amino acid Asp6897 and Asp6928. This study thus identifies and attempts to validate for the first time the potential inhibitory activities of C1 and C2 against nsp16, allowing the development of potent anti-COVID-19 drugs and unique treatment strategies.


Introduction
Since December 2019, disease caused by severe acute respiratory syndrome (SARS-CoV-2) virus has been declared as an accessible disease by the World Health Organization (WHO) [1,2]. As of April 26, 2021, the virus has had an extremely dangerous impact on more than 200 countries, with more than 140 million confirmed cases and more than 3 million deaths [3]. However, few drugs are available for the treatment of COVID-19 infection, and the maximum efficacy of these drugs may only reduce the symptoms of infection [4]. Therefore, there is an urgent need to develop effective antiviral therapies to combat the pandemic.
The 5 ′ end of the coronavirus RNA is modified by a cap structure. Modification of the 5 ′ end of RNA virus genome is very important for the replication of the virus genome [5], the translation and synthesis of virus protein, the assembly of virus particles and the elusion of the host innate immune system from the recognition of virus pathogen molecular pattern [6,7]. Therefore, 5 ′ terminal modification of the RNA virus genome serves as a potential target for antiviral drug design. Because the replication life cycles of coronavirus takes place in the cytoplasm, the capping enzyme of the host cannot complete the capping of coronavirus itself [8,9], so coronavirus completes the capping of viral RNA cap1 through its own protein [10][11][12]. Additionally ribose 2 ′ -O-methylation plays an important role in the evasion of innate immune recognition. In SARS-CoV-2 virus, 2 ′ -O-methyltransferase (2 ′ -O-Mtase) consists of two subunits: catalytic subunit nsp16 and stimulating subunit nsp16 [13,14]. SARS CoV-2 2 ′ -O-methyltransferase (nsp16) is an important enzyme for virus survival [15][16][17]. The role of nsp16 is to protect viral RNA from the cellular innate immunity by participating in the formation of a special arrangement at the 5 ′ end of RNA molecule. Therefore, nsp16 inhibitor can effectively inhibit the infection of COVID-19 [18,19].
The use of computer-aided drug design (CADD) experiments can accelerate the process of drug development. Currently, a multitude of simulation experiments have been conducted to develop inhibitors for nsp16. Mahmoud et al. [20] screened the zinc database by constructing the nsp16-sinefungin complex pharmacophore model, and evaluated the candidate compounds by molecular dynamics and binding free energy in order to indentify potential nsp16 inhibitors. Aldahham et al. [21] screened naphthyridine and quinoline derivatives as potential nsp16-nsp10 complex inhibitors by similar methods. In addition to these, other screening efforts were directed towards natural product development studies. Anuradha et al. [22] screened 10 natural products for nsp16 inhibitors by molecular docking and concluded that theaflavins and catechins have the potential to inhibit the interaction of viral molecules. In another study, Fatoki et al. [23] used gene network analysis, molecular docking, and sequence and structural dynamics simulations to determine potential medicinal phytochemicals to modulate gene expression networks associated with the pathology of SARS-CoV-2 in the human host. Malik et al. [24] identified the plant compound Withanolide and the anti-HIV drug Dolutegravir as Nsp16 inhibitors after analyzing 128 plant compounds and 10 FDA-approved anti-HIV drugs through AutoDock Vina docking and Gromacs dynamics simulations analysis.
In this study, based on the reported nsp16 and sinefungin crystal complex structure (PDB ID: 6WKQ) [25], we applied a virtual screening workflow based on a combination of structure-based pharmacophore modeling with molecular dynamics simulation techniques in order to discover potent SARS-CoV-2 inhibitors.

Protein preparation
The crystal structure of the nsp16 in complex with the pan-methyltransferase inhibitor sinefungin was retrieved from the protein data bank (PDB ID: 6WKQ, resolution:1.98 Å). The crystal structure of nsp 16 in complex with sinefungin was prepared by Schrödinger 2018. First, water molecules and other cofactors were removed. Hydrogen atoms and amino acid residues of the crystal structure were subsequently added via "Protein Preparation Wizard" in the Schrödinger [26]. At pH 7.4, the protonated and tautomeric of amino acids were adjusted to the corresponding state. Finally, through the atomic force field OPLS 2005, the hydrogen atom energy of the crystal structure was minimized and the heavy atom converges to an RMSD of 0.3 Å [27].

Molecular dynamics simulations and clustering
Using the "System Builder" module in the Schrödinger, the solvent system was constructed for the crystal structure of nsp 16 in complex with sinefungin. The system was solvated in an orthorhombic box, with SPC as the solvent system, and then two Cl − were added for neutralization. The 0.15 M NaCl acts to reproduce physiological conditions. At the same time, the system was given OPLS_2005 force field. At 310K and normal pressure (NPT) for 50ns simulation, the prepared system was simulated by "Desmond 3.7" [28]. The energy and trajectory atomic coordinate data were recorded at 1.2ps and 100ps, respectively, and the data were used for statistical analysis. The most stable path in simulation was clustered. According to RMSD, "Desmond trajectory clustering" was used to cluster frames that were collected, and the frequency value was set to 1, resulting in a maximum of 6 clusters. The MD simulation for the complexes generated by docking studies were executed employing the aforementioned protocol.

Pharmacophore construction and virtual screening
The "Phase" in the Schrödinger was used for pharmacophore construction and virtual screening [29][30][31]. The representative conformations from 6 clusters were selected to construct the pharmacophore model, and the exclusion volume was adjusted according to the conformation of the binding pocket. Combined with molecular dynamics analysis, several representative pharmacophore models were merged and simplified to obtain the screening pharmacophore model.
The molecules in Bionet [35] and Chembiv (https://www.chemdiv. com/) databases were matched with the screening pharmacophores. In the screening process, up to 50 conformations were generated for molecules and the energy minimization was utilized when the small molecules were exported. The compound matched at least three of the characteristics in the pharmacophore. The screening results were ranked by phase screen score. The pharmacophore-based virtual screening was performed by "Phase Ligand Screening" in the Schrödinger.

Ligand Preparation
In the "LigPrep", the molecules matched with the screening pharmacophore were optimized using OPLS3 force filed. All molecules retain specified chirality, generate tautomers, and result in the corresponding ionized state at pH 7.0 ± 2.0 by using Epik (Empirical pKa Prediction) [32]. The co-crystalline compound performs the same treatment as described above.

Molecular docking
The docking simulation was performed in Glide. Re-docking was implemented to verify the accuracy of Glide docking protocol. The two docking precisions "standard precision (SP)" and "extra precision (XP)" were used in the docking simulation [33]. Therefore, both docking methods were used in the re-docking of the exfoliated co-crystalline compound. The root mean square deviation (RMSD) between crystal and re-dock pose was calculated. The RMSD value less than 2 Å can confirm that the conformation produced by the docking can reproduce the X-ray results accurately [34]. The receptor grid file was defined according to the sinefungin position, and the parameter coordinates are shown in Table 1. Docking score and glide score were used to evaluate the screening results.

MM/GBSA calculation
The stability of the ligand-protein complex system generated by the docking simulation was verified by the calculation of the free binding energy in Prime MM/GBSA. The free binding energy calculation abide by the following formulas: ΔG bind is the binding energy; ΔG complex is the free energy of the complex system; ΔG protein and ΔG ligand are the free energy of protein and ligand in the complex system; ΔG GB and ΔG SA are the contribution of polarity and nonpolarity to solvent free energy in solvent environment. ΔS represents the entropy change of ligand structure during sampling. ΔE MM is the Gas phase free energy. The parameters in MM/GBSA are the default settings: solvation model is VSGB, and sampling method is Minimize. Check "use constraints on flexible residues" to set the residues around the receptor pocket as flexible conformation during calculation.

MM/PBSA and energy decomposition calculation
Protein-ligand complexes were then subjected to a molecular dynamics study and MM-PBSA binding free energy calculations by Gro-macs2021. The RESP charges were calculated for ligands by Gaussian16, and the GAFF force field was subsequently assigned by the Bio2Byte online server (https://www.bio2byte.be/). The AMBER99SB force field was assigned to the protein in Gromacs. The position of the ligand in the complex system was restricted. The complex system was then placed in a dodecahedral volume (the edge of the volume box was 10 Å from the molecular edge of the complex system), and the TIP3P solvent model was added to the volume box. The "genion" module in Gromacs adds counter ions (Na + or Cl-) for the system to maintain electrical neutrality. Before the molecular dynamics simulations simulation, the energy of the system was minimized, 100ps equilibration was then carried out under the NVT and NPT systems respectively. The equilibrated system was finally subjected to 100ns dynamics simulation, and the trajectories were sampled every 100ps. The "trjconv" module in Gromacs removes the coordinates for the trajectory and corrects the periodicity. The binding free energy of the ligand-protein complex system and the energy contribution of the amino acids in the binding pocket were analyzed from the four terms ΔG van der Waals , ΔG Electrostatic , ΔG Polar Solvation and ΔG Non-Pol Solvation by the MMPBSA.py script in the Amber Tool 21 software.

Molecular dynamics simulation
We analyzed the interaction between sinefungin and nsp16. Crystal structure 6WKQ interactions between sinefungin and nsp16 in PDB database revealed that sinefungin forms multiple hydrogen bonds with amino acid. As shown in Fig. 1a, the carboxylic acid on sinefungin forms hydrogen bonds with Asn6841 and Gly6879, Asp6897 as the hydrogen acceptor forms hydrogen bonds with two hydroxyl groups. The amino group and the nitrogen on the purine act as the donor and acceptor to form hydrogen bonds with Asp6928, Gly6869, Leu6898, Asp6912 and Cys6913. In addition, Asp6928 forms a salt bridge with amino cations.
In order to verify the stability of the interaction between sinefungin and nsp16, a 50ns molecular dynamics simulation was performed for the crystal structure of the nsp16-sinefungin complex, and the stability of the ligand in the binding pocket and the stability of the protein were evaluated by RMSD during the simulation ( Fig. S1 in Supplementary Materials). The results showed that sinefungin maintained a relatively stable state in the binding pocket after 30 ns. In order to explore the main interaction between the ligand and receptor in the simulation process, the 50% of ligand-amino acid interaction was analyzed. As shown in Fig. 1b, the main interaction is dominated by hydrogen bonds. The proportions of hydrogen bonds between sinefungin and Asp6912, Asp6897, Cys6913 and Asp6928 are 98%, 97%, 86%, and 78%, respectively.
Based on molecular dynamics simulation analysis, sinefungin can stably exist in the binding pocket. The analysis of the interaction of the whole trajectory revealed that the hydrogen bonds formed by Asp6912, Cys6913, Asp6897 and Asp6928 with sinefungin were quite conservative ( Fig. 1b; Fig. S2 in Supplementary Materials), and these amino acids were critical for the stable binding of sinefungin as an inhibitor.

Construction of pharmacophore
The key feature of sinefungin stable in the binding pocket can be characterized by pharmacophores based on structure, which can be used in screening simulation to greatly improve the probability of screening  compound with inhibitory activity. In the 30-45ns interval of molecular dynamics simulation, sinefungin maintains a relatively stable state at the binding pocket ( Fig. S1 in Supplementary Materials). The 300 frame structure conformations generated by the selected trajectory was used to cluster according to RMSD value, which is helpful for us to investigate the representative pose that maintains the ligand stability in the pocket and construct the screening pharmacophore with key characteristics. A total of 6 clusters were obtained by cluster analysis. The representative conformations of these 6 clusters were selected for pharmacophore construction. The pharmacophore characteristics included hydrogen bond receptor (A), hydrogen bond donor (D), hydrophobic (H), aromatic ring (R), positive charge (P) and negative charge (N). Based on the structural characteristics of receptor-ligand complexes, 6 pharmacophore models as shown in Fig. 2 were obtained. The characteristic number of 6 pharmacophores varied from 4 to 7. In order to screen common characteristics of six pharmacophores, 6 pharmacophore models were superimposed and placed in binding pocket. As shown in Fig. 3a, the common features are D10, A2, D8, P14 and N13. Among them, the superposition outcomes of D10, A2, D8 and P14 shows that the deviation of the four features in both displacement and direction is very small. These results indicate that the interaction corresponding to the four features (hydrogen bonds with Asp6912, Cys6913, Asp6897 and Asp6928) is highly conservative in the 6 clusters which was consistent with the results of molecular dynamics simulation analysis. N13 has an obvious position shift, N13 is mapped on the carboxylic group at the end of sinefungin structure. The last frame of the dynamics simulation was superimposed with the initial frame. The deflection amplitude and direction during the dynamics simulation are shown in Fig. 3b. Compared with the initial conformation, the overall deviation of the stabilized sinefungin is smaller, but the terminal carboxylic acid part has a significant rotational displacement change. In order to systematically analyze the existing angle deviation, the change of the distance between the hydroxyl group on carboxylic acid and Asn6841 in the molecular dynamics simulation trajectory was monitored. The data (Fig. S3 in Supplementary Materials) showed that the distance between the monitored objects changed memorably (>8 Å), indicating that the electrostatic interaction between the carboxylic acid and Asn6841 corresponding to N13 was extremely unstable. Therefore, we removed the N13 pharmacophore characteristics. In addition, we analyzed the exclusion volume. The flexibility of the amino acid structure around the binding pocket is poor. During the molecular dynamics simulation, the amino acid around the pocket has a micro displacement, so the amino acid that affects the compound binding effect is no longer flexible. Finally, a simplified pharmacophore model retaining the key characteristics was obtained after screening (Fig. 4).
The pharmacophore model has four characteristics, including two hydrogen bond donor D8 and D10, positive charge P14 and a hydrogen bond acceptor A2. These characteristics can reflect the stability of sinefungin and the key factor of its inhibitory activity while also, improving the possibility of screening results with inhibitory activity. At the same time, the simplified pharmacophore model also provides the possibility of diversification of molecular configuration of screening results.

Pharmacophore model screening
The simplified pharmacophore model was used to screen about 1.66 million compounds in Bionet and Chembiv databases. The minimum matching number with pharmacophore model was set to 3. The screening results were sorted according to "phase screen score" and a total of 543 compounds were retained. These candidate compounds were then screened by molecular docking and the interaction was analyzed.

Docking simulation
In total 543 compounds screened by the pharmacophore model could theoretically form the interaction similar to sinefungin. In order to examine whether these compounds fit the binding pocket and form the critical interaction mode with the active site, molecular docking method was used to further screen the pharmacophore screening results. Docking simulation contains Glide SP and XP methods. Before docking experiments, SP and XP were used for re-docking the co-crystalline compounds. As shown in Fig. 5, the RMSD values deviation between the conformation generated by the two docking methods and the original crystal conformation is not greater than 2 Å. In addition, the conformation generated by XP docking result in a more reductive crystal conformation. Because the reproducibility of XP docking is better, the First, we used SP to dock 543 compounds screened by pharmacophores model, of which 527 molecules were docked to the binding pocket. The docking results were ranked by "Docking score". We selected the top 100 docking compounds for XP simulation to reduce the false positive rate of molecular docking. Finally, the top 10 compounds generated from XP simulation were selected to analyze their hydrogen bonding interaction. The docking data and structural formulas of 10 compounds were recorded in Table 2. C2, C3, C4, C7, C8 and C10 can form hydrogen bonds with key amino acids Asp6928 and Asp6897. In addition, the nitrogen ions in these compounds can also form salt bridges with Asp6928. Furthermore ionic interaction can be observed in the π-cation interactions of C3, C4 and C7 with Lys6844. As show in Table 2, the structural difference between C1 and C2 is only in the substituents on the benzene ring, but there is evident difference in the interaction between the C1 and C2 in the binding pocket (Fig. S4 in Supplementary Materials). Moreover, the   distance parameter of the hydrogen bond between the compound and the amino acid raises questions as to whether some of the hydrogen bonds can exist the same time, as shown in Fig. 6. Hence the C1 and C2 conformations produced by the docking simulation were superimposed at the docking pocket. The amino on the double bond side of C2 forms a hydrogen bond with Asp6928, while the position of C1 in the bonding pocket is lower than C2. We observed that there is sufficient space at the end of the binding pocket, which makes it possible for C1 to shift to the position that coincides with the C2. The position change increases the possibility that the amino on the C1 double bond forms a hydrogen bond with Asp6928 (the yellow cutting head direction in Fig. 6 represents this hypothesis). In addition, the amino connected to the two segments of the double bond can undergo a certain degree of rotation, so the amino on the opposite side can be changed by rotation and shorten the distance between Asp6928 form hydrogen bond. Furthermore, the pockets where the two amino linking groups are deeper, which can eliminate the negative effect of steric hindrance when the two amino end of the double bond rotate (the red rotating cutting head in Fig. 6 represents hypothesis). The default protein structure of Glide docking simulation is rigid and insufficient sampling of the ligand would result in a docking pose that lacks accuracy. Therefore, the interaction of ligands at the binding pocket was further analyzed by MM/GBSA and molecular dynamics.

MM/GBSA calculation
The thermodynamic concept was introduced to simulate the process of ligand binding to the protein pocket, which can effectively determine whether the interaction mode between receptor and ligand is stable. Prime MM/GBSA was used to scan and analyze the receptor-ligand binding mode generated by XP docking simulation, and to determine the binding stability of the complex through calculating the free binding energy. As shown in Table 3, the free binding energy (ΔG bind ) of sinefungin was − 51.318 kcal/mol. Among the candidate compounds, the free binding energy of C1, C2, C3 and C7 were above this value, and C2, C3 and C7 can form hydrogen bonds with key amino acid Asp6898 and Asp6928. The similar structures of C1 and C2 results in similar free binding energy values (ΔG bind C1 = − 59.983 kcal/mol; ΔG bind C2 = − 70.916 kcal/mol) which are better than those of the other two compounds. The free binding energy results of C2 show that ΔG bind Coulomd, ΔG bind vdW, ΔG bind Lipo (nonpolar solvation energy) and ΔG bind Hbond provide positive contributions with the main contributions from ΔG bind vdW (− 51.997 kcal/mol) and ΔG bind Coulomd (− 24.058 kcal/ mol). The calculated results of C1, C3 and C7 are similar to C2 and sinefungin. The results for MM/GBSA show that the binding of C1, C2, C3 and C7 with the protein was relatively stable, and the van der Waals energy and Coulomb energy are the main contributors to the free binding energy.

Molecular dynamics simulation
The stability of the interaction between the four compounds and nsp16 was investigated systematically by 150ns molecular dynamics simulation. We expect to use molecular dynamics simulation to investigate whether the four compounds can interact with key amino acids during the simulation, and whether these interactions can persist similar to sinefungin.
As shown in Fig. 7a, during molecular dynamics simulation of the four complexes, the structure of nsp16 showed a stable trend after 105ns, and the fluctuation of RMSD value for C2-protein was significantly weaker than other three groups of compounds. At the same time, the complex system composed of the four compounds and the protein was compared with sinefungin. The deviation amplitude of the protein structure relative to the initial structure at 50ns of molecular dynamics simulation is about the same, and the RMSD value maintained at 2.0-3.5 Å. The reasonable fluctuation range proved that the proteins in the four complexes could maintain a stable state in the simulation process.
In addition, the data of Fig. 7b shows that Stable RMSD curve fluctuations can prove that the four compounds can remain stable in the binding pocket. However, it can be evidently obvious that at the end of  However, the result of C2 is similar to sinefungin with the position at the binding pocket deviating only slightly from the initial docking position after the simulation. In conclusion, at the end of the molecular dynamics simulation, the four compounds maintained a relatively stable state in the pocket, but the interaction of C1, C3 and C7 in the binding pocket could be changed by the obvious position deviation compared with docking pose. By analyzing the histogram of the interaction in the simulation process ( Fig. 8c and d), it evident that the critical hydrogen bonds (interactions fraction was 0.995 and 0.901 respectively) and ionic interactions (0.606 and 0.887 respectively) between C2 and Asp6897, Asp6928 are the main interaction. The RMSD value of C2 at the binding pocket was maintained between 2-3Å, so the critical hydrogen bond with the key amino acid Asp6928 in the docking simulation was retained during the molecular dynamics simulation. In addition, new main hydrogen bonds were formed between the fluorobenzene linked imide group and Asp6897, Tyr6930 (0.754). The change of position in the binding pocket brought a positive effect on C1, which forms critical hydrogen bonds with key amino acids Asp6928 (0.645) and 6897 (0.709). Main hydrogen bonds formed by the carboxylic acid on the C1 benzene ring with Lys6874 (0.577) and the carbonyl oxygen with Asn6899 (0.907) did not appear in the docking simulation ( Fig. 8a and  b). However, for C3 and C7, obvious displacement change brings negative effects. In the docking simulation, both C3 and C7 can form a hydrogen bond with Asp6928 (Table 2), but this hydrogen bond was not stable in the molecular dynamics simulation. Upon stable binding, C3 and C7 deviate from the initial docking position and lose the dominant interaction mode with key amino acids (Fig. 8e, f, g and h). The carboxylic acid at C3 terminal forms strong hydrogen bonds with Leu6855 (0.399) and Thr6854 (1.428). The nitrogen cation on C7 piperazine forms a strong water bridge with Lys7051 (0.904), and the amino groups that linked to two benzene form a strong hydrogen bond with Thr6854 (0.970) and Thr6856 (0.974).
In summary, molecular dynamics simulation verified our doubts about the docking simulation results. In the molecular dynamics simulation, the candidate compounds with superior free binding energy were stable in the binding pocket. The results showed that C1 and C2 could form continuous and stable hydrogen bonds with Asp6897 and Asp6928, while C3 and C7 could not reproduce the hydrogen bonds formed by the docking simulation results with Asp6928. Among the four compounds, the change of RMSD values of protein and ligand during both the molecular dynamics simulation and the interaction mode analysis of C2 was similar to sinefungin C2 also obtained the optimal value in MM/GBSA calculation.

MM/PBSA energy decomposition
The results of Desmond MD simulation showed that the candidate compounds C1 and C2 could form continuous and stable interactions with the key amino acids of nsp16. In order to further explore the energy contribution of these key amino acids in the two candidate compounds and sinefungin, MM/PBSA calculations were performed on the complex system composed of the three compounds and nsp16 after Gromacs MD simulation. The RMSD values of the Cα atoms were monitored along the whole MD trajectory, as shown in Fig. 9 All the trajectories indicated that the MD simulations arrived at equilibrium state. MM/PBSA binding free energy calculation results show that the candidate compound and sinefungin can stably bind to the nsp16 binding pocket (Fig. S5). Fig. 10 shows the results of amino acid energy decomposition. Amino acids whose energy value is lower than − 0.5 kcal/mol were used for analysis. During the binding of sinefungin to nsp16, the polar interaction (ΔG electrostatic and ΔG Polar Solvation ) affects the binding stability to the greatest extent. The histogram shows that the key amino acids (Asp6912, Cys6913, Asp6897 and Asp6928) all have favorable energy contribution, and the electrostatic potential energy contributions are more prominent than other amino acids. This suggests that there are Note: a free binding energy; B the contribution of Coulomb on the free binding energy; c the contribution of covalent on free binding energy; d the contribution of hydrogen bond on free binding energy; e the contribution of lipophilic on free binding energy; f the contribution of van der Waals on the free binding energy; g contribution of stacking on binding energy; h the contribution of polar solvation on free binding energy. strong hydrogen bonds between sinefungin and these amino acids. Asp6897 and Asp6928 also have crucial electrostatic potential energy contributions in the binding process of C1 and C2, and the total energy contribution is greather than other amino acids. Similar to sinefungin, Asp6897 produced the ultimate contribution value. The results of amino acid energy decomposition further verified the rationality involved in the determination of key amino acids. Candidate compounds C1 and C2 and sinefungin can stably bind to the docking pocket, and both form strong hydrogen bonds with the key amino acids Asp6897 and Asp6928.

Conclusion
In this work, we constructed a pharmacophore model that can accurately reflect the crucial interaction mode between sinefungin and nsp16. Molecular dynamic simulation was integrated into the structurebased pharmacophore construction, and the pharmacophore model was used to screen Bionet and Chembiv databases in order to find novel and potent nsp16 inhibitors. Through docking simulation, MM/GBSA, molecular dynamics RMSD stability analysis and main interaction analysis, it is concluded that C1 and C2 can stably bind to the binding pocket and form strong and continuous hydrogen bonds with the key amino acids Asp6897/Asp6928. MM/PBSA energy decomposition of sinefungin showed that Asp6897 and Asp6928 contributed major energy values during C1 and C2 binding, which was similar to sinefungin. Therefore C1 and C2 can be used as nsp16 inhibitors for further design and development.