Repurposing of FDA-approved antivirals, antibiotics, anthelmintics, antioxidants, and cell protectives against SARS-CoV-2 papain-like protease.

SARS-CoV-2 or Coronavirus disease 19 (COVID-19) is a rapidly spreading, highly contagious, and sometimes fatal disease for which drug discovery and vaccine development are critical. SARS-CoV-2 papain-like protease (PLpro) was used to virtually screen 1697 clinical FDA-approved drugs. Among the top results expected to bind with SARS-CoV-2 PLpro strongly were three cell protectives and antioxidants (NAD+, quercitrin, and oxiglutatione), three antivirals (ritonavir, moroxydine, and zanamivir), two antimicrobials (doripenem and sulfaguanidine), two anticancer drugs, three benzimidazole anthelmintics, one antacid (famotidine), three anti-hypertensive ACE receptor blockers (candesartan, losartan, and valsartan) and other miscellaneous systemically or topically acting drugs. The binding patterns of these drugs were superior to the previously identified SARS CoV PLpro inhibitor, 6-mercaptopurine (6-MP), suggesting a potential for repurposing these drugs to treat COVID-19. The objective of drug repurposing is the rapid relocation of safe and approved drugs by bypassing the lengthy pharmacokinetic, toxicity, and preclinical phases. The ten drugs with the highest estimated docking scores with favorable pharmacokinetics were subjected to molecular dynamics (MD) simulations followed by molecular mechanics/generalized Born surface area (MM/GBSA) binding energy calculations. Phenformin, quercetin, and ritonavir all demonstrated prospective binding affinities for COVID-19 PLpro over 50 ns MD simulations, with binding energy values of -56.6, -40.9, and -37.6 kcal/mol, respectively. Energetic and structural analyses showed phenformin was more stable than quercetin and ritonavir. The list of the drugs provided herein constitutes a primer for clinical application in COVID-19 patients and guidance for further antiviral studies.Communicated by Ramaswamy H. Sarma.

Introduction

SARS-CoV-2 or Coronavirus disease-19 (COVID-19) has recently emerged in China and rapidly spread throughout many other countries (WHO, 2020). The new SARS-CoV-2 was classified as a Beta-coronavirus (Liu et al., 2020). As recently evolved, there is a great demand for the discovery of new drugs against SARS-CoV-2.

The new SARS-CoV-2 is classified as a beta-coronavirus (Liu et al., 2020), characterized by a large genome of approximately 30 kb, and comprised of genes for four structural proteins—spike (S), envelop (E), membrane (M), and nucleocapsid (N)—and non-structural (NS) genes (Shi et al., 2015). In CoVs, there are 14–16 NS genes encoding proteins used in the vital processes of virus genomic transcription, protein processing, virion assembly, and replication (Brian & Baric, 2005). CoV infections caused mild respiratory symptoms until more dangerous forms evolved approximately two decades ago, noted by the emergence of SARS CoV. Three highly fatal CoVs have been of global health concern: SARS CoV, MERS CoV, and SARS-CoV-2 (Guarner, 2020). The SARS-CoV-2 genome composition is highly similar to that of SARS CoV, yet the former possesses inherently more efficient gene expression due to more biased codon usage and lower effective codon numbers than SARS CoV, pangolin Beta-CoV, bat SARS CoV, and MERS CoV (Kandeel et al., 2020; Kandeel & Altaher, 2017).

The CoV NS polyprotein is processed by two proteases, the main protease (Mpro) (a strict protease) and the multifunctional PLpro (Harcourt et al., 2004). Inhibition of PLpro impacts virus replication via deficient viral protein processing and can also affect remote PLpro activities, including deubiquitination, de-ISGylation, and innate, anti-host immunity reactions (Báez-Santos et al., 2014; Mielech et al., 2014). PLpro comprises a central catalytic domain and the ubiquitin-like domain (Ubl) (Bailey-Elkin et al., 2014; Lei et al., 2014; Yang et al., 2014), which interferes with the host IRF3 and NF-kB antiviral pathways and modulates the host immune response (Ratia et al., 2014). PLpro subdomains showed different reaction patterns in the presence or absence of ubiquitin (Alfuwaires et al., 2017b). Furthermore, PLpro was an important drug target against several CoVs, including MERS CoV (Kandeel et al., 2019) and SARS CoV (Báez-Santos et al., 2015).

FDA-approved drugs are the first potential treatments to use in emergent cases, due to their commercial availability, proven safety, and direct applicability without the need for any preclinical and clinical testing. This study aims to identify the top candidates among FDA-approved drugs with strong binding potential for SARS-CoV-2 PLpro, with the ultimate goal of repurposing drugs for inhibition of SARS-CoV-2. The most promising drugs based on the predicted docking scores for COVID-19 PLpro will be further investigated using molecular dynamics simulations over 50 ns. A molecular mechanics/generalized Born surface area (MM/GBSA) approach will be applied to evaluate drug- PLpro binding energies (ΔGbinding) throughout the MD simulated time. The stability, binding affinity, and interactions of the drug-PLpro binding will also be estimated.

Methods

Construction of FDA-approved drugs dataset

A dataset of FDA-approved drugs was obtained as an SDF file from Selleckchem Inc. (WA, USA). The drugs dataset was imported to the ligprep module of the Maestro package (Schrodinger Inc., LLC, NY, USA). The structures of drugs were 3 D optimized at neutral pH and desalted.

Preparation of SARS-CoV-2 PLpro structure

The structure of SARS-CoV-2 PLpro was retrieved from the protein data bank (PDB ID 6w9c) and imported to the protein preparation module in the Maestro software package. Protein preparation was as previously described (Alfuwaires et al., 2017a; Kandeel et al., 2020; Kandeel & Al-Nazawi, 2020). Briefly, the protein was preprocessed to assign the substructure components and detect potential errors. The structure was protonated and optimized at neutral pH. The protein structure was energy minimized using an OPLS2005 force field. The docking grid was generated by selecting the co-crystalized ligand, and the grid box extended to approximately 20 Å around the ligand.

Virtual screening

Virtual screening was implemented by the glide docking module of the Maestro package using the standard precision (SP docking) mode (Friesner et al., 2004). 6-Mercaptopurine (6-MP) was proved to be a potent inhibitor of SARS PLpro with IC50 = 0.0216 µM (Báez-Santos et al., 2015). The result of the docking run was compared to the values of 6-MP as a reference inhibitor.

Molecular dynamics simulations

Molecular dynamics (MD) simulations were conducted on the most promising drugs to target PLpro using AMBER16 software (Case et al., 2016). In MD simulations, the general AMBER force field 2 (GAFF2) (“GAFF and GAFF2 are public domain force fields and are part of the AmberTools16 distribution, a.f.d.a.h.a.o.i.a.A.M.A.t,”) and AMBER force field (Maier et al., 2015) were utilized to describe drugs and COVID-19 PLpro, respectively. The restrained electrostatic potential (RESP) approach (Bayly et al., 1993) at the HF/6-31G* level was employed to estimate the atomic partial charges of the drugs. The quantum mechanical calculations were conducted using Gaussian09 software (Frisch et al., 2009). The docked drug- PLpro complexes were solvated in a cubic water box with 15 Å distances between the edges of the box and any atom of the complexes. The solvated complexes were then minimized for 5000 steps and gently annealed from 0 K to 300 K over 50 ps. The systems were equilibrated for 1 ns, and production stages were executed over simulation times of 50 ns. Every 10 ps, a snapshot was taken over the production stage, giving 5,000 snapshots for each complex. Periodic boundary conditions and the NPT ensemble were adopted in all explicit MD simulations, including both the equilibration and production stages. Long-range electrostatic interactions under periodic conditions were treated with the Particle Mesh Ewald (PME) method (Darden et al., 1993) with a direct space cut-off of 12 Å. The temperature was held constant at 298 K using Langevin dynamics with the gamma_ln parameter set to 1.0. A Berendsen barostat with a relaxation time of 2ps was employed to control system pressure (Berendsen et al., 1984). All bonds involving hydrogen atoms were constrained using the SHAKE option, and the time step was set to 2fs. The MD simulations were executed with the GPU version of pmemd (pmemd.cuda) in AMBER16 on the CompChem GPU/CPU cluster (hpc.compchem.net).

Binding energy calculations

Coordinates and energy values were collected every 10ps over the production stage for binding energy calculations. Based on the collected coordinates, the binding energies of the drugs with COVID-19 PLpro were evaluated using the molecular mechanics–generalized Born surface area (MM-GBSA) approach (Massova & Kollman, 2000). The MM-GBSA binding free energies were estimated as follows: where the energy term (G) is estimated as: with Evdw, Eele, GGB and GSA as the van der Waals, electrostatic, General Born solvation, and surface area energies, respectively. Entropy contributions were not considered in this study.

Statistical analysis

Several correlation runs were conducted to evaluate the relation between the calculated docking score and the functional parameters of compounds recognition by SARS-CoV-2 PLpro. Pearson’s correlation coefficient and the significance degree were estimated in GraphPad Prism software.

Results

The virtual screening results of a total of 1697 FDA-approved drugs against the SARS-CoV-2 PLpro are provided in Supplementary File 1. The presented parameters include the docking scores, ligand efficiencies, and lipophilic and hydrogen bonding interactions. Table 1 shows the compounds with the highest docking scores, including 26 compounds with a docking score of –7 or higher. The relative docking score and interaction parameters were calculated based on the values of the reference PLpro inhibitor 6-mercaptopurine (6-MP).

Table 1.

The top hits after virtual screening of FDA-approved drugs against SARS-CoV-2 PLpro.

Name	Mol Weight	Docking score	Glide ligand efficiency	Glide lipo	Glide hbond	Clinical application
Fostamatinib disodium	624.423	−8.888	−0.222	−1.508	−0.47	chronic immune thrombocytopenia
Chlorhexidine HCl	578.368	−8.666	−0.255	−1.253	−1.533	antiseptic
Candesartan (Atacand)	440.454	−8.136	−0.247	−1.943	−0.462	angiotensin receptor blocker
NAD+	663.425	−7.983	−0.181	−0.85	−0.696	General tonic
Zanamivir (Relenza)	332.31	−7.812	−0.34	−0.817	−0.86	Influenza neuraminidase inhibitor
Quercitrin	448.377	−7.757	−0.242	−2.31	−0.32	Anticancer, antioxidant
Montelukast Sodium	608.165	−7.572	−0.185	−3.316	−0.599	leukotriene receptor antagonist in asthma
Losartan	462.009	−7.501	−0.25	−2.349	−0.289	angiotensin II receptor blocker
Famotidine (Pepcid)	337.445	−7.411	−0.371	−1.035	−1.034	antacid
Dacomitinib (PF299804, PF-00299804)	469.939	−7.408	−0.224	−1.63	−0.841	anticancer
Ritonavir	720.944	−7.393	−0.148	−1.396	−1.429	Antiviral – HIV protease inhibitor
Flubendazole (Flutelmium)	313.283	−7.359	−0.32	−1.605	−0.342	Anthelmintic
Moroxydine HCl	207.661	−7.334	−0.611	−0.589	−1.184	Broad spectrum antiviral
Tezacaftor VX-661	520.498	−7.333	−0.198	−1.373	−1.195	cystic fibrosis
Phenformin hydrochloride	241.721	−7.232	−0.482	−1.025	−1.148	hypoglycemic agent
Oxiglutatione	612.631	−7.219	−0.18	−0.928	−1.015	Cell protective- antioxidant
Valsartan (Diovan)	435.519	−7.2	−0.225	−2.401	−0.152	angiotensin II receptor blocker
Doripenem Hydrate	438.52	−7.188	−0.266	−1.481	−0.801	carbapenem antibiotic
Guanabenz (WY-8678) Acetate	291.134	−7.186	−0.513	−1.086	−0.948	alpha-2 selective adrenergic agonist - antihypertensive
Dichlorphenamide (Diclofenamide)	305.159	−7.128	−0.446	−0.947	−0.709	carbonic anhydrase inhibitors for glaucoma
LCZ696 Sacubitril	915.979	−7.127	−0.223	−2.514	−0.242	Cardiovascular
Albendazole (Albenza)	265.331	−7.084	−0.394	−1.868	−0.367	Anthelmintic
Fenbendazole (Panacur)	299.348	−7.032	−0.335	−1.709	−0.337	Anthelmintic
Acetazolamide	222.245	−7.024	−0.54	−0.569	−0.456	carbonic anhydrase inhibitors
Sulfaguanidine	214.245	−7.013	−0.501	−0.733	−0.878	Antimicrobial

Among the top hits (Table 1 and Figure 1), there were three general tonics, cell-protective substances, and antioxidants (NAD+, quercetin, and oxyglutathione), three antivirals (ritonavir, moroxydine, and zanamivir), two antimicrobials (doripenem and sulfaguanidine), two anticancer drugs, three anthelmintics from the benzimidazole group, one antacid (famotidine), and other miscellaneous systemically or topically acting drugs.

Figure 1.

Virtual screening and docking of FDA-approved drugs with COVID-19 CoV PLpro. (a) The top 10 hits docked into the PLpro binding site. (b) The molecular surface of SARS-CoV-2 PLpro and its binding site.

As shown in Supplementary File 1, 6-MP occupied the 132nd rank in the virtual screening study. Therefore, the top retrieved hits were more likely to have more potent interactions with SARS-CoV-2 PLpro. The improved docking scores were derived by a corresponding increase in both the hydrophobic scores and hydrogen bonding interactions.

To analyze the interactions of the top hits in comparison with the reference drug, the ligand interaction profiles of drugs and the analyses of the contributing forces were plotted (Figure 1). The binding mode of the compounds showed a U-shaped configuration (Figure 2); it had a U-shaped cavity with the arms of the U hosting TYR261. Zanamivir formed three hydrogen bonds with ASP157, GLU160, and TYR261. NAD + showed a strong network of five hydrogen bonds, one salt bridge, and stacking interactions. Ritonavir showed four hydrogen bonds and a stacking interaction with ARG159 (Figure 3).

Figure 2.

The site of SARS-CoV-2 PLpro binding with the top hits, showing the binding site residues. The binding site is shown in blue mesh, and the binding site residues are colored by element (white for carbon and green for ligands).

Figure 3.

The ligand interactions of zanamivir, NAD+, ritonavir, and 6-mercaptopurine with SARS-CoV-2 PLpro. Charged residue (negative) in pink, positive charged residue in blue, and hydrophobic residues in cyan; hydrogen bonds indicated by purple arrow; stocking interactions indicated by green lines.

Statistical analysis revealed a significant positive correlation between the docking score and the hydrogen bonds, lipophilic interactions, and ligand efficiency scores (p < .001) (Table 2). The molecular weight did not show a significant correlation with the docking score. The positive correlation between both the hydrogen bond and the Glide lipo confirms the role of both hydrogen bonding and lipophilic interactions in the binding to PLpro. This contributed to the improved docking scores of the top hits compared with the score of the previously identified PLpro inhibitor 6-MP.

Table 2.

Pearson’s correlation of the obtained docking score with the drug’s molecular weight and interaction parameters with SARS-CoV-2 PLpro.

	Docking score vs. glide ligand efficiency	Docking score vs. glide lipo	Docking score vs. glide hbond	Docking score vs. MolWeight
Pearson r
r	0.3024	0.2311	0.3602	0.01789
95% confidence interval	0.2585 to 0.3451	0.1855 to 0.2757	0.3181 to 0.401	−0.02978 to 0.06547
R squared	0.09147	0.05341	0.1298	0.00032
P value
P (two-tailed)	<0.0001	<0.0001	<0.0001	0.4620
P value summary	****	****	****	ns
Significant? (alpha = 0.05)	Yes	Yes	Yes	No
Number of XY Pairs	1693	1693	1693	1693

Molecular dynamics simulation

Conformational flexibilities of drug-receptor complexes, solvent effects, and dynamics must be determined to attain dependable drug-receptor–binding affinities (De Vivo et al., 2016; Kerrigan, 2013). Therefore, molecular dynamic (MD) simulations combined with binding energy (MM-GBSA) calculations over 50 ns were conducted for the most promising drugs (i.e. those which had predicted docking scores in association with COVID-19 PLpro). The evaluated average MM-GBSA binding energies over the 50 ns MD simulations are listed in Table 3. According to the calculated MM-GBSA binding energies (ΔG), only 3 of the 10 investigated drugs demonstrated considerable binding energies (ΔG > −37.6 kcal/mol); 7 exhibited relatively weak binding energies in the range of −36.4 ≥ ΔGbinding ≥ −19.5. Interestingly, phenformin, quercetin, and ritonavir showed promising binding affinities toward COVID-19 PLpro, with binding energies (ΔG) of −56.5, −40.9, and −37.6 kcal/mol, respectively (Table 4).

Table 3:

Calculated MM-GBSA binding energies (in kcal/mol) for the top 10 potent drugs against COVID-19 PLpro over 50 ns MD simulations.

No.	Drug name	MM-GBSA Binding energy (kcal/mol)
1	Phenformin	−56.5
2	Quercetin	−40.9
3	Ritonavir	−37.6
4	Montelukast	−36.4
5	Fosmatinib1	−33.5
6	Nadid	−32.9
7	Candestran	−28.9
8	Valsartan	−28.6
9	Zanamivir	−24.7
10	Oxyglutathione	−19.5

Table 4.

Decomposition of MM-GBSA binding energies for phenformin, quercetin and ritonavir in complex with COVID-19 PLpro through 50 ns MD simulations.

Drug Name	Calculated MM-GBSA binding energy (kcal/mol)
	ΔEVDW	ΔEele	ΔEGB	ΔESUR	ΔGgas	ΔGSolv	ΔGbinding
Phenformin	−23.2	−124.1	94.8	−4.0	−147.3	90.7	−56.5
Quercetin	−33.5	−75.0	72.7	−5.2	−108.4	67.6	−40.9
Ritonavir	−45.0	−26.7	40.0	−5.9	−71.7	34.0	−37.6

Post-dynamics analyses

Binding energy decomposition

Decomposition of the average MM-GBSA binding energy over a 50 ns MD simulation was carried out to disclose the nature of dominant interactions for phenformin, quercetin, and ritonavir with COVID-19 PLpro (Table 2). According to energy decomposition results, the electrostatic interactions (Eele) were the dominant forces in the drug–PLpro binding affinity, with values of −124.1, −75.0, and −26.7 kcal/mol for the phenformin–PLpro, quercetin– PLpro, and ritonavir–PLpro complexes, respectively. In addition, the van der Waals (E) contributions were favorable, with values of −23.2, −33.5, and −45.0 kcal/mol for the phenformin–PLpro, quercetin–PLpro, and ritonavir–PLpro complexes, respectively.

Hydrogen bond length

Correlational analyses of the hydrogen bond lengths between phenformin, quercetin, and ritonavir and the key amino acid ASP157 residue were estimated through 50 ns MD simulations and are presented in Figure 4a. What stands out in the figure is the higher stability of phenformin and ritonavir inside the active site of COVID-19 PLpro, with average hydrogen bond lengths of 2.10 and 2.29 Å, respectively. Quercetin formed a stable hydrogen bond through the first 30 ns, with an average length of 1.70 Å. After 30 ns in the MD simulations, there was an increase in the hydrogen bond length, however, so that the average bond length was 2.92 Å throughout the final 20 ns of the MD time.

Figure 4.

(a) Hydrogen bond lengths. (b) CoM distances between the most promising drugs and the key residue amino acid ASP157. (c) Variations in the MM-GBSA binding energies. (d) RMSD of the backbone atoms from the initial structure for phenformin (black), quercetin (red), and ritonavir (blue) in complex with PLpro through 50 ns MD simulations.

Center-of-mass distance

Phenformin, quercetin, and ritonavir inside the active site of COVID-19 PLpro were further scrutinized by measuring the center-of-mass (CoM) distance between the drug and the ASP157 residue through the 50 ns of the MD simulations (see Figure 4b). The CoM distance data in Figure 4b demonstrate that phenformin, quercetin, and ritonavir showed high stability inside the active site of COVID-19 PLpro during the 50 ns of the MD simulations, with average CoM distances of 4.55, 7.57, and 6.05 Å, respectively.

Binding energy per frame

Correlations between the binding energy and time would grant a deeper insight into the stability of the drug–PLpro interactions. The MM-GBSA binding energies were therefore measured per frame for phenformin, quercetin, and ritonavir with PLpro and are displayed in Figure 4c. Data in Figure 4c revealed that the phenformin, quercetin, and ritonavir complexes demonstrated a tight binding affinity over the 50 ns of the MD, with an average ΔG binding of −56.5, −40.9, and −37.6 kcal/mol, respectively. These findings show a promising stability of phenformin, quercetin, and ritonavir in complex with PLpro over the simulated MD time of 50 ns.

Root-mean-square deviation

The root-mean-square deviation (RMSD) was utilized to scrutinize the structural changes in the drug–PLpro complexes over the MD simulated time. The RMSD for the backbone atoms of phenformin, quercetin, and ritonavir in complex with PLpro relative to the starting structures throughout the 50 ns MD simulations were estimated and are summarized in Figure 4d. What stands out in this figure is the backbone of the phenformin–, quercetin–, and ritonavir–PLpro complexes and the demonstrated stability over the 50 ns MD simulations, giving an RMSD of less than 0.35 nm. Overall, these results indicate that the three investigated drugs are tightly bonded in the active site of the PLpro and do not impact the overall topology of the PLpro.

Discussion

This work covers an essential gap in the discovery of new drugs against COVID-19. CoVs encode a limited number of virus-encoded enzymes. Therefore, the therapeutic opportunities are directed to these targets in the hope of finding new anti-SARS drugs. Of these enzymes, PLpro lies at the heart of drug discovery targets owing to its essential function in processing the viral polyprotein, as well as its multiple remote activities that lead to suppression of the host immune response.

In criticizing the top hits, some of these drugs are not suitable for oral administration due to pharmacokinetic restrictions, for example, because of the poor absorption of sulfaguanidine and most of the benzimidazole anthelmintics. Others have important effects on body system functions, such as tezacaftor (cystic fibrosis), guanabenz (adrenoceptor agonist), venetoclax and taxifolin (anticancer drugs), and chlorhexidine (a drug for topical application).

Three previously known antivirals were found among the top hits (ritonavir, moroxydine, and zanamivir). Ritonavir is a viral protease inhibitor approved for treating HIV (Croxtall & Perry, 2010). Ritonavir/lopinavir was successful in treating an index case with severe pneumonia (Lim et al., 2020). Moroxydine is an old antiviral with several new derivatives that are effective against the hepatitis C virus (Magri et al., 2015). It has broad-spectrum antiviral actions involving DNA and RNA viruses (Yu et al., 2016, 2018). Zanamivir is a strong influenza neuraminidase inhibitor (Zhao et al., 2017). It is used in combination with oseltamivir for treating drug-resistant H1N1 influenza (de Mello et al., 2018).

Three previously known antivirals were found among the top hits (ritonavir, moroxydine, and zanamivir). Ritonavir is a viral protease inhibitor approved for treating HIV (Croxtall & Perry, 2010). Ritonavir/lopinavir was successful in treating an index case with severe pneumonia (Lim et al., 2020). Moroxydine is an old antiviral with several new derivatives that are effective against the hepatitis C virus (Magri et al., 2015). It has broad-spectrum antiviral actions involving DNA and RNA viruses (Yu et al., 2016, 2018). Zanamivir is a strong influenza neuraminidase inhibitor (Zhao et al., 2017). It is used in combination with oseltamivir for treating drug-resistant H1N1 influenza (de Mello et al., 2018).

Theoretically, neuraminidase inhibitors are not feasible for application against SARS-CoV-2 due to the absence of a molecular target. However, the general goal in identifying a lead compound is finding a chemical structure that can theoretically match a receptor. Several neuraminidase inhibitors theoretically matched the receptor of SARS-CoV-2 PLpro and could be inhibitors of the viral protease. The presence of ACE receptor blockers among the top hits suggests that they could be used in SARS-CoV-compromised hypertensive patients. In conclusion, after a virtual screening campaign of 1697 FDA-approved drugs, 10 drugs were selected for MD and energy composition analysis. Of these, three drugs showed favorable profiles suggesting that they could be repurposed for COVID-19 treatment. The top three hits, phenformin, quercetin, and ritonavir, can be a basis for future anti-SARS CoV-2 studies.