Discovery of Potential SARS-CoV-2 M Protease Inhibitors by Virtual Screening, Molecular Dynamics, and Binding Free Energy Analyses

1. INTRODUCTION

In December 2019, patients with acute respiratory syndrome of unknown origin were found. The signs of infection that they showed included fever, cough, shortness of breath, and dyspnea^[1]. In severe cases, the lung infections can even lead to pneumonia, kidney failure, and death^{[2, 3]}. The disease caused by the novel coronavirus has been officially named COVID-19 by the World Health Organization. Due to the rapid spread of COVID-19, the World Health Organization has defined COVID-19 as a risk of all indigenous peoples globally.The Centers of Disease Control and Prevention (CDC) and the National Health Commission of the People's Republic of China officially announced that the cause of the pneumonia outbreak was a novel coronavirus that has the characteristics of a typical coronavirus family, particularly of the SARS-CoV^{[4, 5]}.

Coronavirus (CoV) is a positive-sense RNA virus with an envelope that allows the transmission of zoonotic diseases^{[4, 6]}. This virus has a prominent spike protein on its envelope, which has a coronal or halo-like appearance^{[7, 8]}. There are six types of coronaviruses that have been found to infect humans around the world: HCoV-229E, HCoV-OC43, SARS, HCoV-NL63, HKU1, and MERS^[9]. SARS-CoV-2 is a beta coronavirus similar to MERS-CoV and SARS-CoV, which both infect the lower respiratory tract, and have a high degree of infectivity and a certain mortality^{[10, 11]}.

Currently, some diagnostics, vaccines and therapeutics that have therapeutic effects on COVID-19 have been reported one after another^[12], and effective therapeutic drugs are still being explored. Fortunately, the genome of SARS-CoV-2 has been sequence^{[5, 13]}. Moreover, Rao et al. published a high-resolution crystal structure of the 3-chymotrypsin-like hydrolase protease (Mpro) of SARS-CoV-2 in the PDB database; this provides a structural basis of in silico virtual screening of the potential for therapeutic drugs, active natural products and traditional Chinese medicine. In this study, Mpro was used as the target, and molecular docking methods were used to virtually screen candidate drugs from the marketed or clinical drug database (DrugBank). The possible clinical applications of the candidate drug were then evaluated by using molecular dynamics (MD) simulations and binding free energy analysis. Therefore, it is desirable to identify a drug with good anti-SARS-COV-2 activity.

2. MATERIALS AND METHODS

2.1 SARS-CoV-2 Mpro 3D structure preparation

The crystal structure of 3-chymotrypsin-like hydrolase protease (Mpro, PDB ID: 6LU7) of SARS-CoV-2 was first obtained from the protein database (www.rcsb.org). We then removed the co-crystal ligand and all water molecules and added hydrogen atoms to the protein. Finally, SYBYL-X2.1 was used for the residue repair and energy minimization.

2.2 Virtual screening by SARS-CoV-2 Mpro 3D structure

Based on the crystal structure of Mpro, the binding site was defined by co-crystal ligand. The molecular docking program Surflex-Dock GeomX (SFXC) in SYBYL-X 2.1 was utilized to carry out the virtual screening and construction of Mpro/drugs complex in which the docking scores were expressed in –lg(Kd)^[14]. The main protocols or parameters of docking were listed as follows^[15-19]: (a) The "number of starting conformations per ligand" was set to 10, and the "number of max conformations per fragment" was set to 20; (b) The "maximum number of rotatable bonds per molecule" was set to 100; (c) Flags were turned on for "pre-dock minimization", "post-dock minimization", "molecule fragmentation", and "soft grid treatment"; (d) "Activate spin alignment method with density of search" was set to 9.0; and (e) The "number of spins per alignment" was set to 12.

Specifically, the exact protocol employed for the hits' selection of virtual screening was as follows: (a) First, 10 starting conformations for each antiviral drug dataset were generated; (b) Each compound was then docked with the binding pocket in Mpro using the docking algorithm implemented in SYBYL and the top 20 binding poses of each compound with higher docking scores were saved; (c) The best docking score of each peptide was used for the sorting; and (d) The top three compounds with the highest docking scores were selected as the hit compounds for further analyses^[20].

2.3 Molecular dynamics (MD) simulation

The Mpro complexed with drug candidates were set up for MD simulations. For example, the complexes were put into 0.15 M NaCl solution to a cubic water box, which included 29,156 water molecules, 84 Na⁺ ions, and 79 Cl⁻ ions. The initial configurations of the receptors and ligands were taken from docking studies. The sizes of the initial simulation box were ~103Å × 103Å × 103Å. The other systems were set up with the same protocol.

The MD simulations were carried out using the PMEMD.mpi and PMEMD.cuda modules in the AMBER16^{[21, 22]} package. First, five minimization steps were conducted for the systems to avoid possible steric crashes. Then, each system was gradually heated from –273.15 ℃ to 26.85 ℃ during the heating stage and maintained at 26.85 ℃ during the subsequent equilibrium and production stages. A time step of 2 fs was used for the heating stage, equilibrium stage, density adjustment and the entire production stage. A periodic boundary condition was employed to maintain constant temperature and the pressure ensembles. The pressure was set to 1 atm and controlled by the anisotropic (x-, y-, z-) pressure scaling protocol with a pressure relaxation time of 1 ps. The temperature was regulated using Langevin dynamics with a collision frequency of 2 ps^{−1[23, 24]}. The particle mesh Ewald (PME) method^{[25, 26]} was adopted to handle long-range electrostatics and a 10 Å cutoff was set to treat real-space interactions. All covalent bonds involving hydrogen atoms were constrained with the SHAKE algorithm^[27]. Each system was subject to a 100 ns MD simulation and the trajectory of simulated systems was saved every 100 ps^[28].

2.4 Molecular mechanics/generalized Born surface area (MM/GBSA) calculation

For the saved trajectories of MD simulations, the MM/GBSA^[29-35] method was used to calculate the binding energies of receptors treated with different ligands. A total of 200 snapshots were extracted from 80 to 100 ns to calculate the mean binding energy. The formula is as follows:

Where ΔE_bind is the binding energy and ΔE_MM denotes the sum of the molecular mechanical energies in a vacuum and can be further divided into the contributions to electrostatic, van der Waals, and internal energies. This term could be computed using the molecular mechanics method. ΔE_SOL is the solvation energy, which includes the polar solvation energy (ΔE_GB) calculated with the generalized born (GB) approximation model^{[36, 37]} and the non-polar part (ΔE_SA) obtained by fitting solvent accessible surface area (SASA)^[38] with the linear combinations of pairwise overlaps (LCPO) model^{[39, 40]}. Additionally, the energies of each residue were decomposed into the backbone and side-chain atoms. The energy decomposition can be analyzed to determine the contributions to the key residues to the binding^{[28, 41]}.

3. RESULTS AND DISCUSSION

3.1 Potential binding site on the Mpro interface

To verify the reliability and feasibility of the docking program, the co-crystallized ligand was redocked into a binding pocket by using the Surflex-Dock program. As shown in Fig. 1, the redocked conformation and co-crystallized conformation overlap well, and indicate that the docking method can be used virtually to screen the potential drug candidates. Here, the crystal structure of Mpro was obtained from the Zihe Rao and Haitao Yang's team at the ShanghaiTech University. According to the crystal structure of Mpro-N3 complex, the binding pocket of Mpro was composed of the following residues: His41, Met49, Asn142, Gly143, Met165, Glu166, Pro168, Gln189, Thr190, Ala191 and Gln192.

Figure 1

Figure 1.  Mpro structure and binding site

DownLoad: Full-Size Img PowerPoint

3.2 Analysis of the virtual screening results

In this study, the binding pocket of Mpro was used as the active site to which the virtual screening was performed. The molecular docking screening of antiviral drug datasets was performed by using Surflex-Dock GeomX method. Based on the docking scores of candidate drugs and their matching degree with the co-crystallized ligand, 9 drug candidates with high docking scores of SARS-CoV-2 Mpro were screened. The binding free energies between Mpro and candidates were less than –10 kcal·mol^-1 (Table 1). Interestingly, among these drug candidates, Ombitasvir, Elbasivr, Lopinavir, Tipranavir and Darunavir were approved treatments for hepatitis C virus (HCV) and human immunodeficiency virus (HIV), whereas Remdesivir, GS-6620 and Rupintrivir have been investigated for the treatment of HCV, SARS-CoV and human rhinoviral (HRV) infection. Additionally, Adafosbuvir has been used to evaluate the effect of renal impairment on the pharmaco-kinetics of prodrug AL-335 (investigational HCV nonstructural protein 5B inhibitor).

Table 1

Table 1.  Drug Candidates Screened from DrugBank by Molecular Docking and Free Energy

DownLoad: CSV

No.	Name	Structure	Score	ΔG (kcal/mol)	Treatmet	Status
1	Co-crystal ligand N3		9.00	–36.29
2	Ombitasvir		9.66	–31.83	HCV	Approved
3	Lopinavir		11.21	–24.75	HIV	Approved
4	Tipranavir		9.76	–23.62	HIV	Approved
5	Darunavir		9.66	–14.60	HIV	Approved
6	Elbasvir		9.07	–13.98	HCV	Approved
7	Rupintrivir		9.38	–23.65	HRV	Investigational
8	Remdesivir		9.38	–18.57	SARS-CoV	Experimental
9	GS-6620		8.60	–36.93	HCV	Investigational
10	Adafosbuvir		9.55	–25.20	Renal impairment	Investigational

3.3 Molecular dynamics simulations and binding free energy analysis

3.3.1   Interaction between Mpro and Ombitasvir

To further verify the reliability of the drug candidates, we analyzed the interactions between Mpro and the drugs through MD simulation and binding free energy analysis. Among the candidate drugs, Ombitasvir is a drug marketed for the treatment of Hepatitis C. At present, there have been reports^[42] that Ombitasvir also has a potential inhibitory effect on SARS-CoV-2. Therefore, the Mpro/Ombitasvir complex was constructed, and a 100 ns MD simulation was repeated three times. The root mean square deviation (RMSD) of the bipartite molecular dynamics simulations of the Mpro/Ombitasvir complex is shown in Fig. S1 of the supplementary material. Fig. S1a shows that the complex molecule was stable after 60 ns; the Ombitasvir's RMSD fluctuates around 3.7 Å, and the RMSD of the receptor main chain atoms fluctuates around 1.8 Å. In the RMSD of the second molecular dynamics simulation (Fig. S1b), after 30 ns, the RMSD of Ombitasvir stabilized around ~4.2 Å, and the RMSD of the receptor main chains atoms increased from 1.2 to 3.0 Å at 10~30 ns; at 30~50 ns, the RMSD decreased from 3.0 to around 2.5 Å, then stabilized around 3.4 Å until the end. The RMSD of MD simulations was stable to varying degrees and these results indicate that Ombitasvir and the receptor main chain atoms have reached a stable conformation after a large conformational change in the simulation.

To discover the binding mode between Ombitasvir and Mpro, we thoroughly analyzed the interaction mode between Ombitasvir and Mpro key amino acids after bipartite molecular dynamics simulation studies. The bonding modes between Mpro and Ombitasvir were similar (Fig. 2). The formyl group adjacent to the pyrrole rings formed two hydrogen bonds between Gln192, and the bond lengths were 3.3 and 3.2 Å, respectively. Concurrently, the benzene ring formed a strong hydrophobic interaction with Met49 (3.6 Å) through van der Waals forces, and tert-butyl formed a weak hydrophobic bond with Met165 at a distance of 4.0 Å. Interestingly, Fig. S2 shows that the nitrogen and oxygen atoms of the carbamoyl group produced two strong hydrogen bonds with Asn142 (3.1 Å) and Gly143 (2.9 Å), and the formyl group adjacent to the pyrrole ring and the oxygen atom of the carbamate formed two hydrogen bonds with Gln192 (3.2 and 3.2 Å), which could promote the binding of Ombitasvir to Mpro.

Figure 2

Figure 2.  Analysis of the binding mode of Ombitasvir to Mpro (the first molecular dynamics simulation results. Key residues are shown in cyan and Ombitasvir in blue. The hydrogen and hydrophobic bonds are shown as red and blue dashed lines, respectively)

DownLoad: Full-Size Img PowerPoint

The energy contribution to the key residues of Mpro and Ombitasvir might give more information for their binding and guide its structural modification. We calculated the interactional energy for the binding of Ombitasvir to the key residues of Mpro by using molecular mechanics/generalized-Born surface area (MM/GBSA) method. The results are shown in Fig. 3a and Fig. S3 (Data are listed in Supplementary Table S2). The results show that the interactional energy between Ombitasvir and His41, Met49, Asn142, Pro168, Gln189, and Ala191 was low, which was mainly from the Van der Waals force, and thereby indicates that the affinity of Mpro/Ombitasvir was mainly from hydrophobic interactions. The interaction energy between Gln189 residue and Ombitasvir was lower than –10.0 kcal·mol^-1, which implies that it is important to their binding. Aside from van der Waals force, the electrostatic and non-polar interactional energy of Gly143, Met165, Glu166, Thr190, and Gln192 had a moderate contribution to binding affinity. Subsequently, we determined the hydrogen bonds between the last 20ns (200 conformations) of MD simulation; the results are shown in Fig. 3b. Fig. 3b shows that the hydrogen bonds were formed into Asn124, Gly134, Gln192 and Ombitasvir; however, their total interactional energy were relatively low, which validates that the binding energy was mainly from hydrophobic interactions.

Figure 3

Figure 3.  Mpro/Ombitasvir system energy decomposition analysis results and the number of hydrogen bonds in key amino acid residues. (a) Free energies for the key residues of Mpro calculated by using MM/GBSA. (b) Number of hydrogen bonds in the last 20 ns of the 100 ns MD simulation (Mpro/Ombitasvir system)

DownLoad: Full-Size Img PowerPoint

Based on the linear relationship between affinity and free binding energy, the affinity with Mpro/Ombitasvir might be accurately predicted by the ΔE_bind energy. Here, the binding free energy of Mpro/Ombitasvir was calculated by using MM/PBSA method (listed in Table 2). The free energies were –31.8315 and –31.0475 kcal·mol^–1 in the three MD simulations, respectively, with their mean energy to be –31.44 kcal·mol^–1. Although the MD simulation for complex Mpro/Ombitasvirwere was performed in bipartite randomly, the key residues of Mpro that bound Ombitasvir were almost identical, and these residues promoted binding affinity with hydrophilic and hydrophobic interactions. Interestingly, with the contributions of hydrogen bonds and hydrophobicity, the binding energies of Mpro/Ombitasvir in two MD simulations were highly similar. The results showed not only the reliability of the calculation in this study but also the good binding affinity of Ombitasvir and Mpro. Thus, Ombitasvir has good potential anti-SARS-COV-2 activity.

Table 2

Table 2.  Binding Free Energies (ΔG_bind) Calculated by Using MM/PBSA (kcal·mol^-1)

DownLoad: CSV

Compounds	MD1	MD2	Mean
Ombitasvir	–31.83	–31.05	–31.44
Lopinavir	–23.84	–24.75	–24.23
GS-6620	–36.93	–31.59	–34.26
Remdesivir	–18.57	–14.84	–16.71

2.3.2   Interaction between Mpro and Lopinavir

Lopinavir is an antiretroviral protease inhibitor used in combination with other antiretrovirals in the treatment of HIV-1 infection and reported certainly inhibitory effect of SARS-CoV-2^[43]. We explored the binding mode between Mpro Lopinavir through MD simulations. The RMSD plots (Fig. S4a) show the main chain of Mpro and Lopinavir stabilized at 1.8 and 3.6 Å after 20 ns, respectively. Similarly, the repeated MD simulations reached the stable state after 70 and 20 ns (Fig. S4b), which might be used to analyze their binding mode.

The stable conformations of MD simulation (the last 20 ns) were extracted to analyze the binding mode of Mpro/Lopinavir. Lopinavir complex about to Mpro through binding to key residues was identified with Mpro/Ombitasvir. As shown in Fig. 4, the pyrimidine ring and oxygen atom in butyramide of Lopinvar formed hydrogen bonds with Gln189, and the distances are 3.2 Å and 2.9 Å, respectively. Meanwhile, the pyrimidine and benzene rings could interact with Met165 (3.8 Å) and Met49 (3.3 Å) through hydrophobic bonds correspondingly. From repeated MD simulation (Fig. S5b), the nitrogen and oxygen atoms of acetamido could form two hydrogen bonds between Asn142 (3.0 Å) and Gly143 (2.9 Å), respectively, and the benzene ring has hydrophobic interaction with Met49 (3.4 Å), Met165 (3.6 Å) and Glu166 (3.5 Å).

Figure 4

Figure 4.  Analysis of the binding mode of Lopinavir to Mpro (the first MD simulation results. Key residues are shown in cyan and Lopinavir is shown in blue. The hydrogen and hydrophobic bonds are shown as red and blue dashed lines, respectively)

DownLoad: Full-Size Img PowerPoint

Based on the stable conformation, the energy decomposition for key residues was performed to its contribution to Mpro combined with Lopinavir. The energy contribution to key residues is shown in Fig. 5 and Fig. S6 (Data are listed in Supplementary Table S2). Fig. 5 and Fig. S6 show that lopinavir had strong interactions with His41, Met49, Met165, Glu166, Gln189 and Gln192 of Mpro through van de Waals forces. In addition, Asn142, Gly143, Pro168, Thr190 and Ala191 had lower energy contribution for binding. According to the statistics of hydrogen bonds between the last 20 ns of MD simulation (Fig. 5b), although the conformation of Mpro/Lopinavir complex is slightly changed, Lopinavir still interacts with these key residues. Subsequently, we use the MM/PBSA method to calculate its binding free energy (Table 2). The results of the bipartite calculations were not much different (–23.84 and –24.75 kcal·mol^-1), which shows the key residues of Lopinavir and Mpro were tightly bound under the interaction of hydrogen and hydrophobic bonds, so that Lopinavir had a good anti-SARS-COV-2 activity.

Figure 5

Figure 5.  Mpro/Lopinavir system energy decomposition analysis results and the number of hydrogen bonds in key amino acid residues. (a) Free energies for the key residues of Mpro calculated by using MM/GBSA. (b) Number of hydrogen bonds in the last 20 ns of the 100 ns MD simulation (Mpro/Lopinavir system)

DownLoad: Full-Size Img PowerPoint

3.3.3   Interaction between Mpro and GS-6620

We found that among the screened drug candidates, GS-6620, which is still in the investigational state, has a high binding free energy. Therefore, the Mpro/GS-6620 complex was constructed, and a 100 ns molecular dynamics simulation was conducted twice. Fig. S7 (see Supplementary Material Fig. S7) shows the RMSD of the MD simulation of the Mpro/GS-6620 complex. Fig. S7a shows that the complex about Mpro/GS-6620 is stable after 80 ns, and the RMSD of the main chain of Mpro and GS-6620 fluctuates around 2.2 Å and 3.0 Å, respectively. For repeated MD simulation, the complex reached a stable state after 2 ns (Fig. S7b), which might guarantee the effect of the binding mode analysis.

The interaction between GS-6620 and Mpro is shown in Fig. 6; the hydroxyl and amino groups on the GS-6620 phosphoramide bond formed two strong hydrogen bonds between Glu166 (3.8 and 2.9 Å). The isopropyl and benzene rings of GS-6620 also produced a hydrophobic interaction with Met49. In addition, the benzene rings in the molecule also formed a hydrophobic interaction with Met165. The repeated MD simulation showed that the oxygen atoms of GS-6620 might form five hydrogen bonds between His41, Asn142, Gly143, and Gln189, and the benzene ring produced hydrophobic interactions with both Met49 and Met165 (Fig. S8a); this increased the binding affinity with GS-6620 with Mpro.

Figure 6

Figure 6.  Analysis of the binding mode of GS-6620 to Mpro (the first molecular dynamics simulation results. Key residues are shown in cyan and GS-6620 in blue. The hydrogen and hydrophobic bonds are shown as red and blue dashed lines, respectively)

DownLoad: Full-Size Img PowerPoint

We calculated the interactional energy of GS-6620 with the key amino acids of Mpro by using MM/GBSA (Fig. 7a and Fig. S9). The results showed that GS-6620 had strong interaction with the Met49, Asn142, Gly143, and Met165 of Mpro. The statistics of hydrogen bonds are shown in Fig. 7b and supports the binding of GS-6620 binding with His41, Asn142, Gly143 Glu66 and Gln189. Subsequently, we used the MM/PBSA method to accurately calculate the binding free energy of GS-6620 and Remdesivir with Mpro. The results showed that the binding free energy of Remdesivir (–16.71 kcal·mol^-1) is relatively higher and may be related to its initial conformation, or that it may not directly affect the target^[44]; the relevant results need further verification. It is surprising that the binding energy of GS-6620 with Mpro was comparable to that of the Ombitasvir (–34.26 vs. –31.44 kcal·mol^-1). These binding energies indicate that GS-6620 has good potential anti-SARS-COV-2 activity.

Figure 7

Figure 7.  Mpro/GS-6620 system energy decomposition analysis results and the number of hydrogen bonds in key amino acid residues. (a) Free energies for the key residues of Mpro calculated by using MM/GBSA. (b) Number of hydrogen bonds in the last 20 ns of the 100 ns MD simulation (Mpro/GS-6620 system)

DownLoad: Full-Size Img PowerPoint

4. CONCLUSION

The pneumonia caused by SARS-CoV-2 is highly infectious and pathogenic. It has caused two thousand deaths since its outbreak of December 2019. There have not yet been effective preventions or treatment methods for COVID-19. In this study, we combined virtual screening, molecular docking, and MD simulations to explore the potential anti-SARS-COV-2 activity of drug candidates in the DrugBank database. 9 drug candidates that can bind the Mpro of SARS-CoV-2 in silico were screened. Based on the stable conformation determined from MD simulation, the binding mode accompanied with free energy decomposition showed that His41, Met49, Asn142, Met165 and Gln189 played an important role for the binding of Mpro with the drug candidates. Moreover, the binding free energies of Ombitasvir, Lopintavir and GS-6620 to Mpro were relatively low. These findings indicate that Ombitasvir, Lopintavir and GS-6620 might have good anti-SARS-COV-2 activity.

References

1. [1]

   Bassetti, M.; Vena, A.; Roberto Giacobbe Giacobbe, D. R. The novel Chinese coronavirus (2019-nCoV) infections: challenges for fighting the storm. Eur. J. Clin. Invest. 2020, e13209−4.
   

2. [2]

   Chen, N.; Zhou, M.; Dong, X.; Qu, J.; Gong, F.; Han, Y.; Qiu, Y.; Liu, Y.; Wei, Y.; Xia, J. A.; Yu, T.; Zhang, X. X.; Zhang, L. Epidemiological and clinical characteristics of 99 cases of 2019 novel coronavirus pneumonia in Wuhan, China: a descriptive study. Lancet 2020, 395, 507−513.  doi: 10.1016/S0140-6736(20)30211-7
   

3. [3]

   Huang, C. L.; Wang, Y. M.; Li, X. W.; Ren, L. L.; Zhao, J. P.; Hu, Y.; Zhang, L.; Fan, G. H.; Xu, J. Y.; Gu, X. Y.; Cheng, Z. S.; Yu, T.; Xia, J. A.; Wei, Y.; Wu, W. J.; Xie, X. L.; Yin, W.; Li, H.; Cao, B. Clinical features of patients infected with 2019 novel coronavirus in Wuhan, China. Lancet 2020, 395, 497−506.  doi: 10.1016/S0140-6736(20)30183-5
   

4. [4]

   Hui, D. S.; Azhar, E. I.; Madani, T. A.; Ntoumi, F.; Kock, R.; Dar, O.; Ippolito, G.; Mchugh, T. D.; Memish, Z. A. Arabia, H. R. S.; Drosten, C.; Germany, B.; Zumla, A.; Petersen, E.; Zumla, A. The continuing 2019-nCoV epidemic threat of novel coronaviruses to global health – The latest 2019 novel coronavirus outbreak in Wuhan, China. Int. J. Infect. Dis. 2020, 91, 264−266.  doi: 10.1016/j.ijid.2020.01.009
   

5. [5]

   Xu, X.; Chen, P.; Wang, J.; Feng, J.; Zhou, H.; Li, X.; Zhong, W.; Hao, P. Evolution of the novel coronavirus from the ongoing Wuhan outbreak and modeling of its spike protein for risk of human transmission. Sci. China Life Sci. 2020, 63, 457−460.  doi: 10.1007/s11427-020-1637-5
   

6. [6]

   Gallagher, T. M.; Buchmeier, M. J. Coronavirus spike proteins in viral entry and pathogenesis. Virology 2001, 279, 371−374.  doi: 10.1006/viro.2000.0757
   

7. [7]

   Liu, X.; Wang, X. J. Potential inhibitors against 2019-nCoV coronavirus M protease from clinically approved medicines. J. Genet. Genom. 2020, 119–121.
   

8. [8]

   Seah, I.; Su, X.; Lingam, G. Revisiting the dangers of the coronavirus in the ophthalmology practice. Nature 2020, 1−3.
   

9. [9]

   Su, S.; Wong, G.; Shi, W.; Liu, J.; Lai, A. C.; Zhou, J.; Liu, W. J.; Bi, Y. H.; Gao, G. F. Epidemiology, genetic recombination, and pathogenesis of coronaviruses. Trends Microbiol. 2016, 24, 490−502.  doi: 10.1016/j.tim.2016.03.003
   

10. [10]

    Wang, W.; Tang, J.; Wei, F. Updated understanding of the outbreak of 2019 novel coronavirus (2019-nCoV) in Wuhan, China. J. Med. Virol. 2020, 92, 441−447.  doi: 10.1002/jmv.25689
    

11. [11]

    Pang, J.; Wang, M. X.; Ang, I. Y. H.; Tan, S. H. X.; Lewis, R. F.; Chen, J. I. P.; Gutierrez, R. A.; Gwee, S.; X.; W.; Chua, P. E. Y.; Yang, Q.; Ng, X. Y.; Yap, R. K. S.; Tan, H. Y.; Teo, Y. Y.; Tan, C. C.; Cook, A. R.; Yap, J. C. H.; Hsu, L. Y. Potential rapid diagnostics, vaccine and therapeutics for 2019 novel coronavirus (2019-nCoV): a systematic review. J. Clin. Med. 2020, 9, 623−623.  doi: 10.3390/jcm9030623
    

12. [12]

    Andersen, K. G.; Rambaut, A.; Lipkin, W. I.; Holmes, E. C. The proximal origin of SARS-CoV-2. Nat. Med. 2020, 26, 450−452.  doi: 10.1038/s41591-020-0820-9
    

13. [13]

    Jain, A. N. Scoring noncovalent protein-ligand interactions: a continuous differentiable function tuned to compute binding affinities. J. Comput. Aid. Mol. Des. 1996, 10, 427−440.  doi: 10.1007/BF00124474
    

14. [14]

    Feng, Z.; Pearce, L. V.; Xu, X.; Yang, X.; Yang, P.; Blumberg, P. M.; Xie, X. Q. Structural insight into tetrameric hTRPV1 from homology modeling, molecular docking, molecular dynamics simulation, virtual screening, and bioassay validations. J. Chem. Inf. Model. 2015, 55, 572−588.  doi: 10.1021/ci5007189
    

15. [15]

    Chen, J. Z.; Wang, J.; Xie, X. Q. GPCR structure-based virtual screening approach for CB2 antagonist search. J. Chem. Inf. Model. 2007, 47, 1626−1637.  doi: 10.1021/ci7000814
    

16. [16]

    Sheng, S.; Wang, J.; Wang, L.; Liu, H.; Li, P.; Liu, M.; Long, C. F.; Xie, C. S.; Xie, X. Q.; Su, W. Network pharmacology analyses of the antithrombotic pharmacological mechanism of Fufang Xueshuantong Capsule with experimental support using disseminated intravascular coagulation rats. J. Ethnopharmacol. 2014, 154, 735−744.  doi: 10.1016/j.jep.2014.04.048
    

17. [17]

    Feng, Z. W.; Kochanek, S.; Close, D.; Wang, L. R.; Srinivasan, A.; Almehizia, A. A.; Iyer, P.; Xie, X. Q.; Johnston, P. A.; Gold, B. Design and activity of AP endonuclease-1 inhibitors. J. Chem. Biol. 2015, 8, 79−93.  doi: 10.1007/s12154-015-0131-7
    

18. [18]

    Feng, Z. W.; Pearce, L. V.; Zhang, Y.; Xing, C.; Herold, B. K.; Ma, S.; Hu, Z.; Turcios, N. A.; Yang, P.; Tong, Q.; Blumberg, P. M.; McCall, A. K. Multi-functional diarylurea small molecule inhibitors of TRPV1 with therapeutic potential for neuroinflammation. AAPS J. 2016, 18, 898−913.  doi: 10.1208/s12248-016-9888-z
    

19. [19]

    Wang, Y. Q.; Guo, H. Q.; Feng, Z. W.; Wang, S. Y.; Wang, Y. X.; He, Q. X.; Li, G. P.; Lin, W. W.; Xie, X. Q.; Lin, Z. H. PD-1-targeted discovery of peptide inhibitors by virtual screening, molecular dynamics simulation, and surface plasmon resonance. Molecules 2019, 24, 3784−15.  doi: 10.3390/molecules24203784
    

20. [20]

    Salomon-Ferrer, R.; Götz, A. W.; Poole, D.; Le Grand, S.; Walker, R. C. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 2. Explicit solvent particle mesh Ewald. J. Chem. Theory Comput. 2013, 9, 3878−3888.  doi: 10.1021/ct400314y
    

21. [21]

    Gotz, A. W.; Williamson, M. J.; Xu, D.; Poole, D.; Le Grand, S.; Walker, R. C. Routine microsecond molecular dynamics simulations with AMBER on GPUs. 1. Generalized born. J. Chem. Theory Comput. 2012, 8, 1542−1555.  doi: 10.1021/ct200909j
    

22. [22]

    Loncharich, R. J.; Brooks, B. R.; Pastor, R. W. Langevin dynamics of peptides: the frictional dependence of isomerization rates of N-acetylalanyl-N΄-methylamide. Biopolymers 1992, 32, 523−535.  doi: 10.1002/bip.360320508
    

23. [23]

    Izaguirre, J. A.; Catarello, D. P.; Wozniak, J. M.; Skeel, R. D. Langevin stabilization of molecular dynamics. J. Phys. Chem. 2001, 114, 2090−2098.  doi: 10.1063/1.1332996
    

24. [24]

    Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: an N⋅log(N) method for Ewald sums in large systems. J. Phys. Chem. 1993, 98, 10089−10092.  doi: 10.1063/1.464397
    

25. [25]

    Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Phys. Chem. 1995, 103, 8577−8593.  doi: 10.1063/1.470117
    

26. [26]

    Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 1977, 23, 327−341.  doi: 10.1016/0021-9991(77)90098-5
    

27. [27]

    Wang, Y. Q.; Lin, W. W.; Wu, N.; Wang, S. Y.; Chen, M. Z.; Lin, Z. H.; Xie, X. Q.; Feng, Z. W. Structural insight into the serotonin (5-HT) receptor family by molecular docking, molecular dynamics simulation and systems pharmacology analysis. Acta Pharmacol. Sin. 2019, 40, 1138−1156.  doi: 10.1038/s41401-019-0217-9
    

28. [28]

    Wang, J.; Hou, T. Develop and test a solvent accessible surface area-based model in conformational entropy calculations. J. Chem. Inf. Model. 2012, 52, 1199−1212.  doi: 10.1021/ci300064d
    

29. [29]

    Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium. J. Phys. Chem. 1996, 100, 19824−19839.  doi: 10.1021/jp961710n
    

30. [30]

    Wang, W.; Kollman, P. A. Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model. J. Mol. Biol. 2000, 303, 567−582.  doi: 10.1006/jmbi.2000.4057
    

31. [31]

    Sun, H.; Duan, L.; Chen, F.; Liu, H.; Wang, Z.; Pan, P.; Zhu, F.; John Zhang, Z. H.; Hou, T. J. Assessing the performance of MM/PBSA and MM/GBSA methods. 7. Entropy effects on the performance of end-point binding free energy calculation approaches. Phys. Chem. Chem. Phys. 2018, 20, 14450−14460.  doi: 10.1039/C7CP07623A
    

32. [32]

    Chen, F.; Liu, H.; Sun, H.; Pan, P.; Li, Y.; Li, D.; Hou, T. Assessing the performance of the MM/PBSA and MM/GBSA methods. 6. Capability to predict protein-protein binding free energies and re-rank binding poses generated by protein-protein docking. Phys. Chem. Chem. Phys. 2016, 18, 22129−22139.  doi: 10.1039/C6CP03670H
    

33. [33]

    Sun, H.; Li, Y.; Shen, M.; Tian, S.; Xu, L.; Pan, P.; Guan, Y.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 5. Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring. Phys. Chem. Chem. Phys. 2014, 16, 22035−22045.  doi: 10.1039/C4CP03179B
    

34. [34]

    Sun, H.; Li, Y.; Tian, S.; Xu, L.; Hou, T. Assessing the performance of MM/PBSA and MM/GBSA methods. 4. Accuracies of MM/PBSA and MM/GBSA methodologies evaluated by various simulation protocols using PDBbind data set. Phys. Chem. Chem. Phys. 2014, 16, 16719−16729.  doi: 10.1039/C4CP01388C
    

35. [35]

    Tsui, V.; Case, D. A. Theory and applications of the generalized Born solvation model in macromolecular simulations. Biopolymers 2000, 56, 275−291.  doi: 10.1002/1097-0282(2000)56:4<275::AID-BIP10024>3.0.CO;2-E
    

36. [36]

    Bashford, D.; Case, D. A. Generalized born models of macromolecular solvation effects. Annu. Rev. Phys. Chem. 2000, 51, 129−152.  doi: 10.1146/annurev.physchem.51.1.129
    

37. [37]

    Sitkoff, D.; Sharp, K. A.; Honig, B. Accurate calculation of hydration free energies using macroscopic solvent models. J. Phys. Chem. 1994, 98, 1978−1988.  doi: 10.1021/j100058a043
    

38. [38]

    Still, W. C.; Tempczyk, A.; Hawley, R. C.; Hendrickson, T. Semianalytical treatment of solvation for molecular mechanics and dynamics. J. Chem. Sci. 1990, 112, 6127−6129.
    

39. [39]

    Weiser, J.; Shenkin, P. S.; Still, W. C. Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO). J. Comput. Chem. 1999, 20, 217−230.  doi: 10.1002/(SICI)1096-987X(19990130)20:2<217::AID-JCC4>3.0.CO;2-A
    

40. [40]

    Hu, J.; Feng, Z.; Ma, S.; Zhang, Y.; Tong, Q.; Alqarni, M. H.; Gou, X.; Xie, X. Q. Difference and influence of inactive and active states of cannabinoid receptor subtype CB2: from conformation to drug discovery. J. Chem. Inf. Model. 2016, 56, 1152−1163.  doi: 10.1021/acs.jcim.5b00739
    

41. [41]

    Ekins, S.; Mottin, M.; Ramos, P. R.; Sousa, B. K.; Neves, B. J.; Foil, D. H.; Zorn, K. M.; Braga, R. C.; Coffee, M.; Southan, C.; Puhl, A. C. Déjà vu: stimulating open drug discovery for SARS-CoV-2. Drug. Discov. Today 2020, 5, 928−941.
    

42. [42]

    Sk, M. F.; Roy, R.; Jonniya, N. A.; Poddar, S.; Kar, P. Elucidating biophysical basis of binding of inhibitors to SARS-CoV-2 main protease by using molecular dynamics simulations and free energy calculations. J. Biomol. Struct. Dyn. 2020, 0739, 1−21.
    

43. [43]

    Yin, W. C.; Mao, C. Y.; Luan, X. D.; Shen, D. D.; Shen, Q. Y.; Su, H. X.; Wang, X. X.; Zhou, F. L.; Zhao, W. F.; Gao, M. Q.; Chang, S. H.; Xie, Y. C.; Tian, G. H.; Jiang, H. W.; Tao, S. C.; Shen, J. S.; Jiang, Y.; Jiang, H. L.; Xu, Y. C.; Zhang, S. Y.; Zhang, Y.; Xu, H. Y. Structural basis for inhibition of the RNA-dependent RNA polymerase from SARS-CoV-2 by Remdesivir. Science 2020, 6498, 1499−1504.