Computational drug repurposing for the identification of SARS-CoV-2 main protease inhibitors.

Accepted 7 July 2020ABSTRACT Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is the virus responsible for the known COVID-19 disease. Since currently no definitive therapies or vaccines for the SARS-CoV-2 virus are available, there is an urgent need to identify effective drugs against SARS-CoV-2 infection. One of the best-known targets available is the main protease of this virus, crucial for the processing of polyproteins codified by viral RNA. In this work, we used a computational virtual screening procedure for the repurposing of commercial drugs available in the DrugBank database as inhibitors of the SARS-CoV-2 main protease. Molecular docking calculations and molecular dynamics (MD) simulations have been applied. The computational model was validated through a self-docking procedure. The screening procedure highlighted five interesting drugs that showed a comparable or higher docking score compared to the crystallographic compound and maintained the protein binding during the MD runs. Amongst these drugs, Ritonavir has been used in clinical trials with patients affected by COVID-19 and Nelfinavir showed anti-SARS-CoV-2 activity. The five identified drugs could be evaluated experimentally as inhibitors of the SARS-CoV-2 main protease in view of a possible COVID-19 treatment. Communicated by Ramaswamy H. Sarma.

Introduction

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) caused the well-known coronavirus disease of 2019, COVID-19. In a few months since its outbreak in the end of 2019, it caused more than 5 million cases and 340,000 deaths around the world, involving 216 countries, areas or territories (https://www.who.int/emergencies/diseases/novel-coronavirus-2019). First identified in Wuhan (Hubei, China) in December 2019, this infection was declared a pandemic by the World Health Organization in March 2020. USA is the most affected country at the time of writing, with 1.6 million cases and 100,000 deaths (https://www.cdc.gov/coronavirus/2019-ncov/cases-updates/cases-in-us.html). Italy was seriously affected, with more than 239,000 cases and almost 35,000 deaths (http://www.salute.gov.it/portale/nuovocoronavirus/homeNuovoCoronavirus.jsp?lingua=english). The restrictive measures applied by the Italian Government have led to a decrease in the number of new positive cases and deaths.

SARS-CoV-2 belongs to Coronaviridae family, discovered in the 1960s. SARS-CoV-2 is an enveloped positive-sense single-stranded ribonucleic acid (RNA) virus (Kahn & McIntosh, 2005), bearing club-shaped spike peplomers that cover its surface and confer to it the typical crown appearance (Goldsmith et al., 2004).

COVID-19 causes mild to moderate respiratory illness in most of the infected people. These patients are able to recover without requiring special treatment, whilst elderly, or people affected by health problems like cardiovascular diseases, diabetes, chronic respiratory diseases and cancer are more susceptible to develop serious symptoms. The most common symptoms are fever, dry cough and tiredness, whilst less common symptoms include loss of taste or smell, aches, sore throat, diarrhoea, conjunctivitis and headache. The worst-case scenario involves difficulty in breathing or shortness of breath, chest pain or pressure, loss of speech or movement (https://www.who.int/health-topics/coronavirus#tab=tab_1).

No definitive therapies or vaccines for the SARS-CoV-2 virus infection are currently available. However, many clinical trials are ongoing to evaluate potential treatments and several viral targets are under investigation with the aim to identify novel pharmacological approaches. Amongst these, one of the best-characterized drug targets is the main protease (Mpro or 3CLpro) (Anand et al., 2003), an enzyme essential for the processing of the polyproteins that are codified by the viral RNA (Zhang et al., 2020). An additional advantage deriving from the inhibition of this enzyme is that no human protease shows a similar cleavage specificity, therefore Mpro inhibitors are expected to be selectively toxic for the virus and not for the host cell. In this context, evaluation of commercially available drugs that have already passed clinical trials would be a fast way to identify active molecules with no need to invest too much time and money in R&D activities.

Based on these considerations, a structure based virtual screening approach for the repurposing of commercially available drugs was applied, hoping to speed up the discovery of compounds for COVID-19 treatment. In particular, a computational study was performed aimed at identifying Mpro inhibitors amongst FDA approved drugs reported in the DrugBank database (Wishart et al., 2018), using docking calculations and molecular dynamics (MD) simulations.

Materials and methods

Ligand preparation

The DrugBank database was downloaded and prepared using LigPrep with default settings (Schrödinger, 2018; Wishart et al., 2018), using OPLS3e as force field, a ionization pH value of 7.0 ± 2.0 performed through Epik, and desalting the ligand. Tautomers were generated for each ligand retaining the specified chiralities, and the only one solution per ligand was selected.

Protein preparation

The selected crystal structure 6W63 (Mesecarr, 2020) was downloaded from the Protein Data Bank (PDB) (Berman et al., 2000) and prepared using the Protein Preparation wizard tool of the Schrödinger suite with default settings (Sastry et al., 2013). Briefly, the bond orders were assigned, giving a zero-order bond to metal bonds, whilst disulphide bonds and hydrogens were added. A pH value of 7.0 was used both during the ionization step (performed through Epik) and the pKa values predictions (performed through PROPKA). Water molecules more than 5.0 Å away from het groups or with less than three H-bonds to non-water molecules were removed. Finally, a 25 Å grid was generated with GlideGrid, with the ligand positioned in the centre of the box with coordinates x= −20.59, y = 18.08 and z= −27.0 relative to the crystal, automatically obtained by the software selecting the centroid of the ligand.

Docking procedure

To validate the protocol, a self-docking of the crystallographic compound (N-(4-tert-butylphenyl)-N-[(1R)-2-(cyclohexylamino)-2-oxo-1-(pyridin-3-yl)ethyl]-1H-imidazole-4-carboxamide, referred to as X77 in the PDB) was run. To this aim, the ligand X77 was extracted from the crystal, converted to 2-D SMILE format and then to 3-D SDF format. In this way, the ligand conformation and all its atomic coordinates were computed starting from scratch using LigPrep. The ligand, then, was docked back into the crystal. The scoring function chosen for the validation was Glide XP (Version 8.1) (Friesner et al., 2006). The docking calculation on the DrugBank database was performed using Glide SP (Version 8.1), followed by Glide XP docking calculation. Apart from the precision (which is set to standard for Glide SP and extra for Glide XP), both the calculations were performed leaving untouched the default settings, i.e. using a flexible ligand sampling, and sampling nitrogen atom inversion, ring conformations and torsions for amide groups only, penalizing non-planar conformation. Epik state penalties were added to the docking score as well, and the post-docking minimization was performed on 10 poses per ligand with 0.5 kcal/mol.

Molecular dynamics simulations

The compounds obtained from the Glide XP calculation underwent molecular dynamics simulations, performed with Amber 18 software. The complexes were prepared using Amber Tools 18 software. The minimization calculations were performed using Sander software, and the next steps were realized using pmemd software (Case et al., 2018).

Three minimization phases were performed, using General Amber Force Field and reducing the degree of constraints applied to each phase. In the first one, the constraints were limited to the protein and the ligand and the TIP3P water solvent was minimized with steepest descent gradient method for 1000 steps, followed by 1000 steps using conjugate gradient. The second minimization removed the constraint from the residues side chains and the ligand, leaving them only on the protein backbone. The same method was used for the minimization. The last phase was done without constraints for a total amount of 5000 steps, of which 2000 with steepest descent gradient and 3000 with conjugate gradient. The system was then gradually heated from 0 to 100 K using Langevin thermostat for 25,000 steps, setting constant volume. This stage was followed by a second heating up to 300 K using Langevin thermostat at constant pressure, monitored with Berendsen barostat. The two heating steps were carried out using a constraint on the protein backbone. Once the system was at 300 K, two equilibration phases were performed. The first one was constituted by 125,000 steps at temperature and pressure of 300 K and 1 atm, respectively. The second one was constituted by 250,000 steps at 300 K temperature and 1 atm pressure without constraint. Finally, the molecular dynamics production was performed at 300 K temperature and 1 atm pressure without constraint for 10,000,000 steps for a total amount of 20 ns simulation time. Three replicas were run per each molecule. The MMPBSA.py.MPI script was applied to each replica for the generation of mmpbsa and mmgbsa values (Miller et al., 2012).

Sequences alignment

The primary sequences alignment of the SARS-CoV-2 and HIV protease was performed using the Muscle web server (Edgar, 2004).

Results and discussions

Self-docking procedure

Recently, the first released COVID-19 main protease crystal structure was used for drug repurposing using a docking calculation study (Kandeel & Al-Nazawi, 2020). During the pandemic, several COVID-19 main protease crystal structures in complexes with molecular fragments or covalent inhibitors were released and reported in the PDB. Amongst them, 6W63 was co-crystallized with the broad-spectrum, non-covalent inhibitor X77 which occupied the catalytic site Cys145, coinciding with the binding site of several other inhibitors. For these reasons, 6W63 was chosen for this work. The protocol was validated using a self-docking procedure. The crystallographic compound X77 was extracted from the complex and docked into the protein, using a minimized 3-D structure generated from a monodimensional Smile file of the ligand. The score obtained from the self-docking was −5.8 kcal/mol and the pose was very similar to the crystallographic one. Both compounds formed H-bond interactions with Gly143, Ser144 and His163 and are in front of the catalytic residue Cys145 (Figure 1). Thus, although X77 was able to generate an interaction with Glu166, that the self-docking pose was not able to predict due to the rotation of the amide group, these results allowed us to be confident with the reliability of the docking procedure.

Figure 1.

Crystallographic (cyan) and self-docking (magenta) poses comparison. Catalytic residue Cys145 is highlighted in orange. Images realized using the PyMol software, version 1.8 (Schrodinger, 2015).

Docking of DrugBank compounds

The prepared DrugBank database was filtered through a fast docking of the ligands into the binding pocket, performed using Glide SP, to reduce the number of compounds and, consequently, the calculation time of later stages. This filtration reduced the number of compounds from an initial value of around 10,000 to 673. Then, these compounds underwent the Glide XP docking calculation, which resulted in 36 docked and scored compounds, whilst the remaining were rejected. 14 drugs showed an improved XP score, lower than the crystallographic reference. The XP docking score ranged from −11.6 to −2.7 kcal/mol, with the crystallographic reference compound X77 showing a docking score of −5.8 kcal/mol. The compounds showed a binding mode comparable to that of inhibitor X77, despite the structural differences. As an example, the binding mode of Ritonavir is highlighted in Figure 2. Ritonavir showed a docking score of −6.8 kcal/mol, a slightly better score compared to X77. In particular, the pyridine ring and the tert-butyl group of X77 were overlapped, respectively, with a phenyl ring and the isopropyl group in Ritonavir, whilst the amide oxygen that gave an H-bond interaction with Gly143 is replaced by a hydroxy group in Ritonavir.

Figure 2.

Superposition of docking result of Ritonavir (yellow) and crystallographic compound X77 (cyan).

Molecular dynamics and trajectories analysis

Due to these findings, all the 36 XP scored compounds and the inhibitor X77 underwent MD simulation in order to evaluate their behaviour in the binding pocket. Three replicas of 20 ns duration time were run per each molecule. All the trajectories were analysed through mmpbsa calculation to rank molecules. The inhibitor X77 showed a mmgbsa value of −47.8 ± 0.8 kcal/mol and mmgbsa value of −35.9 ± 0.6 kcal/mol. The compounds, instead, presented mmgbsa values ranging from −69.2 ± 9.7 to −13.3 ± 16.3 kcal/mol and mmpbsa values ranging from −55.6 ± 12.8 to 0.6 ± 25.0 kcal/mol. The compounds were ranked using both mmpbsa and mmgbsa values using a consensus scoring rank by rank procedure that allowed us to obtain the final ranking list (Table 1).

Table 1.

Consensus scoring rank by rank, mmpbsa and mmgbsa results (reported as mean values of three replicas).

Compound	Clinical use	Docking score	mmgbsa	DevST	mmpbsa	DevST	Rank
Leuprolide	Anticancer agent	−9.8	−69.2	9.7	−55.6	12.8	1
Atracurium besylate	Muscle relaxant	−4.1	−64.9	6.3	−44.5	4.4	2
Ritonavir	Antiviral agent	−6.8	−54.5	4.5	−36.5	2.9	3
Deferoxamine	Iron poisoning	−8.1	−44.4	2.2	−36.8	1.1	4
Ligand	–	−5.8	−47.8	0.8	−35.9	0.6	5
Valrubicin	Anticancer agent	−7.2	−49.6	3.1	−29.7	2.8	6
Teniposide	Anticancer agent	−9.4	−49.9	2.9	−28.6	1.4	7
Aprepitant	Prevent nausea and vomiting	−7.4	−44.4	5.7	−30.2	4.8	8
Nelfinavir	Antiviral agent	−6.5	−39.5	5.6	−32.3	4.7	9
Cefpiramide	Antibiotic	−7.3	−38.3	7.6	−35.5	3.2	10
Ergotamine	Migraine	−7.5	−41.1	6.9	−28.8	4.0	11
Pentagastrin	evaluation of gastric function	−9.8	−47.0	7.7	−26.8	9.5	12
Calcitriol	Vitamin	−6.6	−40.2	7.5	−28.0	3.6	13
Trazodone	Treatment of depression	−2.7	−41.8	2.8	−25.9	3.7	14
Amikacin	Antibiotic	−8.6	−34.9	14.4	−31.7	17.5	15
Montelukast	Anti-asthma	−6.5	−37.2	3.8	−27.0	4.7	16
Raloxifene	Osteoporosis	−3.6	−40.8	3.4	−24.4	3.7	17
Verapamil	Hypertension	−3.1	−39.0	8.6	−25.2	4.8	18
Remikiren	Hypertension	−8.6	−36.9	3.5	−26.7	1.8	19
Trimethobenzamide	Prevent nausea and vomiting	−5.7	−34.6	5.5	−25.5	3.0	20
Pravastatin	Hyperlipidemia	−9.0	−32.7	8.7	−21.5	7.5	21
Raltitrexed	Anticancer agent	−6.7	−34.4	7.0	−21.2	5.3	22
Fluphenazine	Antipsychotic	−3.3	−36.3	0.8	−18.1	2.3	23
Lisinopril	Hypertension	−6.4	−29.5	0.5	−24.3	1.5	24
Riboflavin	Vitamin	−9.2	−35.2	0.4	−19.8	3.0	25
Epirubicin	Anticancer agent	−4.2	−32.4	9.6	−21.0	4.5	26
Acarbose	Type 2 diabetes	−11.6	−37.6	2.2	−15.5	2.5	27
Tamsulosin	Benign prostatic hypertrophy	−3.5	−30.4	5.0	−21.0	3.2	28
Capreomycin	Antibiotic	−6.1	−13.3	16.3	−23.4	22.3	29
Celecoxib	NSAID	−6.7	−25.4	2.5	−20.0	2.7	30
Chlorthalidone	Hypertension	−5.0	−27.7	0.4	−17.4	0.2	31
Ertapenem	Antibiotic	−4.8	−26.3	7.9	−17.0	2.2	32
Mefloquine	Antimalarials	−4.9	−22.7	1.8	−17.5	0.9	33
Nabumetone	NSAID	−5.9	−21.7	3.8	−17.7	2.7	34
Cerivastatin	reduction of LDL cholesterol	−7.7	−26.0	4.7	−15.7	1.5	35
Dyphylline	Anti-asthma	−7.0	−15.3	5.7	−11.6	4.3	36
NADH	Parkinson and Alzheimer	−9.0	−17.2	15.1	0.6	25.0	37

Note. Crystallographic compound is highlighted in grey.

As expected, compounds showing lower docking score compared to the reference inhibitor gave rise to unstable MD trajectories associated to low mmgbsa and mmpbsa values. The RMSD analysis showed that top ranked drugs presented convergent RMSD in the three replicas and a visual analysis highlighted that molecules were stable into the binding pocket. Looking at the final rank by rank score, four compounds achieved a higher score if compared to X77. It is worth noting that amongst the top 10 ranked structures, two different HIV protease inhibitors (Ritonavir and Nelfinavir) and three antitumor agents (Leuprolide, Valrubicin and Teniposide) are represented (Figure 3).

Figure 3.

Antiviral and anticancer agents 2-D structures.

These drugs displayed both docking and MD analysis scores comparable to or higher than the crystallographic compound X77. An extended MD analysis was performed on these compounds, to better understand their behaviour inside the binding pocket of the enzyme. Three replicas of 100 ns were run for each compound. The visual analysis of the extended simulations (Figure 4) highlighted that in these three replicas of MD Leuprolide and Nelfinavir maintained the same binding mode throughout all the simulations, with all the replicas overlapping for the whole timeframe. In particular, Leuprolide showed a strong binding with the active site formed through the portion between the central tyrosine and the terminal proline, whilst the remaining part of the molecule floats in the solvent. This could indicate that the latter is not crucial for the binding in this complex. Ritonavir and Teniposide maintained the binding pose during two of the three replicas. Ritonavir, during Replica 1, presents a high value of RMSD in frames ranging from 550 to 850, related to a 360° rotation of the terminal thiazole around the carbamate group. Valrubicin, instead, during the visual analysis of its trajectories, showed the movement of the tetracyclic portion towards the solvent, in spite of RMSD values comparable to the other compounds. This shift causes a drifting of the 2-tetrahydropyranyl ether out from the Cys145 pocket. Overall, there is correlation between RMSD, mmgbsa and mmpbsa scores (reported in Table 2) and the visual analysis of the compounds. The only exceptions are Nelfinavir, which showed slightly lower mmgbsa and mmpbsa values despite its remarkable RMSD and visual analysis results, and Valrubicin, which presented a mirrored trend, with better mmgbsa and mmpbsa scores and suboptimal RMSD and visual analysis results.

Figure 4.

RMSD analysis.

Table 2.

100 ns molecular dynamic simulations of the selected compounds, mmpbsa and mmgbsa results (reported as mean values of three replicas).

Compound	mmgbsa	DevST	mmpbsa	DevST
Leuprolide	−59.6	7.5	−46.1	8.1
Nelfinavir	−40.1	4.7	−32.2	5.3
Ritonavir	−43.1	8.6	−26.9	6.7
Teniposide	−46.3	4.7	−27.9	4.7
Valrubicin	−47.2	5.3	−31.6	5.7

Primary sequences alignment

Ritonavir has been used in clinical trials with patients affected by COVID-19 (Cao et al., 2020) and Nelfinavir showed anti-SARS-CoV-2 activity (Musarrat et al., 2020; Yamamoto et al., 2020; ). This independent result is a further confirmation of the reliability of our computational procedure. Based on these findings, we evaluated the similarity between SARS-CoV-2 main protease and HIV protease using the Muscle server. We noticed that the percent of identity amongst the two sequences was 27.55%, a low value indicating that the two proteins are structurally different, and hence can be considered to be different targets. Consequently, the identification through our virtual screening of two HIV protease inhibitors, Ritonavir and Nelfinavir, as potential Mpro inhibitors that could be used to combat SARS-CoV-2 infection is far from being an obvious result and adds credibility to this approach aimed at paving the way for new therapeutic treatments of COVID-19.

Conclusions

In summary, in this work we applied a structure based virtual screening of the DrugBank database for the repurposing of commercial drugs as SARS-CoV-2 main protease inhibitors. We identified two antiviral agents, Ritonavir and Nelfinavir, whose mode of action entails the inhibition of HIV-1 protease, and three anticancer agents, Leuprolide, Valrubicin and Teniposide. After the experimental evaluation of their anti-SARS-CoV-2 activity, these compounds could represent a good starting point for structural optimization and synthesis of new analogues within a medicinal chemistry programme aimed at defeating COVID-19.