key: cord-277870-o79wph9r
authors: Han, Yanqiang; Wang, Zhilong; Ren, Jiahao; Wei, Zhiyun; Li, Jinjin
title: Potential inhibitors for the novel coronavirus (SARS-CoV-2)
date: 2020-09-18
journal: Brief Bioinform
DOI: 10.1093/bib/bbaa209
sha: 
doc_id: 277870
cord_uid: o79wph9r

The lack of a vaccine or any effective treatment for the aggressive novel coronavirus disease (COVID-19) has created a sense of urgency for the discovery of effective drugs. Several repurposing pharmaceutical candidates have been reported or envisaged to inhibit the emerging infections of severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2), but their binding sites, binding affinities and inhibitory mechanisms are still unavailable. In this study, we use the ligand-protein docking program and molecular dynamic simulation to ab initio investigate the binding mechanism and inhibitory ability of seven clinically approved drugs (Chloroquine, Hydroxychloroquine, Remdesivir, Ritonavir, Beclabuvir, Indinavir and Favipiravir) and a recently designed α-ketoamide inhibitor (13b) at the molecular level. The results suggest that Chloroquine has the strongest binding affinity with 3CL hydrolase (Mpro) among clinically approved drugs, indicating its effective inhibitory ability for SARS-CoV-2. However, the newly designed inhibitor 13b shows potentially improved inhibition efficiency with larger binding energy compared with Chloroquine. We further calculate the important binding site residues at the active site and demonstrate that the MET 165 and HIE 163 contribute the most for 13b, while the MET 165 and GLN 189 for Chloroquine, based on residual energy decomposition analysis. The proposed work offers a higher research priority for 13b to treat the infection of SARS-CoV-2 and provides theoretical basis for further design of effective drug molecules with stronger inhibition.

The recent outbreak of the novel coronavirus disease caused by severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a serious global risk for health and economy. Unfortunately, no drug or vaccine has yet been approved to prevent the infection or to treat the patients, and more time is required to develop more specific strategy against this pathogen [1, 2] . Given the urgency of the situation, researchers have been trying a variety of clinically approved drugs in the hope that Chloroquine and Hydroxychloroquine are being reported to be potential inhibitors of SARS-CoV-2 in vitro [13, 14] . However, Remdesivir failed in a clinical trial against Ebola virus in 2019 [15] , and there has been no clinical trial supporting its efficacy against SARS-CoV-2. Empirical trials often involve substantial trial-and-error costs, which require inevitable time consuming in drug development. To prevent the spread of COVID-19 timely, it is necessary to apply all approaches, including computational analysis, to prioritize the pharmaceutical candidates, so as to rapidly develop a strategy for COVID-19 treatment.

The SARS-CoV-2 genome is comprised of approximate 30 000 nucleotides, where two overlapping polyprotein genes (pp1a and pp1ab) encode replicases essential for viral replications and transcriptions. The 3C-like main protease (3CL Mpro) is a main protease of SARS-CoV-2, essential for the immune regulation and the releases of the polyproteins pp1a and pp1ab [16] . 3CL Mpro controls the functional polypeptides releases from the polyproteins by extensive proteolytic processing [17] . Its functional importance in the viral life cycle and the absence of closely related homologs in human make it an attractive target for anti-SARS-CoV-2 drugs. Therefore, the discovery or design of drug molecules that can bind to 3CL Mpro and inhibit its function is an effective way to combat the virus. In this study, we chose 3CL Mpro as the therapeutic target to ab initio investigate its inhibition mechanism and binding ability of these most promising drug molecules by ligand-protein docking program (Rosetta) and molecular dynamics (MD) simulations. The experiments and high-throughput screening researches have been going on around the clock in hope of finding a cure for the COVID-19. Based on pioneers' studies, we selected seven most promising drugs (i.e. Chloroquine, Hydroxychloroquine, Remdesivir, Ritonavir, Beclabuvir, Indinavir and Favipiravir) and conducted highprecision quantum calculations, so as to find the most effective inhibitor. Recently, Zhang et al. [18] pursued a structure-based design of peptidomimetic α-ketoamides as inhibitors of 3CL Mpro, and demonstrated that the inhibitor 11r exhibited activity against coronavirus. Soon later, based on 11r, Zhang et al. [19] reported three structure modified inhibitors (i.e. 13a, 13b and 14b), and demonstrated that 13b is the most effective designed inhibitor, which showed good tropism to the lungs with a welltolerated administration through inhalation in mice. Therefore, in addition to the seven clinically approved drug molecules, we also included 13b, as a candidate to investigate its binding ability with 3CL Mpro.

The 3CL Mpro structure (PDB ID: 6 LU7) was obtained from the RCSB Protein Data Bank, which was recently released by Rao's group [20, 21] , while the 3D structures of the seven drug molecules were obtained from the PubChem database. Through molecular docking and kinetic analysis, we found that for repurposed drugs, the Chloroquine molecule has the strongest interaction with the 3CL Mpro, indicating that Chloroquine is the best potential inhibitor for SARS-CoV-2, followed by Hydroxychloroquine, Remdesivir, Ritonavir, Beclabuvir, Indinavir and Favipiravir. Although Chloroquine and Hydroxychloroquine share similar chemical structures and mechanisms of acting as a weak base and immunomodulator, we demonstrated that Chloroquine molecule has a stronger binding with 3CL Mpro than Hydroxychloroquine, based on the high-precision calculations of enthalpy and entropy. However, the newly designed inhibitor 13b represents slightly stronger binding affinity with 3CL Mpro compared with Chloroquine, suggesting that 13b may be more effective than the selected clinically approved drugs. We further analyzed the specific residues of 3CL Mpro that contribute significantly to the binding energy, and demonstrated that the residues MET165 and GLN189 in 3CL Mpro have the strongest affinities with Chloroquine molecule, with decomposition energies of −2.50 and −2.34 kcal/mol, respectively. For 13b, the residues MET165 and HIE163 in 3CL Mpro contribute most to the binding energy, with decomposition energies of −3.23 and −2.65 kcal/mol, respectively. The determination of dominant residues in 3CL Mpro provides a theoretical guidance for further drug molecule design with larger binding energies and stronger inhibitory abilities with target protein.

Since the complex structures between 3CL Mpro and inhibitors have not been identified, we used the ligand-protein docking program (Rosetta3) to determinate their 3D binding conformation. Based on the total energy of complex and interaction energy between ligand and protein calculated by Rosetta, we generated 10 000 conformations for each inhibitor and selected the most likely combination based on Rosetta score function. The total Rosetta score, the binding energy between protein and ligand and whether the ligand is touching the protein are considered for the conformational selection. Table S8 in Supplementary Information shows the total Rosetta energy score and important energy terms of the selected conformations for the seven clinically approved drugs from docking studies. The results show that the selected conformation of Hydroxychloroquine has the lowest Rosetta score, while Indinavir has the lowest binding energy. The docking binding energies of all seven clinically approved drugs with 3CL Mpro are lower than −10 kcal/mol, with the total energy score lower than −642 kcal/mol. The virtual docking results show that all seven inhibitors have high binding affinities to 3CL Mpro.

To investigate the binding ability and the interactions between inhibitors and 3CL Mpro, the most favorable complex conformations were used for molecular dynamic simulations. After the equilibrium run of 1.2 ns, the ligand-protein complex system became stable with the convergence of temperature, density and total energy. Then a 40 ns production run was performed while keeping the system stable. The root-meansquare deviation (RMSD) was analyzed to detect the stability of the ligand-protein system during simulation. Figure 1a shows the RMSD of all non-hydrogen atoms for the complex of 3CL Mpro and 13b over the entire production simulation. The result indicates that the ligand-protein complex mostly remains stable from 15 to 40 ns at 2.5-3 Å. Besides, the root means square fluctuation (RMSF) was analyzed to detect the fluctuation and flexibility of protein residues during the MD simulation. Figure 1b shows the RMSFs of all non-hydrogen atoms of 3CL Mpro over the entire simulation, where the LEU50, PRO52, ASN51, ASN53, MET49, GLN189, ARG222, ASP48, SER46, ASN277, GLU55 and TYR154 are the most flexible residues in the protein middle region while 13b binds to 3CL Mpro. The RMSDs and RMSFs of seven clinically approved drug molecules are plotted in Figures S3, S6 , S9, S12, S15, S18 and S21 of supporting information, respectively. From these figures, the RMSFs of the 3CL Mpro are similar when the eight molecules (Chloroquine, Hydroxychloroquine, Remdesivir, Ritonavir, Beclabuvir, Indinavir and Favipiravir and 13b) bind to the protein pocket. Table 1 shows the binding free energies of eight potential inhibitors to 3CL Mpro, where the total enthalpies (ΔE TOT ) were calculated by MD simulation and the total binding free energy (ΔG bind ) were obtained by the Molecular Mechanics-Generalized Born Surface Area (MM/GBSA)-binding free energy. The total enthalpies include the components of gas-phase and solvation energies but without the entropy, which sometimes plays an important role for ligand-protein interaction. So far, many effective methods have been applied to calculate the binding free energy of proteins and ligands. MM/GBSA is a time-saving and highly accurate method, which is often used to compute the absolute binding free energy and shows great success in elucidating the protein-protein and protein-inhibitor interactions [22] [23] [24] [25] . Based on the MM-GBSA method, we calculated the binding energy G bind between inhibitors and 3CL Mpro, consisting of van der Waals energy, electrostatic energy, polar solvation energy nonpolar solvation energy and entropic contribution. From Table 1 , the calculated total enthalpy E TOT by MD simulation (the green column) shows that for clinically approved drugs, Remdesivir has the strongest binding energy as −39.33 kcal/mol. However, the total enthalpy does not consider the contribution of entropy, which depends on the alternation of motional freedom induced by inhibitor binding and contributes partly to the total binding energy [26] . After taking into account of entropy contribution, Table 1 shows that Chloroquine exhibits the lowest binding free energy (−16.119 kcal/mol) among the seven clinically approved drugs (orange column), followed by Hydroxychloroquine (−13.7215 kcal/mol), Remdesivir (−10.876 kcal/mol), Ritonavir (−4.78 kcal/mol), Beclabuvir (−7.75 kcal/mol), Indinanir (−6.52 kcal/mol) and Favipiravir (4.21 kcal/mol). This finding is consistent with the recent in vitro study [13, 14] showing that Chloroquine, Hydroxychloroquine and Remdesivir are the most effective in controlling SARS-CoV-2 infection among a number of clinically approved drugs. For patients with acute porphyria or porphyria cutanea tarda, Chloroquine has some side effects of causing fever, serum aminotransferase elevations and even jaundice. Hydroxychloroquine, a less toxic metabolite of Chloroquine, does not cause these adverse side effects and was found to have partial beneficial effects in porphyria [27] . However, we demonstrated that Chloroquine has a stronger affinity and inhibitory ability to 3CL Mpro than Hydroxychloroquine. In addition, from Table 1 , the enthalpy and binding free energy of the newly proposed inhibitor 13b are higher than those of all selected clinically approved drug molecules, indicating the strongest binding affinity of 13b with 3CL Mpro. Figure 2a shows the complex structure of the binding between 3CL Mpro and 13b, with similar binding position to that of the analyzed clinically approved drugs. The complex structure of 3CL Mpro bounded by clinically approved drugs with the highest Rosetta score are shown in the Supplementary Information. According to the calculated results, 13b molecule as well as the clinically approved drugs binds tightly to 3CL Mpro, with the binding position close to that of N3 inhibitor proposed by Rao group [28] . Such binding pocket is also similar to the previous report for SARS 3CL Mpro with inhibitors [29] . Figure 2b shows the ball-and-stick model for the binding pocket of five residues with dominant binding contributions of 3CL Figure S17 and Table 2 ) very different from those of Chloroquine (MET 165, GLN 189, MET 49, ASP 187 and ARG 188 as shown in Figure S20 and Table 2 ). To reasonably design drugs and select appropriate ligands, we further analyzed the contribution of each residue of 3CL Mpro when binding with ligands. Figure 3a shows the top 15 residues with dominant binding attractions for 13b, while Figure 3b shows the contributions of the binding free energy components for the top five residues, which are decomposed into van der Waals energy (vdW), electrostatic energy (Ele), polar solvation energy (polar) and nonpolar solvation energy (nonpolar). In Figure 3a , the major binding attractions of 3CL Mpro with 13b molecule come from the residues MET 165, HIE 163, HIE 164, CYS 145 and ASN 142, which contribute more than −1.5 kcal/mol to the binding energy. The decomposed energy distribution in Figure 3b shows that the dominant favorable interactions of the above-mentioned residues come from the electrostatic energy (the purple bars). For Chloroquine (Figure S20 ), the top five residues (MET 165, GLN 189, MET 49, ASP 187 and ARG 188) also contribute more than −1.5 kcal/mol and the dominant favorable interactions of these residues are van der Waals interactions.

Furthermore, MET 165 also forms strong hydrophobic interaction with 13b and Chloroquine molecules, but not with Hydroxychloroquine ( Figure S17 ), which makes the former two bind stronger to Mpro than the latter. When 13b binding with 3CL Mpro, the electrostatic interaction contributes the most for residues HIE 163 and HIE 164, while the van der Waals interaction is beneficial to residues MET 165, CYS 145 and ASN 142, as shown in Figure 3b . In summary, the configuration and environment of ligands impact significantly on the interaction of ligand and residue pairs at the binding site, which leads to the binding free energy difference. In this work, residues MET 165 and HIE 163 in 3CL Mpro make significant contributions when binding with α-ketoamide inhibitor 13b, in which the van der Waals interaction and electrostatic interaction play the most important roles. The discussions of binding free energy decomposing of key residues for Chloroquine, Hydroxychloroquine, Remdesivir, Ritonavir, Beclabuvir, Indinavir and Favipiravir are shown in Supplementary Information. Such results can be used for inhibitor molecule design with stronger binding affinity with 3CL Mpro of SARS-CoV-2.

In conclusion, based on the molecular dynamic simulation, we ab initio determine that for clinically approved drugs, Chloroquine, a widely used antimalarial and autoimmune disease drug possesses the strongest binding affinity with 3CL Mpro. However, the newly designed inhibitor 13b exhibits stronger binding ability with 3CL Mpro, thus it is more likely to become a pharmaceutical inhibitor against SARS-CoV-2. The proposed work suggests that the α-ketoamide inhibitor 13b is expected to be a promising candidate to treat the infection of SARS-CoV-2 and provides theoretical guidance for further design of drug molecular structure with stronger inhibitory effect with target protein.

In this study, 3CL Mpro was selected as the target of docking research and binding analysis, due to its functional importance in the life cycle of SARS-CoV-2. The 3D crystal structure of 3CL Mpro (PDB ID: 6 LU7) was taken from RCSB Protein Data Bank, which was resolved by Rao's group [20] with a resolution of 2.16 Å.

For docking and binding analysis, seven clinically approved inhibitors (Chloroquine, Hydroxychloroquine, Remdesivir, Ritonavir, Beclabuvir, Indinavir and Favipiravir) were selected from previous reported virtual screening works or were found to be highly effective in the control of SARS-CoV-2 infection in vitro [6, 12] . And a newly designed α-ketoamide inhibitor (13b) was selected for comparative analysis, which was a co-crystallized ligand of the protein-ligand complex (PDB ID: 6Y2F) [19] .

In this study, we applied ligand-protein docking by Rosetta program [30] [31] [32] [33] , which attempted to find the conformation and relative orientation of each ligand to minimize the Rosetta score function. To explore the binding pocket, the initial position of ligand was set near the protein target. After Monte Carlo minimization, we obtained 1000 conformations for 3CL Mpro and ligand and selected the 5% candidates with the lowest energies according to their energy scores from the ligand-protein contact structures. The Rosetta total energy score is the weighted sum of energy terms, including physical forces (e.g. electrostatics), and statistical terms (e.g. possibility of finding torsion angles in the Ramachandran space). Then, according to the sorting of binding energy of ligand and protein, the structure with the lowest binding energy was selected for MD simulation. The presented binding energy calculation included van der Waals energy, electrostatic energy, polar solvation energy nonpolar solvation energy and entropic contribution; therefore, it considered the core differences between the ligand and protein target and eliminated some irrelevant noise in the protein variation.

Amber16 software was used to perform MD simulation and binding energy calculation between drug molecules and 3CL Mpro [34] . The AM1-BCC charges of drug molecules were computed by sqm program in amber16. The ff14SB force field [35] and General Amber force field (GAFF) [36] were used for parameters' generation of proteins and drug molecules, respectively. The ligand-protein complexes were immersed in a rectangular TIP3P water box (12A), with counterions (Na + and Cl − ) to neutralize the systems. The energy minimization and MD simulation were performed with sander program in amber16, during which a cutoff for nonbonded interactions of 8 Å was used and the time step was set to 2 fs. After 30 000 steps of energy minimization, the solvated ligand-protein complex was gradually heated from 0 to 300 K in the NVT ensemble using Langevin thermostat during 100 ps, followed by 100 ps of density equilibration with weak restraints on the complex (10 kcal mol −1 Å −2 ). Then 1 ns equilibration run was performed in NPT ensemble with temperature of 300 K and pressure pf 1 atm. At last, we performed a 40 ns production run in the NPT ensemble for each solvated complex system with a collection interval of 100 ps, and 300 frames were collected for further energy calculation. All hydrogen atoms were constrained with the SHAKE algorithm.

The binding free energy for small drug molecules and protein was calculated by the Molecular Mechanics-Generalized Born Surface Area (MM-GBSA) [22] [23] [24] [25] method using MMPBSA.py script based on the MD trajectory. The MM-GBSA free energy G bind was calculated by the following equation:

where E TOT is the total enthalpy from the generalized Born model, including internal and solvation energies, and S nmode is the entropy calculated by normal mode [37, 38] . Here, 150 frames were used for E TOT calculation while 50 frames were used for S nmode calculation. The enthalpy E TOT can be calculated by:

E TOT = E TOT−gas + E TOT−solv = E vdW + E ele + G P + G np (2) where E TOT−gas is the gas phase free energy, E TOT−solv is the solvation free energy, E vdW is the van der Waals energy, E ele is the electrostatic energy, G P and G np are the polar and nonpolar solvation free energies, respectively. The binding free energy G bind was determined by:

where G com , G lig and G rec were MM/GBSA free energies, corresponding to complex, ligand and protein target, respectively.

• Ab initio simulation is used to investigate the inhibitory ability between eight potential clinical pharmaceutical molecules and 3CL hydrolase in SARS-CoV-2. • Theoretical calculation demonstrates that the αketoamide 13b is the most promising inhibitor for SARS-CoV-2, while Chloroquine is the best among the clinically approved drugs.

• The MET 165 and HIE 163 contribute the most based on residual energy decomposition analysis.

Supplementary data are available online at Briefings in Bioinformatics Journal..

The data that support the findings of this study are available from the corresponding author, Professor Jinjin Li (email: lijinjin@sjtu.edu.cn) upon reasonable request.

From SARS to MERS, thrusting coronaviruses into the spotlight

Genomic and protein structure modelling analysis depicts the origin and infectivity of 2019-nCoV

Design of efficient computational workflows for in silico drug repurposing

Computational drug repurposing: current trends

Potential 2019-nCoV 3C-like protease inhibitors designed using generative deep learning approaches

Potential therapeutic agents for COVID-19 based on the analysis of protease and RNA polymerase docking

Drug repurposing for coronavirus (COVID-19): in silico screening of known drugs against coronavirus 3CL hydrolase and protease enzymes

Comparative therapeutic efficacy of remdesivir and combination lopinavir, ritonavir, and interferon beta against MERS-CoV

Broad-spectrum antiviral GS-5734 inhibits both epidemic and zoonotic coronaviruses

Coronavirus susceptibility to the antiviral remdesivir (GS-5734) is mediated by the viral polymerase and the proofreading exoribonuclease

Prophylactic and therapeutic remdesivir (GS-5734) treatment in the rhesus macaque model of MERS-CoV infection

First case of 2019 novel coronavirus in the United States

Hydroxychloroquine, a less toxic derivative of chloroquine, is effective in inhibiting SARS-CoV-2 infection in vitro

Remdesivir and chloroquine effectively inhibit the recently emerged novel coronavirus (2019-nCoV) in vitro

New treatments for Ebola virus disease

A pneumonia outbreak associated with a new coronavirus of probable bat origin

Conservation of substrate specificities among coronavirus main proteases

α-Ketoamides as broadspectrum inhibitors of coronavirus and enterovirus replication: structure-based design, synthesis, and activity assessment

Crystal structure of SARS-CoV-2 main protease provides a basis for design of improved αketoamide inhibitors

Structure of Mpro from COVID-19 virus and discovery of its inhibitors

The protein data bank

Semianalytical treatment of solvation for molecular mechanics and dynamics

Py: an efficient program for end-state free energy calculations

Exploring protein native states and large-scale conformational changes with a modified generalized born model

Performance comparison of generalized born and Poisson methods in the calculation of electrostatic solvation energies for protein structures

A comparative insight into amprenavir resistance of mutations V32I, G48V, I50V, I54V, and I84V in HIV-1 protease based on thermodynamic integration and MM-PBSA methods

Clinical and Research Information on Drug-Induced Liver Injury. Bethesda (MD): National Institute of Diabetes and Digestive and Kidney Diseases

The crystal structure of 2019-NCoV main protease in complex with an inhibitor N3. RCSB Protein Data Bank

Molecular dynamic simulations analysis of ritronavir and lopinavir as SARS-CoV 3CLpro inhibitors

Blind docking of pharmaceutically relevant compounds using RosettaLigand

RosettaLigand docking with full ligand and receptor flexibility

Status of GPCR modeling and docking as reflected by community-wide GPCR dock 2010 assessment

Computational characterization of the glutamate receptor antagonist perampanel and its close analogs: density functional exploration of conformational space and molecular docking study

ff14SB: improving the accuracy of protein side chain and backbone parameters from ff99SB

Development and testing of a general amber force field

Fast and accurate computation schemes for evaluating vibrational entropy of proteins

The normal-mode entropy in the MM/GBSA method: effect of system truncation, buffer region, and dielectric constant

Research reported in this publication was partially supported by the National Natural Science Foundation of China (Nos 51672176 and 21901157), and the SJTU Global Strategic Partnership Fund (2020 SJTU-HUJI).

The authors declare no competing interests.