In Silico Evaluation of Prospective Anti-COVID-19 Drug Candidates as Potential SARS-CoV-2 Main Protease Inhibitors.

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a recently emanating human infectious coronavirus that causes COVID-19 disease. On 11th March 2020, it has been announced as a pandemic by the World Health Organization (WHO). Recently, several repositioned drugs have been subjected to clinical investigations as anti-COVID-19 drugs. Here, in silico drug discovery tools were utilized to evaluate the binding affinities and features of eighteen anti-COVID-19 drug candidates against SARS-CoV-2 main protease (Mpro). Molecular docking calculations using Autodock Vina showed considerable binding affinities of the investigated drugs with docking scores ranging from - 5.3 to - 8.3 kcal/mol, with higher binding affinities for HIV drugs compared to the other antiviral drugs. Molecular dynamics (MD) simulations were performed for the predicted drug-Mpro complexes for 50 ns, followed by binding energy calculations utilizing molecular mechanics-generalized Born surface area (MM-GBSA) approach. MM-GBSA calculations demonstrated promising binding affinities of TMC-310911 and ritonavir towards SARS-CoV-2 Mpro, with binding energy values of - 52.8 and - 49.4 kcal/mol, respectively. Surpass potentialities of TMC-310911 and ritonavir are returned to their capabilities of forming multiple hydrogen bonds with the proximal amino acids inside Mpro's binding site. Structural and energetic analyses involving root-mean-square deviation, binding energy per-frame, center-of-mass distance, and hydrogen bond length demonstrated the stability of TMC-310911 and ritonavir inside the Mpro's active site over the 50 ns MD simulation. This study sheds light on HIV protease drugs as prospective SARS-CoV-2 Mpro inhibitors.

Introduction

The coronavirus 2019 (COVID-19) pandemic that has emerged from severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) is a universal crisis [1, 2]. According to the World Health Organization (WHO), there were > 1,500,000 accumulative cases universally, with a ~ 5.9% case mortality reported in April 2020 [3]. One of the most significant challenges is the absence of a specific vaccine and treatment for COVID-19 [4]. Repositioning of approved drugs is therefore required as the fast track in the fight against COVID-19 infection [5-7]. Since January 2020, several clinical trials have tested antiviral, anti-inflammatory, and anti-malarial drugs for treating COVID-19 [8-10]. Among these drugs, remdesivir (GS-5734) is the most promising drug candidate to combat COVID-19 [11, 12]. Consequently, remdesivir has been recently declared to have a positive influence in a clinical trial for the treatment of COVID-19 [13]. Remdesivir, a nucleotide analog antiviral drug, targets the RNA-dependent RNA polymerase (RdRp) enzyme of SARS-CoV-2 and, in turn, prevents viral replication [11]. Prevention of the SARS-CoV-2 replication could also be achieved by targeting the viral main protease (Mpro) (also called 3-chymotrypsin-like protease (3CLpro)) and papain-like protease (PLpro). Because of the essential role of Mpro in the viral life cycle, numerous experimental and in silico studies have attempted to identify small molecules, natural products, and repurposed drugs as potential SARS-CoV-2 Mpro inhibitors [14-23]. Until now, the outcomes of the initiated clinical trials for protease inhibitors have not yet been released [18, 24, 25]. Therefore, the present study was set to investigate eighteen drugs in clinical development as potential Mpro inhibitors using in silico drug discovery techniques. The investigated drugs included anti-influenza drugs, anti-malarial drugs, anti-HIV drugs, anti-inflammatory drugs, anti-HCV drugs, anti-Ebola drugs, and anti-genitourinary infections drugs. Molecular docking calculations were first performed to predict the binding poses of the investigated drugs against SARS-CoV-2 Mpro. The binding affinities and features were highlighted. Molecular dynamics (MD) simulations for 50 ns were executed on the predicted docked drug-Mpro complexes, and a molecular mechanical-generalized Born surface area (MM-GBSA) approach was utilized to evaluate the binding energies (∆Gbinding). Post-dynamics analyses, including structural and energetic aspects, were executed to obtain insight into the stability and affinity of the drugs as potential Mpro inhibitors. The presented results are promising and proposing prospect inhibition for obtainable therapeutics against COVID-19.

Computational Methodology

Drug Preparation

Eighteen repurposed drugs for COVID-19 treatment were retrieved in SDF format from the DrugBank database before molecular docking calculations [26, 27]. Omega software was utilized to generate the 3D structures of the drugs [28, 29], and their geometrical structures were then minimized with the MMFF94S force field utilizing SZYBKI software [30]. The 2D chemical structures of the investigated anti-COVID-19 drug candidates are depicted in Table 1.

Table 1

DrugBank code, 2D chemical structure, chemical description, original usage, and mechanism of action of the repurposed drugs in clinical investigation to combat COVID-19

Main Protease Preparation

For molecular docking and molecular dynamics calculations, the crystal structure of SARS-CoV-2 main protease (Mpro; PDB code: 6LU7 [18]) was taken as a template. Ions and water molecules were removed. H++ server was employed to investigate the protonation state of Mpro. As well, all missing hydrogen atoms were added [31]. For molecular docking calculations, the Autodock protocol was utilized to prepare the pdbqt file of Mpro [32]. The pKa for SARS-CoV-2 Mpro amino acid residues was estimated under physical conditions of salinity = 0.15, pH 6.5, external dielectric = 80, and internal dielectric = 10.

Molecular Docking

Autodock Vina software was used to perform all molecular docking calculations [33]. All docking parameters were conserved to the default, except the exhaustiveness parameter was set to 200. The binding site was realized by a docking box around the active site with XYZ dimensions of 25 Å × 25 Å × 25 Å and a spacing value of 1.00 Å. The center of the grid was positioned at − 13.069, 9.74, 68.49 (XYZ coordinates) for SARS-CoV-2 Mpro.

Molecular Dynamics

AMBER16 software was utilized to perform molecular dynamics (MD) simulations on the docked structures of the investigated drugs inside the active site of Mpro [34]. The minutiae description of the applied MD simulations was demonstrated in Ref. [7, 14]. Concisely, Mpro was described using AMBER force field FF14SB [35], while the general AMBER force field (GAFF2) was utilized to describe the drugs [36]. The restrained electrostatic potential (RESP) approach with the assistance of Gaussian09 software was applied to assign the atomic partial charges of the investigated drugs [37, 38]. The drug-Mpro complexes were solvated in a cubic water box with 15 Å distances between the edges of the box and any atom of drug or drug-Mpro complexes. The solvated drug-Mpro complexes were minimized for 5000 steps, gradually heated from 0 to 300 K over 50 ps, and equilibrated for 1 ns. The NPT ensembles were adopted, and drug-Mpro complexes were simulated for 50 ns. Pmemd.cuda implemented in AMBER16 was utilized to perform all molecular dynamics simulations. CompChem GPU/CPU cluster (hpc.compchem.net) was applied to carry out molecular docking calculations and molecular dynamics simulations. 2D and 3D of the drug-Mpro interactions were visualized using the Discovery studio module of Biovia software (Dassault Systemes of France).

MM-GBSA Binding Energy

The binding energies of the investigated anti-COVID-19 drug candidates with SARS-CoV-2 Mpro were estimated using molecular mechanics-generalized Born surface area (MM-GBSA) approach with modified GB model (igb = 2) implemented in AMBER16 software [39]. For the MM-GBSA calculations, binding energy was calculated according to uncorrelated snapshots collected every 10 ps over the production run. The binding energy (ΔG) was evaluated as follows: where the energy term (G) is estimated as:E and E are van der Waals and electrostatic energies, respectively. G is the electrostatic solvation free energy calculated from the generalized Born equation and G is the nonpolar contribution to the solvation free energy from the solvent-accessible surface area (SASA). A single-trajectory approach was used, in which the coordinates of each drug-receptor, receptor and drug were extracted from a single trajectory. For all investigated drugs, entropy contributions were neglected.

Results and Discussion

Several repurposed drugs alone or in combinations have been subjected to clinical trials to treat COVID-19 [8, 40]. The mechanism of action of most of these drugs is to target the viral replication process or block viral entry into the host cell. Among these drugs, sixteen approved and investigational drugs were originally developed to treat influenza, HIV, HCV, and other respiratory infections. Besides, two approved anti-malarial drugs were proposed for clinical investigation as prospective anti-COVID-19 drugs. The chemical structure, original usage, and mechanism of action of these eighteen anti-COVID-19 drug candidates are summarized in Table 1.

The efficacy and safety of these drugs as promising anti-COVID-19 drugs are still under evaluation. In the present study, binding affinities, features, and stabilities of these anti-COVID-19 drug candidates were evaluated against SARS-CoV-2 main protease (Mpro) using in silico drug discovery techniques.

Molecular Docking Calculations

Molecular docking is used as a fundamental tool in the drug discovery pipeline [41]. In the present study, Autodock Vina software was applied to perform all molecular docking calculations and predict the binding modes of the repurposed drugs with SARS-CoV-2 Mpro. The predicted binding affinities and features of the investigated drugs towards Mpro are listed in Table 2. The 2D and 3D representations of interactions of the inspected drugs with the key amino acid residues of SARS-CoV-2 Mpro are illustrated in Fig. S1.

Table 2

Calculated docking scores (in kcal/mol) and binding features for the COVID-19 drug candidates against SARS-CoV-2 main protease (Mpro)

Drug	Docking score (kcal/mol)	Binding Features(Hydrogen bond length in Å)
Umifenovir	− 6.2	GLU166 (2.42 Å)
Favipiravir	− 5.3	HIS163 (2.22 Å), CYS145 (2.64 Å), SER144 (2.26 Å), GLY143 (2.08 Å), LEU141 (1.89 Å)
Triazavirin	− 5.8	HIS163 (2.12 Å)
Oseltamivir	− 5.8	GLN189 (2.56 Å), HIS41 (2.75 Å), CYS145 (2.06 Å), HIS163 (2.46 Å)
Baloxavir marboxil	− 6.1	HIS41 (2.92 Å), GLY143 (2.52 Å)
Hydroxychloroquine	− 6.2	HIS164 (2.62 Å)
Chloroquine	− 5.8	HIS164 (2.57 Å)
TMC-310911	− 8.3	CYS145 (2.34 Å), SER 144 (3.02 Å), GLU166 (1.89 Å)
Ritonavir	− 7.7	GLN189 (2.63 Å), HIS163 (2.01 Å), GLU166 (2.38 Å), GLY143 (1.83 Å)
Lopinavir	− 8.0	GLY143 (2.20 Å), GLU166 (2.93 Å)
Darunavir	− 7.5	THR24 (2.44 Å), THR25 (2.67 Å), GLY143 (2.34 Å), GLN189 (2.19 Å), GLU166 (2.13 Å)
Tenofovir alafenamide	− 7.5	SER144 (2.49 Å), CYS145 (2.50 Å), GLU166 (2.47 Å), THR190 (2.18 Å)
Emtricitabine	− 5.8	GLU166 (2.65 Å), PHE140 (2.22 Å), HIS163 (2.22 Å), SER144 (2.74 Å), CYS145 (2.65, 2.99 Å)
Methylprednisolone	− 6.7	CYS145 (2.48 Å), LEU141 (2.21 Å)
Ribavirin	− 6.3	LEU141 (1.83 Å), SER144 (2.76 Å), CYS145 (2.57 Å), HIS164 (2.42 Å)
Remdesivir	− 8.2	GLY143 (2.38 Å), HIS163 (2.28 Å), GLU166 (2.46 Å)
Galidesivir	− 7.1	LEU141 (1.93 Å), SER144 (2.25 Å), HIS163 (2.17 Å), GLU166 (2.77 Å)
Azithromycin	− 6.4	GLU166 (2.25 Å)

As can be seen from data in Table 2, the docking scores of the investigated drugs with SARS-CoV-2 Mpro ranged from − 5.8 to − 8.3 kcal/mol. Inspecting the drug-Mpro interactions revealed that most of the drugs share a similar binding pose inside the active site of Mpro, exhibiting a fundamental hydrogen bond with GLU166. Further interactions, including hydrogen bonds, hydrophobic interactions, and pi-based interactions, were also noticed between the drugs and the proximal amino acids inside the active site of SARS-CoV-2 Mpro (Fig. S1). It is also worth noting that a small group of the examined drugs could not form such a fundamental hydrogen bond with GLU166. This is returned to the preferability of these drugs to form a parallel hydrogen bond with HIS164 or HIS163, as in the case of chloroquine and triazavirin (Table 2). Several studies have indicated the vital role of the catalytic CYS145 residue in suppressing the protease activity [42, 43]. As shown in Table 2, seven out of the eighteen examined drugs exhibited a potential hydrogen bond with the catalytic CYS145 residue with bond lengths in the range of 2.06 − 2.99 Å.

Among the examined drugs, HIV protease inhibitors demonstrated higher binding affinities against SARS-CoV-2 Mpro, compared to the other antiviral and anti-malarial drugs. For instance, TMC-310911 showed the highest binding affinity against Mpro with a docking score of − 8.3 kcal/mol. Structural insights into the binding mode of the TMC-310911 with the Mpro demonstrated that the sulfone moiety forms a hydrogen bond with the backbone NH of GLU166 residue with a bond length of 1.89 Å (Fig. 1 and Table 2). While, the carbamate moiety of TMC-310911 interacts with the backbone atoms of SER144, and CYS145 via hydrogen bonds of lengths 3.02 and 2.34 Å, respectively (Fig. 1 and Table 2). The strong interaction of TMC-310911 with the conserved residue GLU166 of Mpro has been recently reported with a docking score of − 7.1 kcal/mol [25].

Fig. 1

2D representations of the predicted binding poses of the scrutinized HIV protease drugs inside the active site of SARS-CoV-2 main protease (Mpro)

The other four examined HIV protease drugs ─namely lopinavir, ritonavir, darunavir, and tenofovir alafenamide─ demonstrated promising binding affinities against Mpro with docking scores of − 8.0, − 7.7, − 7.5, and − 7.5 kcal/mol, respectively, forming two, four, five, and four hydrogen bonds with the key amino acids of Mpro, respectively (Fig. 1). Inspecting the binding mode of ritonavir with the Mpro revealed that (i) the C=O of the carbamoyl moiety demonstrates a hydrogen bond with the backbone NH of GLU166, with a bond length of 2.38 Å, (ii) the thiazole ring interacts with the imidazole ring of HIS163 via a hydrogen bond of length 2.01 Å, (iii) the C=O of the urea group forms a hydrogen bond with the backbone NH of GLY143 with a bond length of 1.83 Å, and (iv) the OH group participates in a hydrogen bond with GLN189 with a bond length of 2.63 Å (Fig. 1 and Table 2).

In agreement with the current results, the higher potency of ritonavir, compared to lopinavir, has been recently reported with SARS-CoV-2 Mpro [44].

Interestingly, remdesivir, an RNA-dependent RNA polymerase (RdRp) inhibitor developed for Ebola-virus, demonstrated the second-highest binding affinity towards SARS-CoV-2 Mpro with a docking score of − 8.2 kcal/mol. The other examined anti-Ebola drug ─namely galidesivir─ showed a satisfactory binding affinity with a docking score of − 7.1 kcal/mol with Mpro.

For anti-influenza drugs, relatively weak binding affinities were observed with SARS-CoV-2 Mpro, with docking scores of − 6.2, − 6.1, − 5.8, − 5.8, and − 5.3 kcal/mol for umifenovir, baloxavir marboxil, triazavirin, oseltamivir, and favipiravir, respectively. Besides, the anti-malarial, HCV, and anti-informatory drugs showed relatively weak potency towards Mpro (docking scores of − 6.2, − 5.8, − 6.7, − 6.3, and − 6.4 kcal/mol for hydroxychloroquine, chloroquine, methylprednisolone, ribavirin, and azithromycin, respectively).

In summary, the results of molecular docking calculations suggested HIV protease drugs and remdesivir as promising SARS-CoV-2 Mpro inhibitors.

Molecular Dynamics Simulations

It has been demonstrated that solvent effects, conformational flexibilities of drug-receptor complexes, and dynamics must be applied to improve the reliability of the predicted ligand–protein binding energies. Therefore, molecular dynamics (MD) simulations combined with MM-GBSA binding energy calculations were carried out for the repurposed drugs complexed with SARS-CoV-2 Mpro.

10 ns MD simulations

Eighteen repurposed drugs in clinical development to treat COVID-19 were further inspected using MD techniques over a simulation time of 10 ns, followed by MM-GBSA binding energy calculations. The evaluated MM-GBSA binding energies for the investigated drugs with Mpro are listed in Fig. 2.

Fig. 2

Evaluated MM-GBSA binding energies for the investigated COVID-19 drug candidates as SARS-CoV-2 main protease (Mpro) inhibitors

What is interesting about the data in Fig. 2 that five out of the studied drugs demonstrated considerable binding energies (ΔGbinding <  − 40.0 kcal/mol), while the rest of the drugs were noticed with relatively weak binding energies ranged from − 10.5 to − 36.0 kcal/mol. The single most striking observation to emerge from the data comparison was the significant-high binding affinities of TMC-310911 and ritonavir towards SARS-CoV-2 Mpro with binding energies (ΔGbinding) of − 59.6 and − 54.4 kcal/mol, respectively.

Overall, the MM-GBSA//10 ns MD binding energies demonstrated the promising potentiality of TMC-310911, ritonavir, remdesivir, lopinavir, and darunavir to inhibit SARS-CoV-2 Mpro.

50 ns MD simulations

To increase the reliability of the predicted binding affinities, the top potent drugs (ΔGbinding <  − 40.0 kcal/mol with Mpro) were further subjected to longer MDs of 50 ns, and the corresponding MM-GBSA binding energies were consequently estimated (Fig. 2).

According to the calculated MM-GBSA//50 ns MD binding energies, TMC-310911 and ritonavir showed promising binding affinities against SARS-CoV-2 Mpro with ΔGbinding of − 52.8 and − 49.4 kcal/mol, respectively. The calculated values of Gdrug, Greceptor, and Gcomplex for the TMC-310911 and ritonavir in complex with Mpro are listed in Table S1. Compared to TMC-310911 and ritonavir, remdesivir showed a lower binding affinity with Mpro (ΔGbinding of − 42.7 kcal/mol). The estimated ΔGbinding over the 50 ns MD for lopinavir and darunavir with Mpro were observed to be higher than − 40.0 kcal/mol (ΔGbinding of − 35.6 and − 34.8 kcal/mol, respectively).

To reveal the nature of interactions of TMC-310911 and ritonavir with Mpro, MM-GBSA binding energy decomposition was performed over the 50 ns MD simulations. The evaluated energy components are listed in Table 3.

Table 3

MM-GBSA binding energies decomposition for TMC-310911 and ritonavir complexed with SARS-CoV-2 main protease (Mpro) over the 50 ns MD simulations

Drug name	Calculated MM-GBSA binding energy (kcal/mol)
	∆EVDW	∆Eele	∆EGB	∆ESUR	∆Ggas	∆GSolv	∆Gbinding
TMC-310911	− 65.8	− 23.1	44.2	− 8.1	− 88.8	36.1	− 52.8
Ritonavir	− 67.5	− 10.5	36.2	− 7.6	− 78.1	28.7	− 49.4

It is apparent from data in Table 3 that the Evwd interactions were predominant forces in the TMC-310911-Mpro and ritonavir-Mpro complexes with average values of − 65.8 and − 67.5 kcal/mol, respectively. As well, the Eele interactions were favorable in both TMC-310911-Mpro and ritonavir-Mpro complexes with average values of − 23.1 and − 10.5 kcal/mol, respectively.

Considering the promising binding affinities of TMC-310911 and ritonavir with Mpro, further investigations were carried out to inspect the drug-Mpro stability over the 50 ns MD simulation.

Post-Dynamics Analyses

Post-dynamics analyses were performed on the collected trajectories over the 50 ns MD simulations for TMC-310911-Mpro and ritonavir-Mpro complexes. The post-dynamics analyses involved root-mean-square deviation (RMSD), binding energy per-frame, center-of-mass (CoM) distance, and hydrogen bond length.

Binding Energy Per-Frame

To get a more in-depth insight into the stability of the drug inside the active site of Mpro, the correlation between MM-GBSA binding energy and time was investigated. The MM-GBSA binding energy per-frame for TMC-310911-Mpro and ritonavir-Mpro complexes is plotted versus time in Fig. 3.

Fig. 3

MM-GBSA binding energies vs. time for TMC-310911 (in black) and ritonavir (in red) towards SARS-CoV-2 main protease (Mpro)

As shown in Fig. 3, there was outright stability for TMC-310911-Mpro and ritonavir-Mpro over the generated MD trajectories with average values of − 52.8 and − 49.4 kcal/mol, respectively.

Hydrogen Bond Length and Center-of-Mass Distance

The stabilities of TMC-310911-Mpro and ritonavir-Mpro complexes were further investigated by measuring the hydrogen bond lengths and center-of-mass (CoM) distances between the drug and the proximal amino acid GLU166 residue over the 50 ns MD simulations (Fig. 4).

Fig. 4

a Hydrogen bond lengths and b center-of-mass (CoM) distances between TMC-310911 (in black) and ritonavir (in red) and the essential residue amino acid GLU166 of SARS-CoV-2 main protease (Mpro) over the 50 ns MD simulations

As can be concluded from data in Fig. 4a, TMC-310911 and ritonavir demonstrated high stabilities inside the active site of SARS-CoV-2 Mpro with average hydrogen bond lengths of 2.03 and 2.07 Å, respectively. Besides, the measured CoM distances were approximately constant around 7.2 and 10.1 Å over the 50 ns MD simulations for TMC-31091-Mpro and ritonavir-Mpro complexes, respectively (Fig. 4b).

Root-Mean-Square Deviation

The root-mean-square deviations (RMSD) throughout the 50 ns MD simulations were estimated to delineate the structural changes of drug-Mpro complexes. The RMSD's of backbone as a function of time with respect to the starting structure of TMC-310911- and ritonavir-Mpro complexes are plotted in Fig. 5.

Fig. 5

Root-mean-square-deviation (RMSD) of SARS-CoV-2 main protease (Mpro) backbone atoms from the initial structure bound with TMC-310911 (in black) and ritonavir (in red) over the 50 ns MD simulations

It can be seen from data in Fig. 5 that the RMSD values of TMC-310911-Mpro and ritonavir-Mpro stayed beneath 0.35 nm for the course of the MD simulations. The RMSD results indicated that TMC-310911 and ritonavir are tightly bonded in the active site of the SARS-CoV-2 Mpro, and they do not impact the overall topology of Mpro.

Overall, the post-dynamics analyses provided evidence on the high stability of TMC-310911-Mpro and ritonavir-Mpro complexes, demonstrating their potentiality as prospective anti-COVID-19 drugs.

Conclusion

Several clinical trials have recently been launched to evaluate the efficacy and safety of repurposed drugs to treat COVID-19. In the current study, the potencies of eighteen repurposed drugs in clinical development were evaluated against SARS-CoV-2 Mpro using combined molecular docking and molecular dynamics (MD) techniques. Molecular docking calculations revealed the high binding affinities of TMC-310911 and ritonavir towards Mpro with docking scores of − 8.3 and − 7.7 kcal/mol, respectively. The calculated MM-GBSA binding energies for TMC-310911 and ritonavir with Mpro over the 50 ns MD were − 52.8 and − 49.4 kcal/mol, respectively, demonstrating their high potencies as Mpro inhibitors. Post-dynamics analyses over the 50 ns MD confirmed the promising binding affinities and stabilities of TMC-310911 and ritonavir with SARS-CoV-2 Mpro. The present results shed light on TMC-310911 and ritonavir as prospective repurposed drugs for the treatment of COVID-19.

Below is the link to the electronic supplementary material.

Supplementary information 1 (DOC 14266 kb)