Fast Identification of Possible Drug Treatment of Coronavirus Disease-19 (COVID-19) through Computational Drug Repurposing Study.

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, 9478; Wang, Curr. Comput.-Aided Drug Des. 2006, 2, 287; Wang; ; Hou J. Chem. Inf. Model., 2012, 52, 1199). 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.

Introduction

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.

Methodologies

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.

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]

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.

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]

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 (EMM) 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.

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.

Results

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.

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.

MD Simulations

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 (ΔGbind ≤ −5.0 kcal/mol) are shown in Figure . The 2D structures of charged drugs with at least one form achieving ΔGbind ≤ −5.0 kcal/mol are shown in Figure . I also found two bioactives (Figure ) that have excellent binding free energies (section ). It is noted that lopinavir was observed to be effective in treating COVID-19 patients.

Figure 1

Two-dimensional structures of promising drugs for repurposing. All five approved drugs are in neural form under physiological conditions.

Figure 2

Two-dimensional structures of promising drugs for repurposing. All three approved drugs are in charged form under physiological conditions.

Figure 3

Two-dimensional structures of promising bioactives that are structurally similar to lopinavir. PubChem 88143175, although studied in neutral form, bears a −3 charge under physiological conditions.

I explored the MD stability of each MD system. Figure 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 LS-fitting 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.

Figure 4

Plots of root-mean-square deviations of receptor main chain atoms and ligand heavy atoms along the MD simulation time for (A) cocrystal ligand N3, (B) DB08889, (C) DB12329, (D) DB00385, (E) DB01601, and (F) DB11574.

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 for the native ligand N3 and in Figures –8 for the other ligands. As shown in Figure , 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 shown below, both His41 and Met49 were hot spot residues in my MM-GBSA free energy decomposition analysis.

Figure 5

Structural comparison between the crystal structure and a representative MD structure of SARS-CoV-2 main protease bound to the known ligand N3. The crystal structure is shown as blue cartoon with the cocrystal ligand shown as brown sticks, while the representative MD structure is shown in gray cartoon and the ligand as green sticks (panel A). The hot spot residues (ΔGlig–res < −3.0 kcal/mol) revealed by MM-GBSA analysis are shown in panel B; the more bluish a residue is colored, the stronger the interaction between the residue and the ligand.

Figure 6

Structural comparison between the crystal structure and a representative MD structure of SARS-CoV-2 main protease bound to three neutral ligands, DB08889, DB12329, and DB00385. The crystal structure is shown as blue cartoon with the docked ligand shown as brown sticks, while the representative MD structure is shown in gray cartoon and the ligand as green sticks. (A) DB08889, (B) DB12329, and (C) DB00385. The detailed ligand–receptor interactions are shown in the panels D–F. All the hot spot residues (ΔGlig–res < −3.0) revealed by MM-GBSA analyses are labeled and colored by a blue to red spectrum; the more bluish a residue is colored, the stronger the interaction between the residue and the ligand. (D) DB08889, (E) DB12329, and (F) DB00385.

Figure 8

Structural comparison between the crystal structure and a representative MD structure of SARS-CoV-2 main protease bound to three charged ligands, DB01082, DB03147, and DB11184. The crystal structure is shown as blue cartoon with the docked ligand shown as brown sticks, while the representative MD structure is shown in gray cartoon and the ligand as green sticks. (A) DB01082, (B) DB03147, and (C) DB11184. The detailed ligand–receptor interactions are shown in panels D–F. All the hot spot residues (ΔGlig–res < −3.0) revealed by MM-GBSA analyses are labeled and colored by a blue to red spectrum; the more bluish a residue is colored, the stronger the interaction between the residue and the ligand. (D) DB01082, (E) DB03147, and (F) DB11184.

Structural comparison between the crystal and a representative MD structure of SARS-CoV-2 main protease bound to two neutral ligands, DB01601 and DB11574. The crystal structure is shown as blue cartoon with the docked ligand shown as brown sticks, while the representative MD structure is shown in gray cartoon and the ligand as green sticks. (A) DB01601 and (B) DB11574. The detailed ligand–receptor interactions are shown in panels C and D. All the hot spot residues (ΔGlig–res < −3.0) revealed by MM-GBSA analyses are labeled and colored by a blue to red spectrum; the more bluish a residue is colored, the stronger interaction between the residue and the ligand. (C) DB01601 and (D) DB11574.

MM-PBSA-WSAS Binding Free Energy Calculations

I measured the ligand binding affinity using the end point MM-PBSA method. Considering that the ligands of SARS-CoV-2 main protease are flexible molecules with large sizes, the contribution of the binding free energy from the conformational entropy change due to protein–ligand binding cannot be neglected. Instead of applying normal-mode analysis to estimate the entropic effect, I applied an efficient method called WSAS[4] to calculate this energy term. This scoring function is therefore called MM-PBSA-WSAS. The calculated binding free energies and the Glide docking scores are summarized in Table . The calculated entropic term, TΔS, is quite different for different ligands as shown in Table , 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 –3. All the known drugs shown in Figure 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-WSAS binding free energies of other less potent ligands are summarized in Table S1.

Table 1

List of Glide Docking Scores and MM-PBSA-WSAS Binding Free Energies (in kcal/mol) for Potential Inhibitors Binding to SARS-CoV-2 Main Proteasea

compound name	docking score	ΔEEEL	ΔEVDW	ΔGPB	ΔGSA	TΔS	ΔGbind
co-crystal ligand covalently bonds to sulfur of Cys145		–75.1 ± 0.3	–82.0 ± 0.4	84.0 ± 0.2	–6.5 ± 0.0	–29.8 ± 0.1	–38.8 ± 0.6
co-crystal ligand (no covalent bond formed)	–9.4	–52.3 ± 0.3	–26.9 ± 0.3	56.3 ± 0.2	–4.8 ± 0.0	–24.1 ± 0.1	–3.6 ± 0.4
Neutral Approved Drugs
DB08889	–8.6	–75.7 ± 0.4	–40.9 ± 0.1	78.3 ± 0.4	–6.0 ± 0.0	–30.4 ± 0.1	–13.8 ± 0.2
DB12329	–8.8	–45.7 ± 0.4	–25.5 ± 0.6	45.0 ± 0.3	–3.4 ± 0.0	–21.8 ± 0.0	–7.7 ± 0.5
DB00385	–9.2	–59.8 ± 0.2	–21.5 ± 0.2	52.8 ± 0.4	–4.6 ± 0.0	–25.9 ± 0.1	–7.2 ± 0.1
DB01601 (lopinavir)	–9.8	–52.5 ± 0.3	–20.1 ± 0.6	46.6 ± 0.6	–4.6 ± 0.0	–23.9 ± 0.0	–6.6 ± 0.3
DB11574	–9.9	–70.6 ± 0.4	–21.8 ± 0.4	65.6 ± 0.6	–6.4 ± 0.0	–26.6 ± 0.2	–6.5 ± 0.3
Charged Approved Drugs
DB01082 (NC = 0)	–8.6	–46.0 ± 0.4	–71.7 ± 1.0	89.4 ± 0.9	–3.8 ± 0.0	–24.2 ± 0.0	–7.9 ± 0.4
DB01082 (NC = 2)	–6.9	–33.4 ± 0.6	–280.0 ± 1.6	291.7 ± 1.3	–3.4 ± 0.0	–21.2 ± 0.1	–3.8 ± 0.5
DB03147 (NC = 0)	–10.2	–67.8 ± 0.3	–67.6 ± 0.9	106.6 ± 0.6	–5.0 ± 0.0	–26.3 ± 0.1	–7.5 ± 0.5
DB03147 (NC = −2)	–8.3	–52.9 ± 0.1	123.0 ± 0.3	–74.3 ± 0.5	–4.5 ± 0.0	–23.2 ± 0.1	14.6 ± 0.3
DB11184 (NC = 0)	–8.6	–57.4 ± 0.8	–52.2 ± 0.4	82.2 ± 0.3	–4.9 ± 0.0	–25.4 ± 0.1	–6.8 ± 0.3
DB11184 (NC = −4)	–7.4	–50.1 ± 0.4	175.4 ± 1.6	–148.7 ± 0.8	–4.4 ± 0.0	–22.9 ± 0.1	–5.0 ± 0.5
Bioactives Structurally Similar to Lopinavir
23727975	–8.8	–63.8 ± 0.1	–50.1 ± 0.7	77.91 ± 0.1	–5.2 ± 0.0	–28.3 ± 0.1	–12.9 ± 0.5
88143175 (NC = 0)	–10.0	–73.9 ± 0.2	–104.5 ± 1.3	148.9 ± 0.9	–6.6 ± 0.0	–30.5 ± 0.0	–5.6 ± 0.6

PubChem IDs are listed for the two bio-actives. The entropic contribution was estimated using T = 298.15 K.

The values of each energy term, van der Waals (ΔEVDW), electrostatics (ΔEVDW + ΔGPB), nonpolar solvation term (ΔGSA), and entropy (TΔS), vary significantly from one system to another (Table 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.

MM-GBSA Free Energy 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 (ΔGlig–res) and their RMSD values were calculated.

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

Table 2

Ligand–Residue MM-GBSA Interaction Energies (kcal/mol)a

			neutral approved drug					charged approved drug			bioactive
residue ID	residue type	co-crystal ligand	DB08889	DB12329	DB00385	DB01601	DB11574	DB01082	DB03147	DB11184	23727975	88143175
24	Thr	–0.1	–0.2	–0.0	–0.0	–0.1	–0.9	–0.0	–0.1	–0.3	–0.1	–1.4
25	Thr	–1.5	–1.5	–0.2	–0.6	–1.6	–4.0	–0.1	–1.6	–3.2	–0.6	–3.7
26	Thr	–0.6	–0.9	–0.2	–0.2	–0.3	–6.2	–0.1	–3.1	–4.5	–0.2	–4.3
27	Leu	–3.7	–1.2	–0.5	–1.0	–1.3	–2.2	–0.7	–2.6	–2.5	–0.2	–1.6
28	Asn	–4.8	–0.1	–0.0	–0.1	–0.1	–0.2	–0.1	–0.2	–0.1	–0.1	–0.1
39	Pro	–0.8	–0.2	–0.1	–0.1	–0.2	–0.2	–0.3	–0.3	–0.1	–0.1	–0.2
41	His	–6.7	–6.6	–2.9	–5.4	–3.6	–3.8	–5.4	–7.8	–4.2	–3.0	–3.0
44	Cys	–1.0	–0.1	–0.1	–0.1	–1.4	–0.7	–0.0	–2.5	–0.4	–0.1	–0.4
45	Thr	–0.6	–0.3	–0.2	–0.1	–1.1	–0.6	–0.0	–1.8	–0.3	–0.2	–0.5
46	Ser	–0.9	–1.8	–1.5	–0.6	–4.3	–1.4	–0.1	–3.5	–1.3	–1.8	–2.7
49	Met	–4.4	–2.4	–3.1	–2.1	–5.6	–4.2	–1.1	–5.5	–2.2	–2.2	–3.9
52	Pro	–0.8	–0.1	–0.1	–0.1	–0.0	–0.2	–0.0	–0.1	–0.0	–0.1	–0.0
54	Tyr	–0.6	–0.1	–0.4	–0.2	–0.1	–0.5	–0.2	–0.2	–0.1	–0.1	–0.0
119	Asn	–0.2	–0.0	–0.1	–0.1	–0.0	–1.4	–0.0	–0.4	–0.0	–0.0	–0.0
140	Phe	–0.2	–1.0	–0.0	–0.1	–0.0	–0.1	–0.8	–0.1	–0.2	–0.4	–1.8
141	Leu	–0.3	–1.2	–0.1	–0.2	–0.0	–0.1	–1.0	–0.0	–0.3	0.5	–2.6
142	Asn	–3.5	–4.4	–3.5	–3.7	–0.6	–2.9	–6.7	–4.6	–4.6	–5.5	–6.8
143	Gly	–2.2	–0.8	–1.0	–0.4	–0.4	–2.3	–1.9	–3.8	–1.5	–2.5	–2.1
144	Ser	–13.0	–1.4	–0.1	–0.2	–0.1	–0.3	–2.1	–1.2	–0.6	–5.3	–2.3
145	Cys	–66.3	–3.2	–0.7	–1.5	–1.6	–2.1	–4.0	–4.3	–1.5	–3.0	–3.8
146	Gly	–10.2	–0.0	–0.0	–0.0	–0.0	–0.0	–0.1	–0.0	–0.0	–0.1	–0.1
147	Ser	–0.6	–0.1	–0.0	–0.0	–0.0	–0.0	–0.1	–0.0	–0.0	–0.0	–0.1
163	His	–3.9	–1.6	–0.1	–0.4	–0.1	–0.2	–2.8	–0.4	–0.6	–1.2	–1.2
164	His	–3.9	–4.9	–1.4	–3.6	–1.8	–1.6	–5.3	–5.4	–1.3	–2.5	–2.5
165	Met	–5.5	–7.8	–4.2	–4.6	–3.5	–3.6	–5.3	–4.9	–4.7	–6.5	–4.4
166	Glu	–7.3	–6.7	–6.0	–6.8	–0.8	–1.8	–17.1	–6.4	–4.9	–14.4	–23.8
167	Leu	–1.7	–1.4	–1.3	–2.1	–0.6	–1.3	–1.9	–0.8	–1.3	–2.1	–1.2
168	Pro	–2.7	–1.2	–1.1	–2.5	–1.3	–2.7	–2.6	–0.4	–1.2	–2.2	–0.9
170	Gly	–0.0	–0.1	–0.0	–0.0	–0.0	–0.1	–0.5	–0.0	–0.0	–0.1	–0.2
172	His	–0.2	–1.2	–0.1	–0.2	–0.1	–0.1	–0.7	–0.1	–0.3	–1.1	–1.1
186	Val	–0.3	–0.6	–0.3	–0.2	–0.3	–0.4	–0.1	–0.6	–0.3	–0.3	–0.3
187	Asp	–1.5	–1.3	–2.2	–1.4	–3.7	–1.9	–2.8	–2.0	–1.6	–1.3	–1.2
188	Arg	–1.6	–2.0	–1.9	–1.7	–1.7	–2.3	–0.6	–1.8	–1.8	–1.8	–1.2
189	Gln	–5.2	–13.8	–5.9	–11.5	–8.1	–7.7	–2.3	–6.1	–8.8	–10.1	–6.9
190	Thr	–1.4	–1.6	–3.9	–2.4	–1.6	–3.8	–0.0	–1.1	–3.2	–1.0	–0.6
191	Ala	–1.6	–2.5	–0.2	–0.7	–1.5	–1.6	–0.0	–0.2	–0.8	–0.5	–0.1
192	Gln	–2.2	–3.5	–4.0	–2.3	–1.1	–1.8	–0.2	–1.2	–6.3	–1.7	–0.8

PubChem IDs are listed for the two bio-actives. The common hot spots for all the ligands (highlighted in bold) are His41, Met49, Asn142, His164, Met165, Glu166, and Gln189.

Discussion

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 multiple-sequence-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 A). 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 B). I also compared the sequences of the four proteases around the seven hot spot residues, which are colored in red in Table . 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 B, 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.

Figure 9

Structural comparison of proteases among three coronavirus viruses (SARS-CoV-2, SARS-CoV, and MERS-CoV) (A) and between SARS-CoV-2 and hepatitis C NS3/4A proteases (B). The SARS-CoV-2 main protease is colored in gray, and its ligands are shown as green sticks. The following are the color codes for the other proteases: SARS-CoV protease and its cocrystal ligand, brown; MERS-CoV protease and its cocrystal ligand, blue; HCV NS3/4A, blue; cocrystal ligand of HCV NS3/4A, brown. Backbone RMSD between SARS-CoV and SARS-CoV-2 is 0.4711 Å, with 284 residues participating in the least-squares fitting and 22 omitted, and the backbone RMSD between MERS-CoV and SARS-CoV-2 is 0.41 Å, but with 195 residues participating in the least-squares fitting and 104 omitted. In contrast, the backbone RMSD between SARS-CoV-2 main protease and HCV NS3/4A is 2.2632 Å, with 108 residues participating in the least-squares fitting and 43 omitted.

Table 3

Sequence Comparison around Hot Spot Residues for Proteases of Four Types of Viruses

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.

Conclusion

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 are also identified, which can facilitate rational design of novel selective inhibitors targeting SARS-CoV-2 main protease.