key: cord-0761229-64mlr1y4
authors: Ibrahim, Mahmoud A. A.; Abdelrahman, Alaa H. M.; Hegazy, Mohamed-Elamir F.
title: In-silico drug repurposing and molecular dynamics puzzled out potential SARS-CoV-2 main protease inhibitors
date: 2020-07-20
journal: J Biomol Struct Dyn
DOI: 10.1080/07391102.2020.1791958
sha: 7ac80c39d04ff3fe0a174fd0d2848bdd4a5a5e22
doc_id: 761229
cord_uid: 64mlr1y4

Herein, the DrugBank database which contains 10,036 approved and investigational drugs was explored deeply for potential drugs that target SARS-CoV-2 main protease (M(pro)). Filtration process of the database was conducted using three levels of accuracy for molecular docking calculations. The top 35 drugs with docking scores > −11.0 kcal/mol were then subjected to 10 ns molecular dynamics (MD) simulations followed by molecular mechanics–generalized Born surface area (MM-GBSA) binding energy calculations. The results showed that DB02388 and Cobicistat (DB09065) exhibited potential binding affinities towards M(pro) over 100 ns MD simulations, with binding energy values of −49.67 and −46.60 kcal/mol, respectively. Binding energy and structural analyses demonstrated the higher stability of DB02388 over Cobicistat. The potency of DB02388 and Cobicistat is attributed to their abilities to form several hydrogen bonds with the essential amino acids inside the active site of M(pro). Compared to DB02388 and Cobicistat, Darunavir showed a much lower binding affinity of −34.83 kcal/mol. The present study highlights the potentiality of DB02388 and Cobicistat as anti-COVID-19 drugs for clinical trials. Communicated by Ramaswamy H. Sarma.

The case of a recent outbreak of coronavirus infectious disease has been identified and attributed to the Severe Acute Respiratory Syndrome Coronavirus 2 (SARS-CoV-2, namely 2019nCOV), which originated from the meat market in Wuhan, China (Cao et al., 2020) . The World Health Organization (WHO) has announced COVID-19 as a potential global health threat because of high mortality rate, high basic reproduction number (R0), and lack of clinically approved drugs and vaccines for this disease. Due to the increasing number of deaths all over the world caused by COVID-19, there is an urgent demand for potent anti-COVID-19 drugs. Traditional de novo drug discovery is a time-consuming and expensive process. Drug repurposing is an alternative time-efficient and cost-effective strategy, in which new uses for approved or investigational drugs are identified (Elmezayen et al., 2020; Pushpakom et al., 2019) . Experimental and computational drug repurposing have attracted interests from both the pharmaceutical industry and the research community over the last decades Pushpakom et al., 2019) . Several successful attempts of drug repurposing have been reported for treatment of various diseases (Karaman & Sippl, 2019; Yella et al., 2018) .

All recent studies showed that the usage of therapeutic agents as anti-viral (Lai et al., 2020) , anti-malarial (Gao et al., 2020) , and herbal treatments (Luo et al., 2020) to treat COVID-19 infection is unclear. The re-investigation of available drugs is urgently needed for COVID-19 treatment and saving of lives all over the world.

Coronavirus polyproteins include two main proteases, main protease (3C-like protease, M pro ) and a papain-like protease (PL pro ) (Hilgenfeld, 2014) . Both coronaviral proteases play a vital role in the replication and transcription of the virus and, in turn, are considered as charming anti-viral targets for the discovery of novel drugs (Xue et al., 2008; Zhao et al., 2013) . The first crystal structure of SARS-CoV-2 M pro with a peptidomimetic inhibitor (N3) has been recently reported (Jin et al., 2020) , opening a new horizon for the utilization of structure-based drug discovery approaches towards the inhibition of SARS-CoV-2 M pro .

The aim of this research project has therefore been to evaluate the reuse of approved and clinical trial drugs as anti-COVID-19 drugs by combining molecular docking and molecular dynamics simulation methods. 10,036 drug candidates in the DrugBank database were first prepared and their binding affinities were filtered based on their docking scores against SARS-CoV-2 M pro . Three levels of docking calculations were utilized towards verification of docking scores accuracy. The top 35 drugs with highest docking scores with M pro were then subjected to molecular dynamics simulation of 10 ns. Molecular mechanics-generalized Born surface area (MM-GBSA) approach was employed to estimate drug-M pro binding energies (DG binding ) over the simulated time. The stability of the most potent drugs inside the active site of M pro was further investigated over 100 ns molecular dynamics simulations. The study presented sheds light on DB02388 and Cobicistat which should be clinically investigated as anti-COVID-19 drugs.

Prior to molecular docking calculations, the DrugBank database containing 10,036 approved and investigational drugs was prepared (Wishart et al., 2006 (Wishart et al., , 2018 . The database was retrieved in SDF format and the 3D structures of drug molecules were then generated using Omega2 software (Hawkins et al., 2010; OMEGA, 2013) . 3D chemical structures of the top potent eight DrugBank drugs showing high MM-GBSA binding energies with SARS-CoV-2 M pro are illustrated in Figure 1 . The geometrical structures of all drugs were minimized with MMFF94S force field using SZYBKI software (SZYBKI, 2016) and the partial atomic charges were assigned using Gasteiger method (Gasteiger & Marsili, 1980) . Only neutral drug molecules containing C, H, N, O, S, P, F, Cl, Br and I were considered in the current study. The prepared files of DrugBank database may be accessed via www.compchem. net/ccdb.

The crystal structure of SARS-CoV-2 main protease (M pro ) complexed with peptidomimetic inhibitor (N3) (PDB code: 6LU7 (Jin et al., 2020) ) was taken as a template for all molecular docking and molecular dynamics calculations. Water molecules and ions were stripped out. The protonation state of M pro was investigated using Hþþ server, and all missing hydrogen atoms were added (Gordon et al., 2005) . For molecular docking calculations, pdbqt file for M pro was prepared according to AutoDock protocol (Forli et al., 2016) .

All molecular docking calculations were carried out using Autodock4.2 software (Morris et al., 2009) . Three different levels of docking accuracy were conducted, namely: standard, moderate and expensive. In the three calculations, all docking parameters were kept to default values, except the number of genetic algorithm (GA) run and the maximum number of energy evaluation (eval). The latter two variables were set to 50, 100 and 250 and 5,000,000, 10,000,000 and 25,000,000 for standard, moderate and expensive molecular docking calculations, respectively. A schematic representation of the utilized calculations is depicted in Figure 2 . The docking grid was set to 60A˚x 60A˚x 60A˚with a grid spacing value of 0.375 A˚, and the grid center was placed at the center of the M pro active site.

Molecular dynamics (MD) simulations were carried out on top potent drugs in complex with M pro using AMBER16 software (Case et al., 2016) . General AMBER force field (GAFF) (Wang et al., 2004) and AMBER force field 14SB (Maier et al., 2015) were used to describe the studied drugs and M pro , respectively. The atomic partial charges of the studied drugs were calculated using the restrained electrostatic potential (RESP) approach (Bayly et al., 1993) at the HF/6-31G Ã level with the help of Gaussian09 software (Frisch et al., 2009 ). The docked drug-M pro complexes were solvated in a cubic water box with 15 Å distances between the edges of the box and any atom of drug or drug-M pro complexes. The solvated drug-M pro complexes were then minimized for 5000 step and subsequently gently annealed from 0 K to 300 K throughout 50 ps. The systems were equilibrated for 1 ns and production stages were conducted over simulation times of 10, 50 and 100 ns. Periodic boundary conditions and the NPT ensemble were adopted in all MD simulations, including both the equilibration and production stages. Long-range electrostatic interactions under periodic conditions with a direct space cut-off 12 Å was treated with Particle Mesh Ewald (PME) method (Darden et al., 1993) . Langevin dynamics with the collision frequency gamma_ln set to 1.0 was used to keep the temperature constant at 298 K. Berendsen barostat with a relaxation time of 2 ps was employed to control the pressure (Berendsen et al., 1984) . A time step of 2 fs and the SHAKE option to constrain all bonds involving hydrogen atoms were used. Coordinates and energy values were collected every 10 ps over the production stage for binding energy calculations and post-dynamics analyses. All MD simulations were performed with the GPU version of pmemd (pmemd.cuda) in AMBER16 on CompChem GPU/CPU cluster (hpc.compchem.net). Visual Molecular Dynamics (VMD) (Humphrey et al., 1996) and LigPlotþ (Laskowski & Swindells, 2011) software were used to visualize drug-M pro complexes in 3D and 2D, respectively.

Molecular mechanics-generalized Born surface area (MM-GBSA) (Massova & Kollman, 2000) approach implemented in AMBER16 software was used to evaluate the binding energy of the studied drugs with M pro . The MM-GBSA binding free energies were estimated as follows:

where the energy term (G) is estimated as:

with E vdw , E ele , G GB and G SA as the van der Waals, electrostatic, General Born solvation, and surface area energies, respectively. For the drugs, entropy contributions were neglected.

Molecular docking is utilized as a substantial tool in the drug discovery process (Ferreira et al., 2015; Guedes et al., 2014; Meng et al., 2011) . Towards drug repurposing for discovering potent M pro inhibitor, the DrugBank database was deeply screened (Wishart et al., 2006 (Wishart et al., , 2018 . The database which contains 10,036 approved and investigational drugs was first retrieved and prepared as described in Computational Methodology. Due to a large number of investigated drugs, three levels of docking calculations were performed to reduce computing time.

First, standard molecular docking calculations with GA ¼ 50 and eval ¼ 2,500,000 were conducted for 10,036 DrugBank molecules against M pro and their docking scores were evaluated. Drugs were ranked based on the calculated docking scores, and the top 1,000 potent ones were subjected to further calculations. The top 1,000 potent drugs identified from standard docking calculations were redocked against M pro using GA and eval values of 100 and 10,000,000, respectively. Based on their binding scores, the drugs were ranked and the top 250 drugs were selected. From the latter calculations, the top potent 250 drugs were finally re-docked using expensive docking parameters of GA ¼ 250 and eval ¼ 25,000,000. The predicted expensive docking scores for the top 250 potent drugs with M pro are listed in Table S1 .

As seen from data listed in Table S1 , 35 drugs showed docking scores > À11.0 kcal/mol. Accordingly, these 35 drugs were subjected to further investigation using MD techniques. The binding affinities and features for these 35 drugs are summarized in Table S2 . The binding modes of eight of the most promising drugs, selected based on later MM-GBSA calculations, inside the active site of M pro are depicted in Figure  3 , and their binding features are listed in Table 1 .

According to data given in Table S2 and Figure 3 , most of the drugs share the same binding modes, forming an essential hydrogen bond with GLU166 and other hydrogen bonds with various amino acid residues inside the active site. For instance, the high docking score of À11.74 kcal/mol for DB02388 with M pro is attributed to its ability to form three 

short hydrogen bonds with GLU166, TYR54 and ASP187 with bond lengths of 2.10, 2.19 and 2.11 Å, respectively.

Although molecular docking is a fast and widely-used technique, reliability of its predicted drug-target binding affinity is still under debate (Pagadala et al., 2017; Pantsar & Poso, 2018) . To achieve reliable predictions for binding affinities, drug-receptor intermolecular interactions and conformational flexibilities of drug-receptor complexes, solvent effects and dynamics must be considered. This can be accomplished using all-atom molecular dynamics (MD) simulations combined with binding energy calculations over reasonable simulation time. Considering computational costs and time, molecular mechanics combined with Generalized Born and surface area continuum solvation (MM-GBSA) methods is a popular approach for the estimation of drug-target binding energies.

Therefore, the 35 identified drugs retrieved from expensive molecular docking calculations were further investigated using MD techniques. The docked drug-M pro complexes were prepared and subjected to 10 ns molecular dynamics simulations and the corresponding binding energies were evaluated using MM-GBSA approach. The estimated average MM-GBSA binding energies over the 10 ns MD simulations are listed in Table S3 . According to the calculated MM-GBSA binding energies (DG binding ), only eight out of the investigated 35 drugs showed considerable binding energies (DG binding > À40.0 kcal/mol); while 9 and 18 drugs were observed with relatively weak binding energies in ranges of À40.0! DG binding > À30.0 and À30.0! DG binding > À23.0 kcal/mol, respectively. Interestingly, Cobicistat (DB09065) and DB02388 exhibited promising binding affinities towards M pro with binding energies (DG binding ) of À59.84 and À50.72 kcal/mol, respectively. Measuring the essential hydrogen bond length formed between the drug and the key amino acid GLU166 would reflect the drug-M pro binding affinity. Therefore, the average hydrogen bond length for the eight drugs (with DG binding > À40.0 kcal/mol) were calculated and listed in Table 2 . According to the average hydrogen bond length, the eight drugs form a strong hydrogen bond with GLU166 with bond length ranging from 2.17 to 2.86 Å, except for MK-6186 which forms hydrogen bond with HIS163 with a length of 2.33 Å.

In an attempt to achieve reliable results and spot-light the most potent drug, the 8 drugs with DG binding > À40.0 kcal/mol were further investigated over 50 ns molecular dynamics simulations, and their corresponding Table 2 .

According to the results presented in Table 2 , there was no significant difference between MM-GBSA binding energies estimated over 50 ns and those over 10 ns, except in the case of Cobicistat. For the latter drug, the DG binding decreased from À59.84 to À53.34 kcal/mol over 10 ns and 50 ns, respectively. Inspecting the hydrogen bond length between Cobicistat and GLU166 revealed an increase in the bond length (from 2.17 Å to 3.01 Å for 10 ns and 50 ns MD simulations, respectively). This finding raises intriguing questions regarding the stability of Cobicistat inside the active site of M pro . Interestingly, DB02388 maintained its binding affinity with M pro over the 50 ns MD simulations, with an average hydrogen bond length of 2.31 Å with GLU166.

Molecular docking and molecular dynamics results demonstrated promising binding affinity of DB02388 and Cobicistat (DB09065) for M pro . DB02388, the imidazole-pyrimidine compound, is an investigational drug developed for the treatment of neurological disorders via inhibition mitogenactivated protein kinase 10 (MAPK10, known as c-Jun N-terminal kinase 3(JNK3)) (Scapin et al., 2003; Wishart et al., 2006 Wishart et al., , 2018 . Cobicistat, marketed under the name Tybost (formerly GS-9350), acts as a pharmacokinetic enhancer by inhibiting cytochrome P450 3A isoforms (CYP3A) for treating the infection with human immunodeficiency virus (HIV) (Deeks, 2014; Wishart et al., 2006 Wishart et al., , 2018 . It is worth mentioning that Cobicistat, which was identified and highlighted in the current manuscript, has undergone in two clinical trials (ChiCTR2000029544 trial started on the third of February 2020 and ChiCTR2000029548 trial started on the fourth of February 2020) (Harrison, 2020) . The trials' results are not yet available to us. Therefore, a comparison between the affinity and stability of DB02388 and Cobicistat was performed over 100 ns MD simulations. Snapshots for DB02388-and Cobicistat-M pro complexes were collected every 10 ps over the 100 ns MD simulations, giving 10,000 snapshots in total. Based on the extracted snapshots, several energetic and structural analyses were performed, including binding energy per frame, hydrogen bond length, center-of-mass (CoM) distance and rootmean-square-deviation (RMSD).

Correlation between the binding-energy and time would give a deeper insight into the stability of drug-M pro interactions. Therefore, MM-GBSA binding energy was calculated per frame for each of DB02388-and Cobicistat-M pro complexes and plotted vs time in Figure 4 .

According to data in Figure 4 , DB02388 showed narrowfluctuated binding energies with an average value of À49.67 kcal/mol over the 100 ns MD simulations. However, higher fluctuations in the binding energy for Cobicistat with Figure 6 . Root-mean-square-deviation (RMSD) of the backbone atoms from the initial structure for Cobicistat (in black) and DB02388 (in red) with the SARS-CoV-2 main protease (M pro ) through 100 ns MD simulations. M pro were noticed, ranging from À22.12 to À72.23 kcal/mol with an average value of À46.60 kcal/mol. These findings indicated that DB02388 bounds more tightly with M pro complex than Cobicistat. The decomposition of MM-GBSA binding energy was also performed to reveal the nature of dominant interactions in DB02388-and Cobicistat-M pro complexes (Table 3) . Energy decomposition results showed that E vdw is the dominant force in the drug-M pro binding affinity, À63.78 and À61.88 kcal/mol for DB02388-M pro and Cobicistat-M pro complexes, respectively. Besides, E ele interaction of DB02388 and Cobicistat with M pro is favorable (À10.93 and À12.85 kcal/mol, respectively). The hydrogen bond lengths between the DB02388 and Cobicistat and H-N of the key amino acid GLU166 residue were measured through 100 ns MD simulations and depicted in Figure 5(a) . As shown in Figure 5 (a), DB02388 displayed higher stability inside the active site of M pro with an average hydrogen bond length of 2.32 Å. However, Cobicistat showed stable hydrogen bond through the first 33 ns with an average length of 1.95 Å. After 33 ns MD simulations, there was a dramatic increase in the hydrogen bond length, giving an average bond length of 5.35 Å throughout the rest 67 ns MD time. This reflects the downturn of Cobicistat-M pro binding energy from À59.84 kcal/mol over the first 10 ns to À46.60 kcal/mol over 100 ns MD simulation.

The stability of DB02388-and Coblistat-M pro complexes was evaluated through measuring the center-of-mass (CoM) distance between the drug and that of GLU166 amino acid residue over the simulated 100 ns (Figure 5(b) ). According to Figure 5 (b), the average CoM distance between DB02388 and GLU166 was nearly constant around 6.3 Å throughout the 100 ns MD simulations, while the corresponding distance for Cobicistat was much higher indicating instability of the latter complex.

Structural changes in the M pro were analyzed using rootmean-square deviation (RMSD) for the backbone atoms of the complexes relative to the starting structures. The RMSD for DB02388 and Cobicistat-M pro over 100 ns MD simulations were calculated and plotted in Figure 6 (a). As shown in Figure 6 (a), the backbone of DB02388-M pro complex was well stabilized after the first 5 ns period, giving RMSD with less than 0.05 nm over the rest 95 ns period. This confirms that DB02388 is tightly bonded in the active site of the M pro and does not affect the overall topology of the protein. On the other hand, Cobicistat-M pro complex showed considerable stability in the first 45 ns period, followed by a dramatic instability for the rest 55 ns, with RMSD 0.05 and 0.17 nm, respectively. The RMSD results are appropriate with lower stability of Cobicistat-M pro compared with DB02388-M pro complex.

In summary, the energetic and structural analyses showed the high stability of DB02388-M pro complex over Cobicistat-M pro complex.

Darunavir (DB01264) is a HIV protease inhibitor and has been recently subjected to clinical trial as COVID-19 drug (Harrison, 2020) . Therefore, expensive molecular docking and MD simulations followed by MM-GBSA binding energy calculations were conducted for Darunavir with M pro . The docking score and MM-GBSA binding energy over 100 ns MD simulations were estimated and compared to that of DB02388 (Table 4) .

Interestingly, the predicted binding mode of Darunavir with M pro showed the ability of Darunavir to form three short hydrogen bonds with GLU166 and LEU167, compared to DB02388 which forms three hydrogen bonds (Figure 7) . Despite its larger number of formed hydrogen bonds with active site's amino acids, Darunavir-M pro docking score was lower than of DB02388 (À8.19 and À11.74 kcal/mol, respectively). MM-GBSA binding energy decomposition of Darunavir-M pro complex revealed that van der Waals (DE VDW ) interaction between Darunavir and M pro was only À47.37 kcal/ mol, compared to À63.78 kcal/mol for DB02388, giving total binding energy (DG binding ) of À34.83 and À49.67 kcal/mol, respectively.

SARS-CoV-2 main protease (M pro ) is one of the most charming targets for the discovery of potent anti-COVID-19 drugs. In the presented study, in-silico drug repurposing technique was utilized to identify potent M pro inhibitors. Filtration of DrugBank database was performed using molecular docking and molecular dynamics (MD) followed by molecular mechanics-generalized Born surface area (MM-GBSA) binding energy calculations. According to molecular docking calculations, 35 drugs showed docking scores higher than À11.0 kcal/mol with M pro . Investigation of binding features revealed the ability of the identified drugs to form essential hydrogen bond with the key amino acid GLU166 besides other hydrogen bonds with active site's amino acids. MM-GBSA binding energy calculations through MD simulations over 50 ns revealed that DB02388 and Cobicistat exhibited promising binding affinities of À49.67 and À46.60 kcal/mol, respectively, with M pro . Energetic and structural analyses throughout 100 ns MD revealed the promising binding affinity and stability of DB02388 with M pro over Cobicistat. Therefore, tenable conclusions from the study are i) DB02388 forms stable hydrogen bond with GLU166 with an average bond length of 2.32 Å, ii) center-of-mass distance between DB02388 and GLU166 was nearly constant over the 100 ns MD simulations, iii) narrow-fluctuated root-mean-square-deviation (RMSD) below 0.4 nm was observed over the 100 ns MD. The current results establish that DB02388 and Cobicistat hold promise as inhibitors against SARS-CoV-2 M pro and are ready for in vitro inhibition against COVID-19.

No potential conflict of interest was reported by the author(s).

The computational work was completed with resources supported by the Science and Technology Development Fund, STDF, Egypt, Grants No. 5480 & 7972. ORCID Mahmoud A. A. Ibrahim http://orcid.org/0000-0003-4819-2040 Mohamed-Elamir F. Hegazy http://orcid.org/0000-0002-0343-4969

A wellbehaved electrostatic potential based method using charge restraints for deriving atomic charges: The RESP model

Molecular-dynamics with coupling to an external bath

A trial of Lopinavir-Ritonavir in adults hospitalized with severe Covid-19

Particle mesh Ewald: AnNÁlog(N) method for Ewald sums in large systems

Cobicistat: A review of its use as a pharmacokinetic enhancer of Atazanavir and Darunavir in patients with HIV-1 infection

Drug repurposing for coronavirus (COVID-19): In silico screening of known drugs against coronavirus 3CL hydrolase and protease enzymes

Molecular docking and structure-based drug design strategies

Computational protein-ligand docking and virtual drug screening with the AutoDock suite

Breakthrough: Chloroquine phosphate has shown apparent efficacy in treatment of COVID-19 associated pneumonia in clinical studies

Iterative partial equalization of orbital electronegativity -a rapid access to atomic charges

Hþþ: A server for estimating pKas and adding missing hydrogens to macromolecules

Receptor-ligand molecular docking

Coronavirus puts drug repurposing on the fast track

Conformer generation with OMEGA: Algorithm and validation using high quality structures from the protein databank and Cambridge structural database

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

VMD: Visual molecular dynamics

Structure of Mpro from SARS-CoV-2 and discovery of its inhibitors

Computational drug repurposing: Current trends

Targeting SARS-CoV-2: A systematic drug repurposing approach to identify promising inhibitors against 3C-like proteinase and 2 0 -O-ribose methyltransferase

Identification of chymotrypsin-like protease inhibitors of SARS-CoV-2 via integrated computational approach

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): The epidemic and the challenges

LigPlotþ: multiple ligand-protein interaction diagrams for drug discovery

Can Chinese medicine be used for prevention of Corona Virus disease 2019 (COVID-19)? A review of historical classics, research evidence and current prevention programs

ff14SB: Improving the accuracy of protein side chain and backbone parameters from ff99SB

Combined molecular mechanical and continuum solvent approach (MM-PBSA/GBSA) to predict ligand binding

Molecular docking: A powerful approach for structure-based drug discovery

AutoDock4 and AutoDockTools4: Automated docking with selective receptor flexibility

OMEGA (Version 2.5.1.4)

Software for molecular docking: A review

Binding affinity via docking: Fact and fiction

Drug repurposing: Progress, challenges and recommendations

The structure of JNK3 in complex with small molecule inhibitors: Structural basis for potency and selectivity

SZYBKI (Version 1.9.0.3)

Development and testing of a general amber force field

DrugBank 5.0: A major update to the DrugBank database

DrugBank: a comprehensive resource for in silico drug discovery and exploration

Structures of two coronavirus main proteases: implications for substrate binding and antiviral drug design

Changing trends in computational drug repositioning

Recent developments on coronavirus main protease/3C like protease inhibitors