key: cord-288862-upcsvjuo
authors: Wang, Junmei
title: Fast Identification of Possible Drug Treatment of Coronavirus Disease-19 (COVID-19) through Computational Drug Repurposing Study
date: 2020-04-21
journal: J Chem Inf Model
DOI: 10.1021/acs.jcim.0c00179
sha: 
doc_id: 288862
cord_uid: upcsvjuo

[Image: see text] The recent outbreak of novel coronavirus disease-19 (COVID-19) calls for and welcomes possible treatment strategies using drugs on the market. It is very efficient to apply computer-aided drug design techniques to quickly identify promising drug repurposing candidates, especially after the detailed 3D structures of key viral proteins are resolved. The virus causing COVID-19 is SARS-CoV-2. Taking advantage of a recently released crystal structure of SARS-CoV-2 main protease in complex with a covalently bonded inhibitor, N3 (Liu et al., 10.2210/pdb6LU7/pdb), I conducted virtual docking screening of approved drugs and drug candidates in clinical trials. For the top docking hits, I then performed molecular dynamics simulations followed by binding free energy calculations using an end point method called MM-PBSA-WSAS (molecular mechanics/Poisson–Boltzmann surface area/weighted solvent-accessible surface area; Wang, Chem. Rev.2019, 119, 947831244000; Wang, Curr. Comput.-Aided Drug Des.2006, 2, 287; Wang; HouJ. Chem. Inf. Model., 2012, 52, 119922497310). Several promising known drugs stand out as potential inhibitors of SARS-CoV-2 main protease, including carfilzomib, eravacycline, valrubicin, lopinavir, and elbasvir. Carfilzomib, an approved anticancer drug acting as a proteasome inhibitor, has the best MM-PBSA-WSAS binding free energy, −13.8 kcal/mol. The second-best repurposing drug candidate, eravacycline, is synthetic halogenated tetracycline class antibiotic. Streptomycin, another antibiotic and a charged molecule, also demonstrates some inhibitory effect, even though the predicted binding free energy of the charged form (−3.8 kcal/mol) is not nearly as low as that of the neutral form (−7.9 kcal/mol). One bioactive, PubChem 23727975, has a binding free energy of −12.9 kcal/mol. Detailed receptor–ligand interactions were analyzed and hot spots for the receptor–ligand binding were identified. I found that one hot spot residue, His41, is a conserved residue across many viruses including SARS-CoV, SARS-CoV-2, MERS-CoV, and hepatitis C virus (HCV). The findings of this study can facilitate rational drug design targeting the SARS-CoV-2 main protease.

A great application of drug repurposing is to identify drugs that were developed for treating other diseases to treat a new disease. Drug repurposing can be achieved by conducting systematic drug−drug target interaction (DTI) and drug−drug interaction (DDI) analyses. We have conducted a survey on DTIs collected by the DrugBank database 5 and found that on average each drug has 3 drug targets and each drug target has 4.7 drugs. 6 The analysis demonstrates that polypharmacology is a common phenomenon. It is important to identify potential DTIs for both approved drugs and drug candidates, which serves as the basis of repurposing drugs and selection of drug targets without DTIs that may cause side effects. Polypharmacology opens novel avenues to rationally design the next generation of more effective but less toxic therapeutic agents. Computer-aided drug design (CADD) has been playing essential roles in modern drug discovery and development. To balance the computational efficiency and accuracy, a hierarchical strategy employing different types of scoring functions are applied in both the drug lead identification and optimization phases. A docking scoring function, such as the one employed by the Glide docking program, 7 is very efficient and thus can be utilized to screen a large library, but it is not very accurate. On the other hand, the molecular mechanical force field (MMFF)-based scoring functions are physical and more accurate but much less efficient. With the ever increasing computer power, MMFF-based free energy calculation methods, such as the end point MM-PBSA (molecular mechanics Poisson−Boltzmann surface area) and MM-GBSA (molecular mechanics generalized Born surface area) methods 2,3,8−21 and the alchemical thermodynamic integration (TI) and free energy perturbation (FEP) methods, 22, 23 have been extensively applied in structure-based drug discovery projects. We have developed a hierarchical virtual screening (HVS) to balance the efficiency and accuracy and improve the success rate of rational drug design. 8, 24 The newly released crystal structure of SARS-CoV-2 main protease 1 provides a solid structural basis for identification of drugs that might interact with this protein target. In this work, I applied multiscale modeling techniques to identify drugs that may be repurposed to target SARS-CoV-2 main protease. Flexible docking and MM-PBSA-weighted solvent-accessible surface area (WSAS) were applied as the first and second filters, respectively, to improve the efficiency and accuracy of HVS in inhibitor identification for SARS-CoV-2 main protease. Compared to the experimental means, CADD-based approaches are more efficient in providing possible treatment solutions for epidemic disease outbreaks like COVID-19. The detailed ligand−residue interaction profile, as well as the decomposition of binding free energy into different components, provides insight into rationally designing potent and selective inhibitors of SARS-CoV-2 main protease.

I conducted a hierarchical virtual screening (HVS) using the newly resolved crystal structure of SARS-CoV-2 main protease (resolution 2.16 Å). 1 Later on more crystal structures of SARS-CoV-2 main protease were resolved. 25 Two types of HVS filters were employed: Glide 7 flexible docking followed by MM-PBSA-WSAS. 2, 4 Detailed computational methods are described below.

2.1. Docking Screening. The crystal structure was first treated using the protein structure preparation wizard provided by the Schrodinger software, followed by docking grid generation. Glide flexible docking was performed using the default settings except that the formation of intramolecular hydrogen bonds was rewarded and the enhancement of planarity of conjugated π groups was turned on. The cocrystal ligand, N3, was covalently bonded to Cys145. I generated a new version of N3, N3′, by breaking the covalent bond and filling in open valence. I then evaluated whether Glide flexible docking can reproduce the native binding pose. In addition, a data set of approved drugs was prepared using DrugBank, 5 and a set of PubChem compounds that are structurally similar to lopinavir were enriched for docking screenings. Lopinavir, a potent inhibitor of HIV-1 protease, 26 was found to be effective in treating COVID-19 patients. Top hits from the docking screenings were advanced to the next HVS filter, MM-PBSA-WSAS. The 3D-structures of the screening compounds were generated using the OpenBabel software. 27 2.2. System Setup for Molecular Dynamics (MD) Simulation and Free Energy Calculation. MD simulations were first performed for a docking hit for two purposes: (1) studying the relative stability of the ligand residing in the binding pocket; (2) sampling a set of conformations for MM-PBSA-WSAS binding free energy calculations and MM-GBSA residue−ligand binding free energy decomposition analysis. A MD system consisted of one copy of SARS-CoV-2 main protease, one copy of docked ligand, 17597 TIP3P 28 water molecules, and about 50 Na + and Cl − ions depending on the charge state of the ligand. The whole system was neutralized. For the force field parameters, the partial atomic charges of ligands were derived using the RESP 29 program to fit the HF/6-31G* electrostatic potentials generated using the Gaussian 16 software package. 30 The other force field parameters of ligands come from GAFF, 31 and the AMBER FF14SB 32 force field was employed to model the viral protein. The residue topologies for ligands were prepared using the Antechamber module 33 implemented in the AMBER software package. 34 For the covalently bonded N3 ligand, I applied the residuegen program to generate nonstandard amino acid residue topology.

2.3. MD Simulation Protocols. For a protein−ligand complex, the MD system was first relaxed through a series of minimization procedures. The main chain atoms of the receptor and the bound ligand were restrained using a harmonic potential, and its force constant decreased from 20 to 10, 5, 1, and 0 (kcal/mol)/Å 2 , progressively in five 10 000-step minimizations. Note that the last step applied no restraint at all as the force constant is 0. The system was further relaxed by a set of 100 ps atomistic MD simulations with the same restraint setting of minimizations.

There were three phases for a MD simulation: the relaxation phase, the equilibrium phase, and the sampling phase. In the relaxation phase, the simulation system was heated progressively from 50 to 250 K in steps of 50 K. At each temperature, a 1 ns MD simulation was performed without any restraints or constraints. In the next equilibrium phase, the system was equilibrated at 298 K, 1 bar for 10 ns. Finally, a 100 ns MD simulation was performed at 298 K, 1 bar to produce NTP (constant temperature and pressure) ensembles. In total, 10 000 snapshots were recorded from the last simulation; 200 snapshots were evenly selected for the MM-PBSA-WSAS binding free energy calculation, and 5000 were selected for the MM-GBSA ligand−protein binding free energy decomposition analysis. Additional settings for constant pressure MD simulations performed in this work are listed as follows: temperature was regulated using Langevin dynamics 35 with a collision frequency of 5 ps −1 ; pressure was regulated using the isotropic position scaling algorithm with the pressure relaxation time being set to 1.0 ps; integration of the equations of motion was conducted at a time step of 1 fs for the relaxation phase and 2 fs for the equilibrium and sampling phases. The particle mesh Ewald (PME) method 36 was used to calculate the full electrostatic energy of a unit cell in a macroscopic lattice of repeating images. All bonds were constrained using the SHAKE algorithm 37 in both the minimization and MD simulation stages. All MD simulations were performed using the pmemd program in AMBER 2018. 34 2.4. MM-PBSA-WSAS Binding Free Energy Calculation. A total of 200 MD snapshots were evenly selected for the binding free energy calculations. For each selected MD snapshot, the molecular mechanical (MM) energy (E MM ) and the MM-PBSA solvation free energy were calculated without further minimization. 8,10,11,38−40 Key parameters controlling the MM-PBSA-WSAS analyses are listed as follows: external dielectric constant, 80; internal dielectric constant, 1; surface tension for estimating the nonpolar solvation energy by using solvent assessible surface area, 0.054. The Parse radii 41 were used in the MM-PBSA solvation calculation using the Delphi package (http://compbio.clemson.edu/delphi). The entropic term was estimated using a method coined WSAS (weighted solvent accessible surface area) described elsewhere. 4 It is noted that the entropic contribution cannot be neglected for this protein target as most ligands are large and have many rotatable bonds.

2.5. MM-GBSA Ligand−Residue Free Energy Decomposition Analysis. I conducted ligand−residue free energy decomposition analysis for 5000 snapshots evenly selected from the sampled snapshots. Besides the electrostatic and van der Waals interactions, the solvation effect was taken into account using a generalized GB model developed by Hawkins et al. 42 The ligand−residue MM-GBSA interaction energies were calculated using the Sander program in AMBER 2018. 34 Data analysis was performed using an internal program developed by us. A hot spot residue is recognized when its ligand−residue MM-GBSA interaction is stronger than −1.0 kcal/mol.

In this work, I have performed two-step hierarchical virtual screenings to identify drugs for repurposing to target SARS-CoV-2 main protease from a set of 2201 approved drugs downloaded from DrugBank. 5 In step 1, Glide docking, 7 an efficient but less accurate method, was applied to enrich repurposing drug candidates; in step 2, the docking hits were further evaluated using a more accurate but less efficient method, MM-PBSA-WSAS. The final repurposing drug candidates were selected based on the MM-PBSA-WSAS binding free energies.

3.1. Docking Screenings. I first evaluated the docking power of Glide program for the cocrystal ligand of SARS-CoV-2 main protease, N3. The ligand root mean square deviation (RMSD) of the best docking pose based on docking score (−9.4), 3.3 Å, was acceptable for a big ligand of about 100 atoms in flexible docking. I then applied the docking setting to conduct docking screenings. All the drug molecules that had docking scores better than −8.5 kcal/mol, roughly corresponding to 1 μM, were selected as hits and advanced to the next filter, MM-PBSA-WSAS. By utilizing this cutoff, the number of hits accounted for about 1% of total screening compounds.

For the promising docking hits, I conducted molecular dynamics simulations using the AMBER software package. 34 In total, 39 ligands including the cocrystal N3 ligand and 5 bioactives that are structurally similar to lopinavir were studied in the second phase of HVS. The top 5 approved neutral drugs that have excellent MM-PBSA-WSAS binding free energies (ΔG bind ≤ −5.0 kcal/mol) are shown in Figure 1 . The 2D structures of charged drugs with at least one form achieving ΔG bind ≤ −5.0 kcal/mol are shown in Figure 2 . I also found two bioactives ( Figure 3 ) that have excellent binding free energies (section 3.3). It is noted that lopinavir was observed to be effective in treating COVID-19 patients.

I explored the MD stability of each MD system. Figure 4 showed the RMSD fluctuations along the MD simulation time. It is shown that the main chain atoms of the receptor (black curves) and the secondary structures (red curves) reached equilibrium after 20 ns. The least-squares (LS) fitted RMSDs of the ligands (green curves) are around 2 Å, which is reasonable for large ligands like SARS-CoV-2 main protease inhibitors. A ligand's no-fit RMSD was calculated by first performing LSfitting for the main chain atoms of the receptor, and the resulting translation−rotation matrix was applied to the ligand, and then the RMSD was calculated directly. Evidently, the ligand no-fit RMSD measures not only the conformational changes but also its translational and rotational movements inside the binding pocket. The ligand no-fit RMSDs (blue curves) are larger than the LS-fit values; however, those values around 3−4 Å are still acceptable for large ligands. It is pointed out that sometimes the LS-fit and no-fit RMSD values are overestimated due to rotation of symmetrical functional groups of a ligand. In summary, the RMSD fluctuation analysis suggests that the MD trajectories are overall stable during the sampling phase for all the studied MD systems.

I first generated the average structure of the collected snapshots, and then the MD snapshot that had the smallest main chain atom RMSD against the average structure was chosen as the representative conformation. The comparisons between the crystal structure and the representative MD conformations are shown in Figure 5 for the native ligand N3 and in Figures 6−8 for the other ligands. As shown in Figure 5 , the benzene motif located in the dashed red cycle was inserted between two hot spot residues (His41 and Met49) for the MD structure, which is quite distinct from the crystal structure, which shows that the benzene motif has no direct interactions with His41 and Met49. I believe that under physiological conditions, the benzene motif becomes less solvent exposed and has more favorable interactions with His41 and Met49 by inserting itself between the side chains of the two residues. It is known that histidine has versatile roles in protein interactions. 43 After the sulfide bond is formed as a result of nucleophilic attack on catalytic Cys145 by the N3 ligand, His41 can stabilize the ligand−protein interaction by forming π−π stacking interactions between the His41 and the benzene motif of N3. As Table 1 . The calculated entropic term, TΔS, is quite different for different ligands as shown in Table 1 , suggesting the necessity of including this term in binding free energy calculations. The structures of the promising drug repurposing candidates that have both excellent docking scores and MM-PBSA-WSAS binding affinities are shown in Figures 1−3. All the known drugs shown in Figure 1 are neutral and have a better MM-PBSA-WSAS affinity than −5.0 kcal/mol. It should be noted that the cocrystal ligand, N3 is covalently bonded to the receptor; therefore its binding free energy is not directly comparable to those noncovalent ligands. The individual terms of MM-PBSA- Table S1 . The values of each energy term, van der Waals (ΔE VDW ), electrostatics (ΔE VDW + ΔG PB ), nonpolar solvation term (ΔG SA ), and entropy (TΔS), vary significantly from one system to another (Table 1 and Table S1 ), suggesting that there is no single energy term that dominates the protein−ligand interaction.

For the charged drug molecules, caution should be taken in result interpretation. For example, the neutral form of streptomycin (DB01082) has a MM-PBSA-WSAS binding free energy of −7.92 kcal/mol, much better than the charged form (−3.82 kcal/mol). The difference is caused by the distinct electrostatic properties between the neutral and charged molecules. However, the charged form is dominant under physiological conditions; 44 we therefore should use the result of the charged form or take the penalty of protonation into consideration when using the result of the neutral form.

Decomposition. I performed MM-GBSA binding free energy decomposition to identify the hot spot residues that make substantial contributions to the protein−ligand binding. The identified hot spots could enable us to rationally design potent and selective inhibitors of this drug target. To obtain statistically meaningful results, I studied 5000 MD snapshots for each system, and both the average ligand−residue interaction energies (ΔG lig−res ) and their RMSD values were calculated.

A hot spot residue is defined as a residue with ΔG lig−res equal to or smaller than −1.0 kcal/mol. The identified hot spots of each ligand are summarized in Table 2 . The most significant hotspot residues (ΔG lig−res < −3.0) are illustrated in Figures  5−8 . The common significant hot spot residues for most ligands (in bold in Table 2 ) are His41, Met49, Asn142, His164, Met165, Glu166, and Gln189. Asn Gly 

The outbreak of highly infectious diseases such as COVID-19 demands identification of multiple treatment plans as soon as possible. Computational drug repurposing studies can provide treatment options in a short period of time. For this study, amounts of computational time used for individual tasks are as follows. Docking screenings of all 2201 approved drugs with a single CPU core (Intel Xeon CPU E5-2683) took 11 h. For each docking hit, we need to perform ab initio calculations to derive point charges. The ab initio calculation using wB97XD/6-31G*//HF/6-31G* consumed about 1 day using four CPU cores; then it took us about 1.2 days to sample 120 ns using one GTX-1080 ti GPU; the following MM-PBSA-WSAS calculation consumed 1 day. Therefore, equipped with sufficient numbers of CPUs and GPUs and the current hardware, we can finish the drug repurposing screenings within 4 to 5 days using a reliable HVS strategy. Given that the inhibitors of SARS-CoV-2 main protease have relatively large sizes, the screening time can be even shorter for other drug targets with smaller ligands. Another consideration is the availability of high-quality drug target structures. Luckily, a high-resolution crystal structure of SARS-CoV-2 protease in complex with a ligand was resolved quickly, allowing us to conduct this drug repurposing screening. If no high-quality structure is available, one can rely on homology modeling techniques, probably with a reduced success rate of identifying repurposing drugs. Take SARS-CoV-2 main protease as an example: I performed structural alignments using an internal program that takes a multiplesequence-alignment (MSA) as an input. The MSA was generated by using the Promals3D web server. 45 The structure of SARS-CoV-2 main protease is found to be most similar to those of SARS-CoV protease (PDB 3TNT 46 ) and less similar to MERS-CoV protease (PDB 5WKK 47 ) ( Figure 9A ). In comparison, the structure of hepatitis C virus (HCV) NS3/4A (PDB code 3M5L 48 ) is quite different from the main proteases of coronaviruses: the RMSD of 2.26 Å between HCV and SARS-CoV-2 is much larger and with only 108 residues participating the least-squares fitting ( Figure 9B ). I also compared the sequences of the four proteases around the seven hot spot residues, which are colored in red in Table 3 . It is shown that SARS-CoV-2 and SARS-CoV share all seven hot spot residues. MERS-CoV and SARS-CoV-2 have four of the seven hot spot residues in common, while HCV NS3/4A and SARS-CoV-2 have only one hot spot residue (His41) in common. Even though the sequence identity is low between SARS-CoV-2 main protease and HCV NS3/4A, as shown in Figure 9B , the cocrystal ligands, N3 (green sticks) for SARS-CoV-2 and ITMN-191 (brown sticks) for HCV NS3/4A, largely overlap. This suggests that homology models can be constructed using Modeller 49 with SARS-CoV, MERS-CoV, and even HCV NS3/4A as templates.

In this work, we have conducted a computational drug repurposing study for the SARS-COV-2 main protease. To find robust treatments of COVID-19 particularly after the virus has developed different variations, it is necessary to screen repurposing drugs targeting other proteins that are essential in the life cycle of the virus.

In this study, I took advantage of the recently released crystal structure of SARS-CoV-2 main protease and conducted multiscale drug repurposing screenings. Five neutral drugs, namely, carfilzomib, eravacycline, valrubicin, lopinavir, and elbasvir, are identified to have inhibitory activities against SARS-CoV-2 main protease. Streptomycin, a charged molecule may also be an inhibitor of this SARS-CoV-2 main protease. Our study suggests that computational drug repurposing screening is very efficient and it can provide potential repurposing drug candidates in less than 5 days. A set of hot spot residues that make substantial contributions to the protein−ligand binding

Crystal structure of COVID-19 main protease in complex with an inhibitor N3

End-Point Binding Free Energy Calculation with MM/PBSA and MM/GBSA: Strategies and Applications in Drug Design

Recent Advances in Free Energy Calculations with a Combination of Molecular Mechanics and Continuum Models

Develop and test a solvent accessible surface area-based model in conformational entropy calculations

Development and Testing of Druglike Screening Libraries

Glide: a new approach for rapid, accurate docking and scoring. 1. Method and assessment of docking accuracy

Use of MM-PBSA in reproducing the binding free energies to HIV-1 RT of TIBO derivatives and predicting the binding mode to HIV-1 RT of efavirenz by docking and MM-PBSA

Revisiting free energy calculations: a theoretical connection to MM/PBSA and direct calculation of the association free energy

Validation and use of the MM-PBSA approach for drug discovery

Assessing the performance of the MM/PBSA and MM/GBSA methods. 1. The accuracy of binding free energy calculations based on molecular dynamics simulations

Assessing the performance of MM/PBSA and MM/GBSA methods. 3. The impact of force fields and ligand charge models

Improved docking performance using high solute dielectric constant MM/GBSA and MM/PBSA rescoring

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

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

Combined virtual screening, MMPBSA, molecular docking and dynamics studies against deadly anthrax: An in silico effort to inhibit Bacillus anthracis nucleoside hydrolase

Assessing the performance of MM/PBSA and MM/GBSA methods. 8. Predicting binding free energies and poses of protein-RNA complexes

Assessing the Performance of MM/PBSA, MM/GBSA, and QM-MM/GBSA Approaches on Protein/Carbohydrate Complexes: Effect of Implicit Solvent Models, QM Methods, and Entropic Contributions

Assessing the performance of MM/PBSA and MM/GBSA methods. 7. Entropy effects on the performance of endpoint binding free energy calculation approaches

Assessing the performance of the MM/PBSA and MM/ GBSA methods. 10. Impacts of enhanced sampling and variable dielectric model on protein-protein Interactions

Assessing the performance of MM/PBSA and MM/GBSA methods. 9. Prediction reliability of binding affinities and binding poses for proteinpeptide complexes

Toward Fast and Accurate Binding Affinity Prediction with pmemdGTI: An Efficient Implementation of GPU-Accelerated Thermodynamic Integration

Accurate and reliable prediction of relative ligand binding potency in prospective drug discovery by way of a modern free-energy calculation protocol and force field

Hierarchical database screenings for HIV-1 reverse transcriptase using a pharmacophore model, rigid docking, solvation docking, and MM-PB/SA

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

The HNRC Group. Lopinavir concentrations in cerebrospinal fluid exceed the 50% inhibitory concentration for HIV

Open Babel: An open chemical toolbox

Comparison of simple potential functions for simulating liquid water

A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges-the RESP model

Development and testing of a general amber force field

Improving the Accuracy of Protein Side Chain and Backbone Parameters from ff99SB

Automatic atom type and bond type perception in molecular mechanical calculations

Langevin stabilization of molecular-dynamics simulations of polymers by means of quasisymplectic algorithms

New tricks for modelers from the crystallography toolkit: the particle mesh Ewald algorithm and its use in nucleic acid simulations

Settle -an analytical version of the Shake and Rattle algorithm for rigid water models

Can MM-PBSA calculations predict the specificities of protein kinase inhibitors?

An improved method to predict the entropy term with the MM/PBSA approach

Assessing the performance of the molecular mechanics/Poisson Boltzmann surface area and molecular mechanics/generalized Born surface area methods. II. The accuracy of ranking poses generated from docking

Accurate Calculation of Hydration Free Energies Using Macroscopic Solvent Models

Parametrized models of aqueous free energies of solvation based on pairwise descreening of solute atomic charges from a dielectric medium

The multiple roles of histidine in protein interactions

Dihydrostreptomycin Directly Binds to, Modulates, and Passes through the MscL Channel Pore

PROMALS: towards accurate multiple sequence alignments of distantly related proteins

From SARS to MERS: crystallographic studies on coronaviral proteases enable antiviral drug design

Structure-guided design of potent and permeable inhibitors of MERS coronavirus 3CL protease that utilize a piperidine moiety as a novel design element

Drug resistance against HCV NS3/4A inhibitors is defined by the balance of substrate recognition versus inhibitor binding

Comparative Protein Structure Modeling Using MODELLER