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

HE Qing-Xiu LI Gung-Ping GUO Hi-Qiong WANG Yu-Xun CHU Hn HU Yong SHEN Yn LIN Zhi-Hu, b③ WANG Yun-Qing, b③

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

HE Qing-Xiua②LI Guang-Pinga②GUO Hai-QiongaWANG Yu-XuanaCHU HanaHU YongaSHEN YanaLIN Zhi-Huaa, b③WANG Yuan-Qianga, b③

a(400054)b(400054)

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) gained tremendous attention due to its high infectivity and pathogenicity.The 3-chymotrypsin-like hydrolase protease (Mpro) of SARS-CoV-2 has been proven to be an important target for anti-SARS-CoV-2 activity.To better identify the drugs with potential in treating coronavirus disease 2019 (COVID-19) caused by SARS-CoV-2 and according to the crystal structure of Mpro, we conducted a virtual screening of FDA-approved drugs and chemical agents that have entered clinical trials.As a result, 9 drug candidates with therapeutic potential for the treatment of COVID-19 and with good docking scores were identified to target SARS-CoV-2.Consequently, molecular dynamics (MD) simulation was performed to explore the dynamic interactions between the predicted drugs and Mpro.The binding mode during MD simulation showed that hydrogen bonding and hydrophobic interactions played an important role in the binding processes.Based on the binding free energy calculated by using MM/PBSA, Lopiravir, an inhibitor of human immunodeficiency virus (HIV) protease, is under investigation for the treatment of COVID-19 in combination with ritionavir, and it might inhibit Mpro effectively.Moreover, Ombitasvir, an inhibitor for non-structural protein 5A of hepatitis C virus (HCV), has good inhibitory potency for Mpro.It is notable that the GS-6620 has a binding free energy, with respect to binding Mpro, comparable to that of ombitasvir.Our study suggests that ombitasvir and lopinavir are good drug candidates for the treatment of COVID-19, and that GS-6620 has good anti-SARS-CoV-2 activity.

severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), virtual screening, molecular docking, molecular dynamics (MD) simulation, binding free energy;

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, Raopublished 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 silivirtual 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 3Dstructure 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 minimi-zation”, “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 candidateswere 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 wasemployed to maintain constant temperature and the pressure ensembles.The pressure was set to 1?atm and controlled by the anisotropic (-,-,-) 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:

Δbind=ΔMM+ΔSOL=ΔMM+ΔGB+ΔSA

Where Δbindis the binding energy and ΔMMdenotes 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.ΔSOLis the solvation energy, which includes the polar solvation energy (ΔGB) calculated with the generalized born (GB) approximation model[36,37]and the non-polar part (Δ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 de-composition 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 screenthe 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.

Fig.1.Mpro structure and binding site

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 forthe 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 nonstruc-tural protein 5B inhibitor).

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

To be continued

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 energyanalysis.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.Theroot 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.

Fig.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 redandblue dashed lines, respectively)

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 byusing 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 importantto 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.

Based on the linear relationship between affinity and free binding energy, the affinity with Mpro/Ombitasvir might be accurately predicted by the Δbindenergy.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–1in 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.

Fig.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)

Table 2. Binding Free Energies (ΔGbind) Calculated by Using MM/PBSA (kcal×mol-1)

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 complexabout 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 distancesare 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 ?).

Fig.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)

Based on the stable conformation, the energy decompo-sition 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 hydrogenand hydrophobic bonds, so that Lopinavir had a good anti-SARS-COV-2 activity.

Fig.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)

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 MDsimulation, 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.

Fig.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)

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.71kcal×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–31.44kcal×mol-1).These binding energies indicate that GS-6620 has good potential anti-SARS-COV-2 activity.

Fig.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)

4 CONCLUSION

The pneumonia caused by SARS-CoV-2 is highly infec-tious andpathogenic.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.

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

(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.2020, 395, 507-513.

(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.2020, 395, 497-506.

(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..2020, 91, 264-266.

(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.2020, 63, 457-460.

(6) Gallagher, T.M.; Buchmeier, M.J.Coronavirus spike proteins in viral entry and pathogenesis.2001, 279, 371-374.

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

(8) Seah, I.; Su, X.; Lingam, G.Revisiting the dangers of the coronavirus in the ophthalmology practice.2020,1-3.

(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.2016, 24, 490-502.

(10) Wang, W.; Tang, J.; Wei, F.Updated understanding of the outbreak of 2019 novel coronavirus (2019-nCoV) in Wuhan, China.2020, 92, 441-447.

(12) 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.2020, 9, 623-623.

(13) Andersen, K.G.; Rambaut, A.; Lipkin, W.I.; Holmes, E.C.The proximal origin of SARS-CoV-2..2020, 26, 450-452.

(14) Jain, A.N.Scoring noncovalent protein-ligand interactions: a continuous differentiable function tuned to compute binding affinities..1996, 10, 427-440.

(15) 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..2015, 55, 572-588.

(16) Chen, J.Z.; Wang, J.; Xie, X.Q.GPCR structure-based virtual screening approach for CB2 antagonist search.2007, 47, 1626-1637.

(17) 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..2014, 154, 735-744.

(18) 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..2015, 8, 79-93.

(19) 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.2016, 18, 898-913.

(20) 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.2019, 24, 3784-15.

(21) Salomon-Ferrer, R.; Go?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..2013, 9, 3878-3888.

(22) 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..2012, 8, 1542-1555.

(23) Loncharich, R.J.; Brooks, B.R.; Pastor, R.W.Langevin dynamics of peptides: the frictional dependence of isomerization rates of N-acetylalanyl-N?-methylamide.1992, 32, 523-535.

(24) Izaguirre, J.A.; Catarello, D.P.; Wozniak, J.M.; Skeel, R.D.Langevin stabilization of molecular dynamics..2001, 114, 2090-2098.

(25) Darden, T.; York, D.; Pedersen, L.Particle mesh Ewald: an N?log(N) method for Ewald sums in large systems..1993, 98, 10089-10092.

(26) Essmann, U.; Perera, L.; Berkowitz, M.L.; Darden, T.; Lee, H.; Pedersen, L.G.A smooth particle mesh Ewald method..1995, 103, 8577-8593.

(27) 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..1977, 23, 327-341.

(28) 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..2019, 40, 1138-1156.

(29) Wang, J.; Hou, T.Develop and test a solvent accessible surface area-based model in conformational entropy calculations..2012, 52, 1199-1212.

(30) 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.1996, 100, 19824-19839.

(31) Wang, W.; Kollman, P.A.Free energy calculations on dimer stability of the HIV protease using molecular dynamics and a continuum solvent model.2000, 303, 567-582.

(32) 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..2018, 20, 14450-14460.

(33) 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..2016, 18, 22129-22139.

(34) 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..2014, 16, 22035-22045.

(35) 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..2014, 16, 16719-16729.

(36) Tsui, V.; Case, D.A.Theory and applications of the generalized Born solvation model in macromolecular simulations.2000, 56, 275-291.

(37) Bashford, D.; Case, D.A.Generalized born models of macromolecular solvation effects..2000, 51, 129-152.

(38) Sitkoff, D.; Sharp, K.A.; Honig, B.Accurate calculation of hydration free energies using macroscopic solvent models..1994, 98, 1978-1988.

(39) Still, W.C.; Tempczyk, A.; Hawley, R.C.; Hendrickson, T.Semianalytical treatment of solvation for molecular mechanics and dynamics.1990, 112, 6127-6129.

(40) Weiser, J.; Shenkin, P.S.; Still, W.C.Approximate atomic surfaces from linear combinations of pairwise overlaps (LCPO).1999, 20, 217-230.

(41) 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.2016, 56, 1152-1163.

(42) 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.2020, 5, 928-941.

(43) 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..2020, 0739, 1-21.

(44) 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.2020, 6498, 1499-1504.

28 August 2020;

29 October 2020

① This project was supported by the National Natural Science Foundation of China (31400667), Chongqing Municipal Education Commission Science and Technology Research Project (KJZD-K201801102), Chongqing Research Program of Basic Research and Frontier Technology (cstc2018jcyjAX0683), and Opening Foundation of State Key Laboratory of Silkworm Genome Biology (sklsgb1819-2).Computational support from the Information Center of Chongqing University of Technology

② These authors have equal contribution to the study

Wang Yuan-Qiang, E-mail: wangyqnn@cqut.edu.cn, Lin Zhi-Hua, E-mail: zhlin@cqut.edu.cn

10.14102/j.cnki.0254–5861.2011–2966