Discovery of New Hydroxyethylamine Analogs Against 3CL pro Protein Target of SARS-CoV-2: Molecular Docking, Molecular Dynamics Simulation and Structure-Activity Relationship Studies
A novel coronavirus, SARS-CoV-2 has caused a recent pandemic called COVID-19 and a severe health threat around the world. In the current situation, the virus is rapidly spreading worldwide, and the discovery of vaccine and potential therapeutics are critically essential. The crystal structure for main protease (M pro ) of SARS-CoV-2, 3-chymotrypsin-like cysteine protease (3CL pro ) was recently made available and is considerably similar to previously reported SARS-CoV. Due to its essentiality in viral replication, it represents a potential drug target. Herein, computer-aided drug design (CADD) approach was implemented for the initial screening of 13 approved antiviral drugs. Molecular docking of 13 antivirals against 3chymotrypsin-like cysteine protease (3CL pro ) enzyme was accomplished and indinavir was described as a lead drug with a docking score of -8.824 and a XP Gscore of -9.466 kcal/mol. Indinavir possesses an important pharmacophore, hydroxyethylamine (HEA), and thus a new library of HEA compounds (>2500) was subjected to virtual screening that led to 25 hits with a docking score more than indinavir. Exclusively, compound 16 with docking score of -8.955 adhered to drug like parameters, and the Structure-Activity Relationship (SAR) analysis was demonstrated to highlight the importance of chemical scaffolds therein. Molecular Dynamics (MD) simulation studies carried out at 100ns supported the stability of 16 within the binding pocket. Largly, our results supported that this novel compound 16 binds to the domain I & II, and domain II-III linker of 3CL pro protein, suggesting its suitablity as strong candidate for therapeutic discovery against COVID-19. Lead compound 16 could pave incredible directions for the design of novel 3CL pro inhibitors and ultimately therapeutics against COVID-19 disease.
Introduction.<!>Material and Method.<!>Preparation of Protease structure and active site identification.<!>Preparation of antiviral, other approved drugs and analogs library for virtually screening.<!>Molecular docking of chemical libraries against target protease active site.<!>ADME Calculations.<!>Molecular Dynamics (MD) Simulation.<!>Result and Discussion.<!>Virtual screening of antiviral drugs and new HEA analogs against SARS-CoV-2 protein 3CL pro<!>Molecular Dynamics (MD) Simulation of indinavir, compound 16.<!>Conclusion.<!>*Corresponding Author
<p>Coronaviruses (CoV) are a large group of viruses consisting of a core of genetic material enveloped with protein spike appearing like a crown, which means "corona" in Latin. 1 Diverse variety of coronaviruses is known, which causes mild respiratory diseases and sometimes gastrointestinal symptoms in different animal species. In humans, four CoV (229E, NL63, OC43 and HC HKU1) are endemic and cause respiratory diseases that can range from a common cold to lung failure, but the disease remains mild in most of the cases. 2,3 However, other types of CoV i.e., SARS-CoV (Severe Acute Respiratory Syndrome), and MERS-CoV (Middle East Respiratory Syndrome) can cause severe respiratory diseases as identified in China (2002) and Saudi Arabia (2012), respectively. 4 These viruses are bat-borne in nature and circulate in a range of animals and sometimes transmitted from animals to humans (i.e., SARS-CoV transmitted to humans from civet cats 5 and MERS-CoV from dromedary camels). 6 In December 2019, a cluster of pneumonia cases were observed in a group of people associated with seafood and animal market in China. 7 Subsequently, the outbreak was attributed to a novel CoV related to the SARS virus based on the genetic similarities. 8 Later, the disease was named as COVID-19 by World Health Organisation (WHO) caused by the virus SARS-CoV-2. 9 The outbreak started in Wuhan, China and rapidly escalated to other countries. The USA, Italy, Spain and China lead with the highest number of COVID-19 cases officially reported. COVID-19 caused approximately 2,924 deaths, and 85,403 confirmed cases were identified with a mortality ratio of 3.42% till the end of February 2020 across the countries. 10 The number of cases increased suddenly, and the disease was declared a pandemic by WHO on 11 March 2020. 11 As of 30 March 2020, a total of 745,308 confirmed cases of COVID-19 have been reported in 177 countries, resulting in 35,307 deaths. 12 A sharp increase in the number of confirmed cases and the number of deaths due to COVID-19 can be clearly observed across the world during the month of March, 2020 as shown in Fig. 1. 13 Major symptoms of this disease includes high fever, and respiratory symptoms (i.e. cough and shortness of breath) whereas in severe cases pneumonia and kidney failure are a major cause for the death. 14 Although, aseries of rapid and specific diagnostic tools are available for COVID-19 disease 15 , vaccines and SARS-CoV-2 specific therapeutic treatments are not available. Recently, the 3D crystal structure of SARS-CoV-2 main protease, M pro , also called 3CL pro , complexed with an inhibitor N3 was released. 16 3CL pro enzyme of SARS-CoV-2 processes polyproteins by proteolytic action of replicase enzyme (pp1a and pp1ab) to release the functional polypeptide (Fig. 2B). It is a dimeric protein, which contains two asymmetric units designated as protomers. Each protomer is composed of three domains, namely domain I (residues 8-101), domain II (residues 102-184), and domain III (residues 201-303). Domain III contains five α-helices and connected with domain II through an extended loop region (residues 185-200) (Fig. 2A). The 3CL pro has a Cys_145 and His_41 catalytic dyad, and the substrate-binding site is located in a cleft between domain I and II. These features are similar to previously reported protein structures of SARS-CoV. [17][18][19][20][21][22][23] The 3CL pro of SARS-CoV-2 is an attractive drug target, and CADD is considered as an undisputable and significant approach to discover antiviral drug candidates. Identically, studies have been adopting a fast identification of drug candidates against SARS-CoV-2 through virtually screening. 33,34 Antiviral drugs such as HIV protease inhibitors, 24 approved drugs such as ribavirin, oseltamivir and few more have been studied against earlier reported for SARS-CoV. 25 In this paper, we aim to perform the fast discovery of the potential drug candidates against 3CL pro of SARS-CoV-2 through virtual screening of approved antiviral therapeutics and build a focused library of new potential compounds. Our screening results indicated that indinavir, a known inhibitor of HIV protease, has high potential against main protease, 3CL pro of SARS-CoV-2 based on docking scores. The presence of a high-valued pharmacophore, hydroxyethylamine (HEA) in indinavir remains an incredible core moiety for the design of new antiviral drugs. Encouraged with this findings, we designed a new library of compounds possessing this crucial backbone, HEA, and followed extensive in silico studies towards the discovery of potential drugs candidates against SARS-CoV-2.</p><!><p>The protease structure, 3-chymotrypsin-like cysteine protease (3CL pro ) enzyme of SARS-CoV-2 (PDB ID: 6LU7) with 2.1Å was obtained from the RCSB site (http://www.rcsb.org). The computational work was performed using Schrodinger software.</p><!><p>The SARS-CoV-2 virus protein structure was prepared in the protein preparation wizard and prime module of Schrodinger suite to remove defects such as missing hydrogen atoms, incorrect bond order assignments, charge states, orientations of various groups and missing side chains suite. 26,27 Removal of steric clashes and strained bonds/angles were achieved through restrained energy minimization, allowing movement in heavy atoms up to 0.3 Å.</p><!><p>The two libraries were created, firstly for antiviral and other approved drugs such as indinavir, saquinavir, nelfinavir, ritonavir, lopinavir, atazanavir, amprenavir, darunavir, nelfinavir, oseltamivir, tipranavir, fosamprenavir, and galidesvir to select the top-ranked approved drug as reference molecule or positive control. Another library was created from HEA based analogs to compare the efficieny of these compounds as potent analogs. The coordinate files of approved drugs were downloaded from Protein Data Bank at RCSB site (http://www.rcsb.org) and PubChem. 28,29 Over thousand analogs were prepared by the maestro tool. The structure of both approved drugs and analogs were prepared prior to docking by Ligprep. The main objective of Ligprep was to take 2D or 3D structures and produce the corresponding lowenergy 3D structures for use by programs such as Glide. All the parameters were kept default except chirality parameters for antivirals and analogs. The chirality was kept default from 3D structure and all combinations of chirality were developed for known antivirals and analogs.</p><p>The next steps included desalting, generation of tautomer, and possible ionization states at pH 7 ± 2. 30 Schrödinger suite inbuilt Epik module was used to predict the ionization states of all compounds. 31</p><!><p>Site specific molecular docking of approved drug (Table 1) and HEA based analogs (Table 2) against protease was performed using the Glide module of Schrödinger suite. Glide tool was used for receptor grid preparation at default parameters. 32 The screening of both libraries was performed by Glide at extra precision (XP). The purposes of the XP method are to weed out false positives and to provide a better correlation between excellent poses and good scores.</p><!><p>Swiss ADME 33 and PKCSM 34 were used to calculate ADMET (i.e. Absorption, Distribution, Metabolism, Excretion, and Toxicity) profile of identified hit compounds (1-25). The predicted ADMET properties includes Molecular weight, H-bond acceptor, H-bond donors, Predicted octanol/water partition coefficient (MLogP), TPSA (Total Polar Surface Area), Lipinski (druglikeness), Rat LD50, and Hepatotoxicity. The results are shown in Table 3.</p><!><p>Extensive 100ns MD simulation was carried out on the complex structure of 6LU7 receptor with the selected approved drug (indinavir) and designed analog (compound 16) using Desmond software (D. E. Shaw Research, New York, NY, 2015) to access the binding stability of Ligand-6LU7 complex. 35 The system was solvated in TIP3P water model and 0.15 M NaCl to mimic a physiological ionic concentration. The full system energy minimization step was done for 100ps. The MD simulation was run for 100ns at 300K temperature, standard pressure (1.01325 bar), within an orthorhombic box with buffer dimensions 10× 10× 10 Å3 and NPT ensemble. The energy (kcal/mol) was recorded at an interval of 1.2 ps. The protein-ligand complex system was neutralized by balancing the net charge of the system by adding Na + or Clcounter ions. The Nose-Hoover chain and Martyna-Tobias-Klein dynamic algorithm was used to maintain the temperature of all the systems at 300 K and pressure 1.01325 barr, respectively.</p><!><p>The current COVID-19 pandemics urges the immediate therapeutic treatments to reduce the morbidity, and to protect the health of people who are at high risk of infection, particularly when vaccines are still under development. Currently, no specific treatment is available and the only treatment remains is supportive care i.e. oxygen therapy, fluid management, use of approved antiviral developed for Ebola (remdesivir) 36 , anti-HIV (lopinavir-ritonavir) 37 and recently apporoved antimalaria drugs (hydroxychloroquine) 38 as a combination therapy. The combination of azithromycin and hydroxychloroquine was also found to be effective as a treatment against COVID-19. 39 In order to develop targeted therapeutics of unknown infection from new virus, intitially, we need to design library of compounds based on the available target protein structure and active binding moiety of the protein residues. Although, it is hard to predict the potency and binding specificity from the structure of a chemical compound. The structure-based design of new bioactive molecules involves, as a crucial step, the development of a reaction strategy affording enough chemical diversity and structural variety of potentially bioactive products. As such, not just a single molecule, but a significant part of the chemical space becomes available for the biological investigation that facilitates SAR analysis. Of particular note, finding a potent molecule possessing the drug-like properties remains challenging, and therefore advanced drug discovery involves CADD approaches that plays a crucial role in the early discovery phases. The main protease, 3CL pro of SARS-CoV-2 was presented as an attractive drug target and a rapid identification of drug candidates against COVID-19 was achieved through virtually screening. 33,34 Encouraged with this and taking the advantage of a recently released crystal structure of 3CL pro from SARS-CoV-2 (PDB code 6LU7) 16 , we performed virtual screening for the 13 antiviral drugs. The list of the compounds are presented in Table 1.</p><!><p>In the begening, molecular docking of the approved antiviral drugs was carried out with protein 3CL pro using the Glide module of Schrödinger suite. 32,40 Docking score, and glide XPGscore (kcal/mol) were considered as the measures of binding affinity to rank the poses of the ligands.</p><p>Among the listed antiviral drugs, indinavir, which is an approved HIV-1 protease inhibitor, was realized as the best inhibitor based on docking score and XP Gscore of -8.824 and -9.466 kcal/mol, respectively (Table 1, Entry 1). Second rank was secured by another HIV-1 protease inhibitor atazanavir with a docking score of -7.912 (Table 1, Entry 2). Interactions of indinavir with 3CL pro protein are presented in the Fig. 3. Indinavir showed H-bond interactions with a residue Gln_189 in domain-I, and Glu_166 in domain-II of 3CL pro protein. It also formed saltbridge interaction with residue Glu_166 in domain-II (Table 4, Entry 1). Therefore, indinavir possessing an important pharmacophore, HEA was selected as a positive control to build a new library of HEA analogs that could effectively inhibit the main protease 3CL pro of SARS-CoV-2.</p><p>Our research studies [41][42][43] , have demonstrated HEA as an important pharmacophore. 44,45 Together, presence of HEA in indinavir 46,47 and our ongoing research towards discovery of HEA analogs, a library of new HEA compounds (>2500) was designed and subjected to virtual screening against 3CL pro protein of SARS-CoV-2. Of the screened >2500 new HEA analogs, 25 displayed a notable docking score over the indinavir as shown in Table 2. Among top-ranked analogs, compound 1 (Table 2, Entry 1) was shown to possess highest docking score (-9.864) and XP Gscore (-10.227 kcal/mol).</p><p>Next, to get an insight about the drug-likness properties, all the 25 top-ranked analogs, were screened for their ADME profiles and the results are depicted in Table 3. Notably, only compound 16 (Table 3, Entry 16) strictly abide to Lipkinski's "rule of five" and theortically met the druglikeness properties. Compound 4 (Table 3, Entry 4) also followed the criteria with a slight deviation in number of rotatable bond (i.e 11). 48 However, compound 4 displayed hepatoxicity. Based on the favorable physicochemical properties displayed by compound 16 was sleceted for further studies.</p><p>Compound 16 also offered an extra H-bond interaction, and also enhanced hydrophobically packed H-bond 32 (Phob EnHB) interaction by -1.5. As depicted in Fig. 4, the phenyl group of HEA moiety of compound 16 interacted with Leu_141 and Phe_140 residue by hydrophobic interaction (Table 4, Entry 2). The 4-fluoro-aniline scaffold was found to be packed deeply inside the pocket by hydrophobic interactions with Cys-44 and Cys-145 and pi-pi-interaction with His-41 in domain-I. Residues Cys_145, Ser_144, and Gly_143 interacted with carboxamide oxygen by H-bond, and His_164 in domain-I with nitrogen of aniline moiety as shown in Fig. 4.</p><p>SAR analysis was perfomed to designate the role of scaffolds present in HEA compounds. The lead compound, 16 contains three main functionalities i.e. HEA, anilines derivatives on pocket 1, and aromatic/aliphatic groups on pocket 2 as shown in Table 5. Firstly, the role of HEA backbone was studied, one compound 33 (Table 5, entry 2 and 5), and was considered for MD simulation study.</p><!><p>MD simulation at 100ns was performed for indinavir-3CL pro and compound 16-3CL pro . Unligated-3CL pro complex was also considered for MD studies to identify the effect on simulations in the absence of ligands. There are several factors such as conformation of ligand, water molecules, ions, cofactors, ligand protonation state, conformational and solvation entropies that can affect in silico predictions in an unexpected pattern. Several reports are available to support the role of MD simulations for the refinement of docking results. 49,50 MD simulations were carried out to determine the stability of the interactions of ligand-protein docked complexes. The final structure of simulated complexes exhibited strong stereochemical geometry of the residues as analysed by the Ramachandran map (Fig. 5).</p><p>The number and percentage of residues in favored, allowed, and outlier regions for these three systems are presented in Table 6. Both, complex of indinavir, and un-ligated 3CL pro possess two (Gly_251, Thr_224) and three (Gly_11, Asn_84, Asp_187) outlier residues, respectively (Table 6, entry 2-3). Interestingly, there was no outlier residues for 16-3CL pro complex (Table 6, entry 1). Further, the stability of compound 16-3CL pro complex was monitored during the entire simulation process where total energy (E), potential energy (E_p), temperature (T), pressure (P), and volume (V) values were computed as shown in Fig. 6. The plots for un-ligated and indinavir complexes are displayed in supporting information (Fig. S1). Notably, no significant change was observed in potential energy and total energy and other parameters for compound 16 and indinavir complexes, as the average values with a standard deviation of these parameters were found in same range (supporting information, Table S1-S2). While entire simulations, Root Mean Square Deviation (RMSD) change for Cα was also monitored. This parameter measures the global deviation of atoms from a reference status i.e., frame 0. The RMSD plot indicated the fluctuations in the initial conformation of the receptor for all three systems till 30 ns (Fig. 7), which later stabilized in production phase with an average values of RMSDbackbone (2.015 Å), RMSDCα (1.993 Å) and RMSDside-chain (2.876 Å), for 16-3CL pro complex. The average RMSD values for backbone, Cα, and side-chain for indinavir-3CL pro and un-ligated-3CL pro complexes are depicted in supporting information (Table S3). The RMSD of the both Cα and backbone for compound 16-3CL pro showed fluctuations between the range of 0.766 -4.304 Å. Protein-RMSF was monitored to assess the local residue flexibility, and ligand-RMSF was examined to study the the atom-wise fluctuations in the ligand. When a dynamical system fluctuates over well-defined average position, the RMSD from the average over time can be referred to as RMSF. The RMSF, 'Fit on Protein' trend (brown line; Fig. 8B) indicated the ligand fluctuations with respect to the binding site residues present on target protein. Whereas 'Fit on Ligand' trend (pink line, Fig. 8B) showed the fluctuations where the ligand in each frame is aligned on the ligand in the first reference frame. Both, compound 16 and indinavir interactions with 3CLpro are shown in Fig. 8. Compound 16 interacted with residues (Gly_143, Ser_144, Cys_145, and His_164) upon docking in absence of water. In the presence of water molecules, only one hydrogen bond interaction was persistent with Gly_143, observed during MD simulation. This observation showed the importance of water molecules within binding pocket of protein 3CL pro for ligand 16. The detailed interaction of compound 16 with binding site residues of 3CL pro were studied and are presented in Fig. 9, where interacting residues include Thr_26, Leu_27, His_41, Ser_46, Met_49, Asn_119, Phe_140, Leu_141, Asn_142, Gly_143, Ser_144, Cys_145, His_163, His_164, Met_165, Glu_166, Pro_168, His_172, and Gln_189. The interacting residue with compound 16 are presented with green color, while the orange and blue bands indicate protein secondary structures i.e. helices and β-strands, respectively (Fig. 9). These binding site residues possessed RMSF value of <2Å.</p><p>The distribution of residues index of protein in secondary structure elements (SSE) was also studied (Fig. S3; supporting information). There are two panels in SSE i.e. top panel and bottom panel, where the top panel reports the percent SSE composition from helices and strands over time, and the bottom panel monitors each residue and its SSE assignment over the time. SSE percentage for compound 16, indinavir, and un-ligated protein with values, 44.45%, 44.74%, and 44.45%, respectively, indicated negligible while 100ns simulation period (Table S4; supporting information). Compound 16 interacted with Leu_27, His_41, Met_49, and Met_165 residues of protein forming hydrophobic interactions, and with Asn_142, Gly_143, Ser_144, Cys_145, His_163, His_164, and Glu_166 formed H-bond-water interactions. There was no detectable ionic interaction for 16-3CL pro complex as depicted in Fig. 10. The total number of contacts formed by protein with the ligand 16 over the course of the trajectory was presented in the top panel and the bottom panel (Fig. 11) that showed the residues interacted with the ligand in each trajectory frame. It was observed that some protein residues formed more than one specific contact with the ligand (presented with a darker shade of orange).</p><p>Overall, six properties were analysed to explain the stability of the compound 16 in the 3CL pro receptor during the simulation of 100ns: a) ligand RMSD -Root mean square deviation of a ligand with respect to the reference conformation (typically the first frame was used as the reference and it is regarded as time t=0); b) radius of Gyration (rGyr) -it is used for measuring the 'extendedness' of a ligand, and is equivalent to its principal moment of inertia; c) NS34 is not capable of intramolecular hydrogen bonding, so this plot appears empty; d) Molecular Surface Area (MolSA) -molecular surface calculation with 1.4 Å probe radius; e) Solvent Accessible Surface Area (SASA) -surface area of a molecule accessible by a water molecule; f) Polar Surface Area (PSA) -Solvent accessible surface area in a molecule contributed only by oxygen and nitrogen atoms as shown in Fig. 12.</p><p>From the graph, it is evident that ligand RMSD remains constant during the simulation process.</p><p>The overall RMSD of the compound 16 was up to 2.771 Å. In the initial stage, fluctuation was observed from 0ns to 18ns afterwards constant RMSD was observed during the entire simulation process. Fluctuation for Radius of Gyration were noted until 40 ns and subsequently, a stable conformation was acquired throughout the simulation. The radius of gyration throughout the 100ns simulation was ranged from 3.275 Å to 4.524 Å for compound 16. The SASA plot revealed a fluctuating pattern for initial 42ns and then became stabilized until simulations were completed. MolSA, and PSA plots also indicated the consistency of the ligand, 16 during the simulation process. Initially, compound 16 showed intramolecular Hbond, which became stabilised after protein interaction and no intramolecular H-bond was observed at the end.</p><p>A 2D schematic representation for compound 16 is presented with color coded rotatable bonds (Fig. 13). The rotatable torsional bond of compound 16 was accompanied by a dial (radial) plot and bar plots of the same color. The radial plot explains the conformation of the torsion angle throughout the time course of the simulation. The simulation starts from the centre of the radial plot and the time progress was plotted radially towards the outside. The probability density of the torsion angle was summarized by bar plots, which were recapitulating the data on the radial plots. The Y-axis of the bar plots revealed the potential of the rotatable bond, which are expressed in kcal/mol. The radial and bar diagram explained the torsion potential relationships and the conformational strain of compound 16 undergoing to maintain a protein-bound conformation.</p><!><p>COVID-19 is continue to become a global burden on human health and economic losses remains unabated. A rapid progress, particularly in vaccine and therapeutic development, is essential to overcome the further explosion in spread and loss of human life from COVID-19.</p><p>The availability of a protein target structure is usually helpful in identifying and designing potential drug candidates. Such approaches usually involve explicit molecular docking of ligands into the receptor-binding site, producing a predicted binding mode for each candidate compound. Virtual screening is widely used for structure-based drug design against a known drug target to identify those chemical structures, which are most likely to bind well within the active site of the target protein. Newly released 3CL pro structure of SARS-CoV-2 was considered as potential target for the design of drug candidates and a rapid virtual screening for the libraries of approved drugs. Such approaches are certainly important to get better insight about the targeted therapeutic development for this novel infectious disease. In this article, we have screened approved antiviral drugs against 3CL pro , and indinavir was identified with highest docking score and XP Gscore and thus selected as a positive control for the design of new drug candidates based on HEA scaffold. Molecular docking of new compounds was carried out and a total of 25 analogs were identified with improved or similar docking score as of indinavir. Both, compound 16 and indinavir bind to the domain I, II, and domain II-III linker of targeted protein with H-bond, pi-pi and hydrophobic interactions. Compound 16 and indinavir were further studied for MD simulation. RMSF for indinavir and 16 complex were below 2 Å. Interestingly, no outlier residue was noted for compound 16-3CL pro complex as compared to two outlier residues found in indinavir-3CL pro complex. Morever, presence of water molecules within binding site of protein suggested stability of 16-3CL pro complex.</p><p>RMSD and Ligand-RMSF percentage for Cα demonstrated that stability of 16-3CL pro complex and a protein-bound conformation was confirmed by torsional analysis, suggesting suitability as a potent candidate for further exploration as a candidacy as an explorative drug for use against SARS-CoV-2.</p><!><p>Brijesh Rathi, PhD brijeshrathi@hrc.du.ac.in; brathi@luc.edu Caption for Figures.     Ramachandran plot for 3CL pro in case of: A) Compound 16-3CL pro complex; B) Indinavir-3CL pro complex; C) un-ligated-3CL pro complex.</p>