Old Arsenal to Combat New Enemy: Repurposing of Commercially available FDA Approved drugs against Main Protease of SARS-CoV2
In lack of vaccination and therapeutic drugs, the ongoing COVID-19 pandemic affected millions of people, causing 1,018,957 deaths worldwide (World health organization; 1 st October 2020). The conventional drug design pipeline for effective and safer drug development is a costly and time-intensive affair. It takes around ten years in general from identifying a clinical candidate to get the approvals for actual applications. An effective way to cut short drug design pipeline in such emergency cases could be the repurposing of already approved drugs against novel targets. Here in this work, we explored the structure-based drug screening approach to find potential inhibitors of SARS-CoV2 main protease (M pro ) from the library of already FDA approved commercially available drugs. The site-specific and blind docking studies, in combination, suggest three potential inhibitors of M pro , Ergotamine (ZINC000052955754), Nilotinib (ZINC000006716957) and Naldemedine (ZINC000100378061). Molecular dynamics (MD) simulations and binding free energy calculations using the MMPBSA method further reinforced the efficiency of the screened M pro inhibitor candidates. The work yields enough evidence to conduct rigorous experimental validation of these drugs before utilizing them for the therapeutic management of SARS-CoV2 infection.
Introduction<!>1-Protein and Ligand Preparation<!>2-Structure-based virtual screening<!>3-Blind docking studies<!>4-Inhibition constant calculation<!>5-Molecular dynamics (MD) simulation studies<!>∆𝐺 = ∆𝐻 -𝑇∆𝑆 = ∆𝐸𝑒l𝑒 + ∆𝐸𝑉DW + ∆𝐺𝑔b + ∆𝐺𝑛p -𝑇∆𝑆<!>… (2)<!>…(3)<!>1-Structure-based virtual screening results<!>2-Molecular Dynamic Simulation Analysis<!>3-Binding Free-energy Calculation<!>Conclusions
<p>As early as 31st December 2019, a cluster of pneumonia of unknown origin has been reported in Wuhan, Hubei Province, Republic of China. Later the causal organism is reported as a newly mutated strain of SARS Coronavirus (SARS-CoV) having higher pathogenicity and termed as SARS-CoV2 (Severe Acute Respiratory Syndrome CoronaVirus-2) causing the CoronaVirus Disease (COVID- 19) 1 . Among all the seven known coronaviruses, SARS-CoV2 is one of the three highly pathogenic human coronaviruses, including the MERS-CoV and SARS CoV. Although SARS-CoV2 exceeds the transmission rate from previously reported outbreaks of human coronaviruses, the virus has a much lower case fatality rate (CFR), i.e., 2-5% 34.4% and 10% reported in the case of MERS-CoV and SARS-CoV respectively [2][3][4] . The worldwide spread of the disease prompted the World Health Organization (WHO) to declare it a pandemic of a significant threat to public health, which needs to be addressed adequately by prevention, and therapeutics. As per WHO COVID-19 Dashboard, there were 31,798,308 confirmed cases of COVID-19 worldwide resulting in 973,653 deaths on 25th september 2020 5 . Clinical manifestation of SARS-CoV2 infection ranges in between being symptomatic to experience severe respiratory failure 6 . Early symptoms of the disease comprise anything or the combination of low-grade fever, fatigue, dry cough, myalgia, dyspnea, lack of taste, and smell 7 . Pneumonia is found to be associated with most of the SARS-CoV2 infections, but cases of having pleuritic chest pain also been reported 7 . Based on the severity of the symptoms, infected patients could be classified into mild, severe, and critical types 8 . However, most infections are self-limiting in nature, leaving approximately 15% cases to reach a severe level, which requires hospitalization and oxygen supplementation. Based on the reports, an additional 5% of such patients become critically ill, showing hypoxemia, acute respiratory distress, and multiple organ failure leading to death in the absence of necessary ventilator support for several weeks 7,[9][10][11][12][13] . Similar to other coronaviruses, SARS-CoV-2 is an enveloped spherical viral particles consisting of a 29903 bp (Accession no. NC_045512) extended positive-sense single-stranded RNA (ssRNA (+)) as genetic material. The genomic features of the SARS-CoV2 are wellconserved, showing 79.5% sequence similarity with early reported SARS-CoV. Despite having high sequence similarity, genetic modifications are highly impactful, resulting in resistance towards those drugs which were previously found effective against SARS-CoV 14 . Virus assembly consists of several protein molecules, viz-spike(S), envelope (E), membrane (M), nucleocapsid (N), major protease (Mpro), replicase-transcriptase, etc. Several of these proteins could serve as essential targets for therapeutic design and development, considering their crucial role in the viral life cycle [15][16][17] . The spike (S) protein of SARS-CoV2 is a trimeric class I viral fusion glycoprotein, which assists in an initial attachment to the angiotensin-converting enzyme II (ACE2) is known to allow fusion with the host cells such as alveolar epithelial cells 8,18,19 . The protein is an essential target for vaccine and inhibitor drug design 3,4 .</p><p>Chymotrypsin, like protease (3CLpro) or SARS-CoV2 major protease (Mpro), is also an attractive target for drug design and development due to high sequence conservation and its critical role viral pathogenesis 20 . The protease process the viral polyprotein 1a/1ab, into mature non-structural proteins (mNSPs), which play an essential role in blocking host immune response 21 . However, the primary function of SARS-CoV-2 Mpro is the regulation of replication and transcription machinery of SARS-CoV2 by controlling the synthesis of nonstructural proteins (NSPs) required for replicase-transcriptase assembly 22 . The protein also shows a high sequence similarity, i.e., 96% with SARS-CoV major protease (Mpro). It makes Mpro an ideal target for structure-based drug design and development. SARS-CoV2 Mpro is a homodimer of 306 amino acid long polypeptide chains. Crystal structure analysis suggests amino acid residues His41, Phe140, Leu141, Asn142, Gly143, Cys145, His163, Met165, Glu166, Leu167, Pro168, Phe185, Asp187, Gln189, Thr190, Ala191, and Gln192 define the substrate-binding pocket of protease 20,22,23 . At the same time, the His41 and Cys145 served as a part of the catalytic dyad. These residues could be regarded as prime targets during structurebased drug development to measure their interaction with potential ligand molecules to obtain SARS-CoV2 Mpro inhibitors [24][25][26] .</p><p>The repurposing of approved drugs provides an alternative approach to develop safe and effective therapeutics in case of rapidly emerging diseases and sudden outbreaks, as in the case of COVID-19 27 . While the conventional drug development pipeline is time-intensive and costs millions, repurposing already available drug molecules could result in the timely development of effective and safe drugs for clinical testing and applications 15,17,24 . In this work, we performed the structure-based virtual screening of the FDA-approved commercially available drugs to find high-affinity inhibitors of SARS-CoV-2 Mpro. First, we estimated the drugbinding affinities with SARS-CoV-2 Mpro and then conducted the interaction analysis to find better hits. Finally, based on the interaction analysis and drug properties, we identified two drugs, viz-Ergotamine and Naldemedine showing the appreciable binding affinity and specific interactions toward the SARS-CoV2 Mpro active site.</p><!><p>SARS-CoV2 crystal structure having PDB ID-5R83 has been studied extensively along with other crystal structures such as 5R7Y, 5R7Z, 5R80, 5R81, 5R84, 5RE4, 5RE9, 5REB, 5RF1, 5EFQ, 6LU7, and 6W63 to identify the critical residues involved in M pro -inhibitor interaction and understand the mechanism of inhibition. For molecular docking-based virtual screening, atomic coordinates of SARS-CoV2 M pro were obtained from the PDB ID 5R83. The crystal structure contained 304 amino acid residue long single chain of SARS-CoV2 M pro in complex with hydrolase inhibitor Z44592329 at 1.58 Å. Water molecules and hydrolase inhibitor were removed, and the geometry of the protein was cleaned through the "Prepare Protein" module of Biovia Discovery Studio 4.5(DS4.5). SDF file representing the library of 1045 FDA approved and commercially available drugs were downloaded from the ZINC15 database (https://zinc15.docking.org/). Energy minimization of all the 1045 ligand molecules was performed in vacuo, using 500 steps of conjugate gradients optimization algorithm with universal forcefield updated at one step interval. The energy minimization process was stopped at an energy difference of 0.1 kcal.mol -1 . Further, the energy minimized ligand molecules were converted to Autodock ligands (.pdbqt) format using PyRx 0.8 (https://pyrx.sourceforge.io/).</p><!><p>PyRx 0.8, along with AutoDock Vina v.1.1 28 , was used for structure-based docking studies of Mpro with FDA approved ligand molecules. The software was selected owing to its rapid yet efficient search capability to find accurate and active conformations 29 . While the SARS-CoV2 Mpro protein molecule was considered as a rigid structure, ligand drug molecules were considered flexible during the docking process. Molecular docking at the Active site of M pro was performed using AutoDock Vina (V1.1) for all the ligand drug molecules with a grid size of 41.67, 33.75, and 49.66 Å, centralized at -29.93, −39.34, and 14.58 Å for X, Y, and Z coordinates as mentioned in the table-1. The RMSD of 2 Å was found in the overlapped structure of the native inhibitor in the crystal form and the virtually docked form, as shown in Figure 1. The assessed binding energy score of the native hydrolase inhibitor (Z44592329) in the crystal structure (PDB ID: 5R83) was found -6.0 kcal.mol -1 . The threshold value for binding energy was set to -8.0 kcal.mol -1 to narrow down the search space to obtain only the highaffinity ligand molecules to SARS-CoV2 M pro as a target. A total of 9 out of 1045 compounds showed binding energy (DG˚) less than -8.0 kcal.mol -1 showing the high affinity of these ligand molecules towards M pro active site (Table -2).</p><!><p>Another round of structurally blind docking was performed with initially screened nine FDA approved drug ligand molecules using AutoDock Vina. The box dimensions for structurally blind docking were 57.68 64.00 and 48.89 Å while the M pro was centralized at -29.93, −39.34, and 14.58 Å for X, Y, and Z coordinates during the docking (Table -1). The grid spacing was set to 1.00 Å. The blind docking was performed to cover the entire protein structure under the grid search space to check whether the primary screened molecules have an affinity to bind to the enzyme other than the active site. Six out of the previously screened nine compounds showed comparable or high affinity to the sites other than the target's active site. Remaining three molecules, Ergotamine, Naldemedine, and Nilotinib, had shown a high affinity to the active site having minimal or no affinity to other sites. The docked poses of which are shown in Figure 2, and selected further for detailed analysis.</p><!><p>The inhibition constant selected ligand molecules has been calculated from the obtained ∆𝐺 parameter from molecular dynamics simulations and putting the values in the following formula-</p><p>The value is represented in terms of the negative decimal logarithm of inhibition constant</p><p>Where ∆𝐺 is binding affinity (kcal.mol -1 ), R (gas constant) is 1.98 cal. (mol.K) −1 , and T (room temperature) is 298.15 Kelvin.</p><!><p>Molecular dynamics (MD) simulation studies were conducted using LiGRO 30 , a GUI based MD simulation software with GROMACS 5.1.5 package 31 to evaluate the metastability of the native inhibitor, and screened drug ligands with SARS-CoV2 M pro . Transferable intermolecular potential water molecules (Tip3p) model for solvation along with Amber FF99sb forcefield was applied. The complex was placed inside a dodecahedron box (60 × 60 × 90 Å) containing 35,897 water molecules, such that the minimum distance between any atom of the proteinligand complex with the box wall was 2nm. The parameters of ligands were calculated using General Amber Force Field (GAFF) and AnteChamber PYthon Parser interfacE (ACPYPE) 32 module in LiGRO. The protein-ligand complex with an overall -3 charge was neutralized with 0.15M concentration of 107 Na+ and 103 Cl-ions. The equilibration steps under NVT and NPT ensemble were run for 1ns each in LiGRO, and the MD production run for each complex and free M pro was run under NPT ensemble at 1.01 bars at 310.15K for 100ns each at High performance cluster (HPC) facility at IIT Delhi. The particle mesh Ewald (PME) method and LINCS algorithm were applied to correct the electrostatic interactions and constrain the bonds with hydrogen atoms. The time step was 2fs, and the frames were stored after every 500 steps. The trajectories and images were visualized through VMD 33 and Chimera Software 34 .</p><p>The obtained data were analysed using GROMACS tools (rmsd, rmsf, and energy) to examine the overall stability of the protein-ligand system, fluctuations in the local residues, and structure over the complete simulation trajectory. In this study, MM-PBSA calculations were done using the g_MMPBSA 35 software, and the ΔG was calculated by equation (1),</p><!><p>(1) Where ΔEele and ΔEVDW are the total gas-phase energy/Molecular mechanics energy (ΔEMM), ΔGgb and ΔGnp refer to polar and non-polar solvation free energies respectively, and TΔS is conformational entropy upon binding. The solute and the solvent dielectric constant values were set to 80. The entropy was calculated as described in the methods reported by Wang et al. 36 given in equation 2 and 3,</p><!><p>Where BSASi = 4π (ri + rprob) 2</p><!><p>Here the N is the number of atoms in a molecule, wi, SASi, BSASi, ri denotes the weight, the solvent-accessible surface area, buried solvent accessible surface area, calculated using Equation 2 and the radius of the i atom respectively. The probe radius, rprob, was set to 0.08 nm. The term k is the adjustable parameter with the values taken from the same work.</p><p>In their approach, the conformational entropy of a molecule, S, can be obtained by summing up all atoms' contributions, no matter whether they are buried or exposed. The two types of surface areas SAS and BSAS, were weighted to estimate the contribution of an atom to S. Atoms having the same atom type shared the same weight. A general parameter k is applied to balance the contributions of the two types of surface areas. Since we obtained the enthalpic and entropic term, it was possible to obtain the complexes' binding free energies. The inner product of the residues present on the interface in the complex enabled us to match if the same residues are involved in interface formation in two different complexes. The values range from 0, which represents no match to 1 that represents the perfect match. Numerical value equal to or greater than 0.5 is usually considered a good match. In the present study, the inner product, γ, have been found between vectors Ac and A1, A2, and A3 respectively where, Ac = residue wise contribution to the interface area of the M pro in M pro -Z4492329 complex. A1 = residue wise contribution to the interface area of the M pro in M pro -Ergotamine complex. A2 = residue wise contribution to the interface area of the M pro in the M pro -Naldemedine complex. A3 = residue wise contribution to the interface area of the M pro in M pro -Nilotinib complex.</p><p>The ith element of vector A is the contribution of the ith residue of the protein to the interface of the complex. As a case study, a γ value of 0.67 suggests a 67 % match in the residues involved in the protein for two different complexes Ac and A1. γ = ( Ac) ⃗⋅(A1) ⃗ =0.67</p><!><p>FDA approved drug molecules were screened primarily based on their calculated binding energy score during docking studies. Three FDA approved drug molecules, viz-Ergotamine, Nilotinib, and Naldemedine were selected based on their binding affinities towards M pro active site. Among all docked drug molecules, the highest binding energy score was of Ergotamine, which was -9.5 kcal.mol -1 . The molecule is a α-1 adrenergic agonist vasoconstrictor, used for treating chronic migraine disorders. The molecule interacted within M pro active site, forming conventional carbon-hydrogen bonds with residues Glu166, His41, and Asn142. A π-Sulfur bond was evident with Met165 along with four alkyl and π-alkyl interactions with Met49 and Pro168 residues. A π-π T-shaped interaction was visualized with His41 residue, and several van der Waals interactions were evident in the docked complex (Figures -3B and 4B).</p><p>The next two drug molecules interacted with M pro of SARS-CoV2 were Nilotinib and Naldemedine having predicted binding energy score of -9.0 kcal.mol -1 and -8.9 kcal.mol -1 . Nilotinib is a tyrosine kinase inhibitor used in the second-line treatment of Chronic Myelogenous Leukaemia (CML). Unlike Ergotamine, the Nilotinib molecule in the docked complex has shown quite a different molecular interactions with SARS-CoV2 f active site residues. There were two conventional hydrogen bonds between Nilotinib and M pro within active site involving Glu166 and Asn119 residues. The docked assembly consisted of three halogen bonds in between Nilotinib molecules and amino acid residues Leu141, His163, and Met165. The three alkyl and π-alkyl bonds with Cys145 and His163 and multiple van der Waals interactions were evident in the docked assembly of M pro with Nilotinib molecule (Figures-3C and 4C). An unfavourable but peculiar donor-donor interaction with Glu143 was seen in the docked molecule. The μ-opioid receptor antagonist-Naldemedine, which is used to prevent opioid-induced constipation in adult patients, formed four in total (conventional and carbonhydrogen) hydrogen bonds with residues Thr26, Asn119 and Gly143. There were three alkyl and π-alkyl interactions found with the residues His41 and Met49. A π-π stacked interaction with Tyr118 and many van der Waals interactions were also evident in the docked assembly (Figures -3D and 4D).</p><p>In contrast to the top three screened FDA drugs, the native hydrolase inhibitor (Z44592329), PDB ID-5R83, was docked on the enzyme's active site. The binding energy was -6.0 kcal.mol -1 , much lower than the screened FDA approved inhibitors. The inhibitor formed five hydrogen bonds, including the conventional and non-conventional carbon-hydrogen bonds in total with the residues Phe140, His163, Met165, and Glu166. The inhibitor formed two alkyl and π-alkyl bonds with Met49 and Met165 residues while interacting with His41, His141, and Met165 through π-π T-shaped interactions. Few van der Waals interactions were also evident in the docked structure (Figures -3A and 4A). More van der Waals interactions were visible in drug docked structure than native hydrolase inhibitor docked assembly, which might be the reason for the higher binding energy score of drug molecules. Considering the substantial difference in the binding energy score and significant interactions between drug molecules and M pro active site, it seems possible for these FDA approved molecules to act as M pro inhibitors and abolish its enzymatic activity.</p><!><p>The stability of the drug molecule in the active site serves as an accurate quantifier to assess the fitness of potential inhibitors. M pro complex with all drug molecules and the native hydrolase inhibitor Z44592329 at the active site were subsequently subjected to molecular dynamic simulation. The protein backbone stability was assessed in terms of root-mean-square deviation (RMSD) values of the protein-drug complexes. The RMSD values were calculated concerning the initial frame of protein-drug complex and plotted, as shown in figure 5A. The RMSD values of M pro and M pro -ligand complexes increased gradually from 0 to 5 ns and remained stable after throughout the trajectory. The stability for residues was assessed in terms of root-mean-square fluctuation (RMSF) values of protein Cα backbone averaging over the 85-95ns interval, as shown in figure 5B. The RMSF values for the active site residues in M prodrug complexes ranged less than 0.7 Å; however, no significant change in RMSF values was observed as compared to that of free M pro except the N-terminal and C-terminal regions, which showed higher flexibility. Also, the mean inner product value was at least 0.5 or more for all the M pro -drug complexes indicating a good match for the residues involved in the M pro -Z44592329 complex, as shown in figure 6A.</p><!><p>The meta-stable region for all the complex, i.e., 85-95ns of the molecular dynamics simulation, was selected to calculate the binding free energy, as mentioned in the table-3. When compared with native inhibitor (Z44592329), Ergotamine and Naldemedine showed lower Gibbs free energy, which was in accordance as per the docking score obtained by PyRx. Nilotinib showed higher Gibbs free energy, which could be explained due to the lower conformational entropy with some fluctuations. The conformational entropy for Ergotamine, Naldemedine, and Z44592329 was constant, as shown in figure 6B. Since MM-PBSA is a computational method based on approximations and shows an intermediate level of accuracy between empirical scoring, and alchemical perturbation methods, there will always be a scope of improvement in the calculated values. The pKi (the negative logarithm of predicted inhibition constant) of the screened drug molecules, Ergotamine, Nilotinib, and Naldemedine, was 6.974, -1.036 and 5.823 micromolar (µM) respectively as calculated during the meta-stable region of the simulation (Table -3).</p><!><p>SARS-CoV2 M pro regulates the replication and transcription of the virus through the generation of non-structural proteins (NSPs) from viral polyprotein, which ultimately forms a replicasetranscriptase assembly. Hence, SARS-CoV2 M pro is an attractive drug target for structure-based drug design and development. Blocking of the M pro enzyme active site of SARS CoV2 will have prophylactic effects on viral infection and progression. Structure-based screening through docking of FDA approved drugs in the active site of SARS-CoV2 M pro helped to identify three potential inhibitors of main protease activity; viz-Ergotamine, Nilotinib, and Naldemedine. MD trajectory analysis and calculated binding energy (MMPBSA) verifies the stability of these drugs in the active site of the enzyme. The obtained result strengthens the idea-the binding of these drugs with SARS-CoV2 M pro could impede the essential proteolytic function of the enzyme causing inhibition of viral replication and disease progression. However, the identified drug molecules need to be tested through in vitro and clinical experiments for their efficacy before being utilized in prophylactic therapy and management of SARS-CoV2 infection.</p>