Discovery of New Hydroxyethylamine Analogs against 3CLpro Protein Target of SARS-CoV-2: Molecular Docking, Molecular Dynamics Simulation, and Structure-Activity Relationship Studies.

The 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 a vaccine and potential therapeutics are critically essential. The crystal structure for the main protease (Mpro) of SARS-CoV-2, 3-chymotrypsin-like cysteine protease (3CLpro), was recently made available and is considerably similar to the previously reported SARS-CoV. Due to its essentiality in viral replication, it represents a potential drug target. Herein, a computer-aided drug design (CADD) approach was implemented for the initial screening of 13 approved antiviral drugs. Molecular docking of 13 antivirals against the 3-chymotrypsin-like cysteine protease (3CLpro) 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 a 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 analysis performed at 100 ns supported the stability of 16 within the binding pocket. Largely, our results supported that this novel compound 16 binds with domains I and II, and the domain II-III linker of the 3CLpro protein, suggesting its suitability as a strong candidate for therapeutic discovery against COVID-19.

Introduction

Coronaviruses (CoV) belong to a group of viruses consisting of a core of genetic material enveloped with a protein spike appearing like a crown, which means “corona” in Latin.[1] A 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 are 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 an animal market in China.[7] Subsequently, the outbreak was attributed to a novel CoV and related to the SARS virus based on its genetic similarities with a previous known coronavirus (SARS-CoV).[8] Later, the disease was named as COVID-19 by the World Health Organization (WHO) caused by a novel coronavirus, SARS-CoV-2.[9]

The outbreak started in Wuhan, China, and rapidly escalated to other countries. The USA, Spain, Russia, UK, Italy, France, and Brazil lead with the highest number of COVID-19 cases officially reported. COVID-19 caused approximately 2924 deaths, and 85,403 confirmed cases were identified with a mortality ratio of 3.42% until 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 10 May 2020, a total of 38,55,788 confirmed cases have been noted all around the world, resulting in 2,65,862 deaths.[12] The number of confirmed cases has increased sharply, and also, the number of deaths due to COVID-19 can be clearly observed across the world during the month of April 2020 as shown in Figure .[13]

Figure 1

Number of COVID-19 cases between 1 and 30 April 2020 globally: (A) confirmed cases and B) deaths reported (Source: ref (13)).

Major symptoms of this disease include high fever, cough, and shortness of breath, whereas in severe cases pneumonia and kidney failure are the major cause for death.[14] However, in many cases, the basic symptoms of this deadly disease were not observed, and asymptomatic transmission of the virus from the infected person could be more dangerous. Although, a series of rapid and specific diagnostic tools are available for the COVID-19 disease,[15] vaccines and SARS-CoV-2 specific therapeutic treatments are not available. Recently, the crystal structure of the SARS-CoV-2 main protease, Mpro, also called 3CLpro, complexed with an inhibitor N3 was released.[16] The 3CLpro enzyme of SARS-CoV-2 processes polyproteins by proteolytic action of replicase enzyme (pp1a and pp1ab) to release a functional polypeptide (Figure A).[16] It is a dimeric protein, which contains two symmetric units designated as protomers. Each protomer has three domains, namely, domain I (residues 8–101), domain II (residues 102–184), and domain III (residues 201–303). Domain III comprises five α-helices and is linked with domain II through an extended loop region (residues 185–200). 3CLpro has Cys_145 and His_41 catalytic dyads, and a substrate-binding site is positioned in the cleft between domains I and II (Figure B). These descriptions match with the previously reported protease enzyme of SARS-CoV.[17−23]

Figure 2

3CLpro enzyme of SARS-CoV-2 displaying (A) proteolytic action of replicase enzyme (pp1a and pp1ab) to release nsp16 and (B) three domains, namely, domain I, domain II, and domain III (residues Cys_145 and His_41 are part of catalytic dyads).

3CLpro of SARS-CoV-2 is an important drug target, and computer-aided drug design (CADD) is considered as an undisputable and significant approach to discover antiviral drug candidates. Identically, studies have been adopting various approaches for therapeutic discovery against SARS-CoV-2 through virtual screening.[33,34] Antiviral drugs such as HIV protease inhibitors[24] and approved drugs such as ribavirin, oseltamivir, and a few more have been studied against earlier reports for SARS-CoV.[25] In this article, we aim to perform a fast discovery of the potential candidates against 3CLpro of SARS-CoV-2 through virtual screening of approved antiviral therapeutics and build a focused library of novel potential compounds.

Materials and Methods

The protease structure, the 3-chymotrypsin-like cysteine protease (3CLpro) enzyme of SARS-CoV-2 (PDB ID: 6LU7) with 2.1 Å, was obtained from the the Protein Data Bank at the RCSB site (http://www.rcsb.org). The computational work was performed using Schrödinger software.

Preparation of Protease Structure

The SARS-CoV-2 virus protein structure was prepared in the Protein Preparation Wizard and Prime module of the Schrödinger suite to remove defects such as missing hydrogen atoms, inappropriate bond order assignments, charge states, alignments of several groups, and missing side chains.[26,27] Steric hindrance and strained bonds/angles were removed through restrained energy minimization, permitting movement in heavy atoms up to 0.3 Å.

Preparation of Antiviral Drugs and HEA Analogs for Virtual Screening

Approved antiviral drugs such as indinavir, saquinavir, nelfinavir, ritonavir, lopinavir, atazanavir, amprenavir, darunavir, nelfinavir, oseltamivir, tipranavir, fosamprenavir, and galidesivir were screened to select the hits as reference molecule or positive control for the design of novel compounds. The coordinate files of approved drugs were obtained from the Protein Data Bank at the RCSB site (http://www.rcsb.org) and PubChem.[28,29] Over 1000 analogs were prepared by the Maestro tool. The structures of both approved antiviral drugs and HEA analogs were prepared prior to docking by Ligprep. Ligprep assisted the production of 2D or 3D structures and their corresponding low-energy 3D structures, which were immediately available for Glide program applications. All the parameters were kept as the default except chirality parameters for the antiviral drugs and new HEA analogs. The chirality was kept as the default for the 3D structures, and all the combinations of chirality were developed for the approved antiviral drugs and new HEA analogs. Next, desalting, tautomer generation, and probable ionization states at pH 7 ± 2 were adjusted.[30] Ionization states of all the molecules were predicted by the Schrödinger suite inbuilt Epik module.[31]

Molecular Docking of Compound Libraries against Target Protease Active Site

Site-specific molecular docking of approved antiviral drugs (Table ) and novel HEA-based analogs (Table ) against protease was accomplished using the Glide module of the Schrödinger suite. The Glide tool was used for receptor grid preparation at default parameters, which include a van der Waals radius scaling factor (1.0) and partial charge cutoff (0.25).[32] The screening of both libraries was performed by Glide at extra precision (XP), which indicates an appropriate connection between excellent poses and good scores.

Table 1

Docking Score and XP Gscore (kcal/mol) for All 13 Antiviral Drugs Considered for the Study

Entry	Drug	Docking score	XP Gscore (kcal/mol)
1	Indinavir	–8.824	–9.466
2	Atazanavir	–7.912	–7.92
3	Remdesvir	–7.804	–7.804
4	Amprenavir	–7.747	–7.747
5	Saquinavir	–7.455	–7.468
6	Ritonavir	–7.422	–7.422
7	Lopinavir	–7.041	–7.041
8	Darunavir	–7.028	–7.028
9	Nelfinavir	–6.73	–6.744
10	Oseltamivir	–5.825	–5.907
11	Tipranavir	–5.64	–5.778
12	Fosamprenavir	–5.309	–6.578
13	Galidesvir	–4.967	–6.322

Table 2

Docking Score and XP Gscore (kcal/mol) for Identified HEA-Based Hit Analogs (1–25)

ADMET Calculations

Swiss ADME[33] and PKCSM[34] were employed to predict the ADMET (i.e., absorption, distribution, metabolism, excretion, and toxicity) profile of identified hit compounds (1–25) and two antiviral drugs (indinavir and atazanavir). 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 (drug-likeness), Rat LD50, and hepatotoxicity. The results are shown in Table .

Table 3

Physicochemical Properties of HEA-Based Analogs (1–25) and Anitiviral Drugs (Indinavir and Atazanavir)a

Entry	Compound	MW	nON	nOHNH	Nrot	TPSA	mLogP	LD50	Hept.	nVio.
1	1	717.9	8	5	18	177.2	0.47	2.88	Yes	2
2	2	935.1	12	4	22	179.9	2.83	2.48	Yes	2
3	3	821.0	6	4	16	105.1	4.74	2.48	Yes	1
4	4	442.5	5	3	11	61.4	4.68	2.32	Yes	1
5	5	1067	18	4	26	198.3	1.96	2.47	Yes	2
6	6	995.1	10	4	22	179.9	3.23	2.48	Yes	2
7	7	1035	16	4	24	179.9	3.17	2.45	Yes	2
8	8	1011	10	4	24	179.9	3.61	2.46	Yes	2
9	9	899.1	10	4	22	179.9	2.41	2.45	Yes	2
10	10	1011	10	4	24	179.9	3.61	2.47	Yes	2
11	11	1003	12	4	22	179.9	3.35	2.48	Yes	2
12	12	715.8	9	2	14	123.9	3.41	2.52	Yes	2
13	13	600.7	5	3	12	98.9	3.18	2.44	Yes	1
14	14	1103	16	4	24	179.9	3.66	2.48	Yes	2
15	15	899.1	10	4	22	179.9	2.41	2.45	Yes	2
16	16	374.4	4	3	10	70.6	3.31	2.79	No	0
17	17	871.0	10	4	20	179.9	2.1	2.44	Yes	2
18	18	1035	16	4	24	179.9	3.17	2.46	Yes	2
19	19	935.1	12	4	22	179.9	2.83	2.46	Yes	2
20	20	683.8	8	2	14	123.9	2.89	2.51	Yes	2
21	21	1135	18	4	16	198.3	2.46	2.48	Yes	2
22	22	927.1	10	4	22	179.9	2.72	2.45	Yes	2
23	23	524.6	5	3	10	98.9	2.3	2.41	Yes	1
24	24	652.8	6	3	15	111.2	3.14	2.35	Yes	1
25	25	1079	10	4	24	179.9	4.08	2.48	Yes	2
26	Indinavir	613.8	7	4	14	118	1.33	2.91	Yes	1
27	Atazanavir	704.9	9	5	22	171.2	1.76	2.67	Yes	2

MW = molecular weight (g/mol); nON = no. of hydrogen bond acceptor; nOHNH = no. of hydrogen bond donors; Nrot = no. of rotatable bonds; TPSA = total polar surface area; MLogP = predicted octanol/water partition coefficient; LD50 = oral rat acute toxicity; Hept. = hepatoxicity; nVio. = no. of Lipinski violation.

Molecular Dynamics (MD) Simulation

An extensive 100 ns MD simulation was performed for the complex structure of the 6LU7 receptor with the selected approved drug (indinavir) and designed novel HEA compound (compound 16) using Desmond software (D. E. Shaw Research, New York) to access the stability of binding for the ligand-6LU7 complex.[35] The system was solvated in a TIP3P water model and 0.15 M NaCl to mimic a physiological ionic concentration. The full system energy minimization step was done for 100 ps. The MD simulation was run for 100 ns at a 300 K temperature and standard pressure (1.01325 bar), within an orthorhombic box with buffer dimensions of 10 Å × 10 Å × 10 Å and NPT ensemble. The energy (kcal/mol) was recorded at an interval of 1.2 ps. Neutralization of the protein-ligand complex system was achieved by adding Na+ or Cl– counterions that balanced the net charge of the systems. The Nosé–Hoover chain and Martyna–Tobias–Klein dynamic algorithm were used to sustain the temperature of all the systems at 300 K and pressure of 1.01325 bar, respectively.

Result and Discussion

The current COVID-19 pandemic urges the immediate therapeutic treatments to protect the health of people who are at high risk of infection, particularly when vaccines are still under development. To date, no specific therapeutic treatment is available, and current treatment includes supportive care, i.e., oxygen therapy, fluid management, use of approved antiviral drugs developed for Ebola (remdesivir),[36] anti-HIV (lopinavir-ritonavir),[37] and recently approved antimalaria drugs (hydroxychloroquine)[38] as a combination therapy. The combination of azithromycin and hydroxychloroquine was also shown as a partial treatment against COVID-19.[39] In order to develop targeted therapeutics, initially, design of compound libraries based on the available target protein structures and active binding moiety of the protein residues is an important step. Although it is hard to predict the potency and binding specificity from the structure of a chemical compound, structure-based design of novel bioactive scaffolds is considered as an important approach to accomplish adequate chemical diversity of compound libraries that facilitate SAR analysis. Of particular note, finding a potent molecule with drug-like properties remains challenging, and therefore, advanced drug discovery involves CADD approaches that accelerates the early discovery phases. The main protease, 3CLpro, of SARS-CoV-2 was presented as an attractive drug target, and a rapid identification of drug candidates was achieved through virtual screening.[33,34] Encouraged with this, and exploiting a recently published crystal structure of 3CLpro from SARS-CoV-2 (PDB code 6LU7),[16] we executed virtual screening for the 13 antiviral drugs. The list of the compounds is presented in Table .

Virtual Screening of Antiviral Drugs and New HEA Analogs against SARS-CoV-2 Protein 3CLpro

In the beginning, molecular docking of the approved antiviral drugs was carried out with protein 3CLpro using the Glide module of the Schrödinger suite.[32,40] Docking score and Glide XP Gscore (kcal/mol) were considered as the measures for binding affinity to rank the poses of the ligands. Among the listed antiviral drugs, indinavir, which is an approved HIV-1 protease inhibitor, was realized as the best inhibitor based on a docking score and XP Gscore of −8.824 and −9.466 kcal/mol, respectively (Table , entry 1). The second rank was secured by another HIV-1 protease inhibitor, atazanavir, with a docking score of −7.912 and XP Gscore of −7.92 (Table , entry 2). Our results were in agreement with the docking score of indinavir (−10.424) reported by Chang[41] and colleagues on AutoDock Vina. Interactions of indinavir with the 3CLpro protein are presented in Figure . Indinavir showed H-bond interactions with a residue Gln_189 in domain I and Glu_166 in domain II of the 3CLpro protein. It also formed a salt-bridge interaction with residue Glu_166 in domain II (Table , entry 1). Therefore, indinavir, possessing an important HEA pharmacophore, was selected as a positive control to build a new library of HEA analogs that could effectively intervene with the functioning of the main protease 3CLpro of SARS-CoV-2.

Figure 3

Indinavir-3CLpro docked complex, showing interaction with crystal structure 6LU7 (brown color) and surrounded binding site residue.

Table 4

List of Interacting Residues of 3CLpro with Indinavir and Lead Compound 16

Our research studies[42−44] have demonstrated HEA as an important pharmacophore.[45,46] With the presence of HEA in indinavir[47,48] and our ongoing research toward the discovery of HEA analogs, a library of new HEA compounds (>2500) was designed and subjected to virtual screening against the 3CLpro protein. Of the screened >2500 new HEA analogs, 25 analogs exhibited a notable docking score over the indinavir score as shown in Table and Table S1.

Among the top-ranked analogs, compound 1 (Table , entry 1) was shown to possess the highest docking score (−9.864) and XP Gscore (−10.227 kcal/mol).

Next, to get an insight about the drug-likeness properties, all the 25 top-ranked analogs along with two antivirals (indinavir and atazanavir) were screened for their ADMET profiles, and the results are depicted in Table . Notably, only compound 16 (Table , entry 16) strictly abides to the Lipinski’s “rule of five” and theoretically met the drug-likeness properties. Compound 4 (Table , entry 4) also followed the criteria with a slight deviation in number of rotatable bonds (i.e., 11).[49] However, compound 4 exhibited hepatoxicity. On the basis of the favorable physicochemical properties displayed, compound 16 was selected for further studies.

Compound 16 also offered an extra H-bond interaction and also enhanced hydrophobically packed H-bond[32] (Phob EnHB) interaction by −1.5. The phenyl group of a HEA moiety of compound 16 interacted with Leu_141 and Phe_140 residues by hydrophobic interactions (Table , 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-bonds and His_164 in domain I with nitrogen of an aniline moiety as shown in Figure .

Figure 4

Compound 16-3CLpro docked complex, showing interaction with crystal structure 6LU7 (brown color) and surrounded binding site residue.

Structure activity relationship (SAR) analysis was performed 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 . First, the role of the HEA backbone was studied; one compound 33 (Table , entry 8) that does not contain HEA and lacking the stereogenic center secured the lowest docking score (i.e., −7.121) highlighting the importance of HEA. Second, variations made on pocket 1 were analyzed keeping the tert-butylacetate group at pocket 2 constant. Compound 16 possessing a docking score of −8.955 contains a 4-fluoroaniline at pocket 1 and found to be the best analog among others (i.e., 29, −8.131; 30, −7.995; and 32, −7.949) that contains 2-fluoroaniline, 2,4-difluoroaniline, and 3-fluoroaniline, respectively (Table , entries 4, 5, and 7). The results clearly indicated the significance of the fluorine group at the para position of aniline as compared to ortho and meta positions. Similarly, chemical diversity at pocket 2 was correlated, while 4-fluoroanline at pocket 1 remains constant. Compound 26 possessing a docking score of −8.835 contains a 4-trifluoromethylbenzyl moiety at pocket 2 and was found to display an improved docking score compared to 27 (−8.55), 28 (−8.543), and 31 (−7.965) that contain similar moieties 3-fluoro-4-trifluorophenyl, 4-fluorophenyl, and 4-trifluorophenyl moiety, respectively (Table , enties 1–3 and 6). However, 26 was shown to possess an inferior docking score in comparison to compound 16 that contains a tert-butylacetate group at pocket 2, indicating more pronounced effects without fluorinated aromatic components. Therefore, on the basis of molecular docking, ADMET properties, and SAR analysis, compound 16 was identified as the most suitable analog and was considered for MD simulation study.

Table 5

SAR Analysis

Molecular Dynamics (MD) Simulation Studies

MD simulation at 100 ns was performed for indinavir-3CLpro and compound 16-3CLpro. The unligated-3CLpro complex was also considered for MD studies to identify the effect on simulations in the absence of ligands. There are several factors such as the conformation of ligands, water molecules, ions, cofactors, ligand protonation state, and 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.[50,51] MD simulations were carried out to determine the stability of the interactions of ligand-protein docked complexes. The final structures of simulated complexes exhibited strong stereochemical geometries of the residues as analyzed by the Ramachandran map (Figure ).

Figure 5

Ramachandran plot depicting stereochemical geometry for (A) compound 16-3CLpro complex, (B) indinavir-3CLpro complex, and (C) unligated-3CLpro complex.

The number and percentage of residues in favored, allowed, and outlier regions for these three systems are presented in Table . The complexes of indinavir and unligated 3CLpro possess two (Gly_251, Thr_224) and three (Gly_11, Asn_84, Asp_187) outlier residues, respectively. Interestingly, there was no outlier residues for the 16-3CLpro complex (Table , entry 1). Further, the stability of the compound 16-3CLpro 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 Figure . The plots for the indinavir-3CLpro complex is displayed in the Supporting Information (Figure 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, Tables S2 and S3). While for 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 until 30 ns (Figure ), which later stabilized in the production phase with average values of RMSDbackbone (2.015 Å), RMSDCα (1.993 Å), and RMSDside-chain (2.876 Å) for the 16-3CLpro complex. The average RMSD values for backbone, Cα, and side chain for indinavir-3CLpro and unligated-3CLpro complexes are depicted in the Supporting Information (Table S4). The RMSD of both the Cα and backbone for compound 16-3CLpro showed fluctuations between the range of 0.766–4.304 Å. The ligand RMSD values for compound 16 and indinavir with 3CLpro are shown in Figure A. Protein-RMSF was monitored to assess the local residue flexibility, and ligand-RMSF was examined to study the atom-wise fluctuations in the ligand. With fluctuation of a dynamical system over the well-defined average position, the RMSD from the average over time may be referred as RMSF.[35] The RMSF values, “Fit on Protein” trend (brown line) and “Fit on Ligand” trend (pink line), are depicted in Figure B for compound 16 and in Figure S2 for indinavir. Compound 16 interacted with residues (Gly_143, Ser_144, Cys_145, and His_164) upon docking in the absence of water. In the presence of water molecules, only one hydrogen bond interaction was persistent with Gly_143, observed during MD simulation (Table , entry 2). This observation showed the importance of water molecules within binding pockets of protein 3CLpro for ligand 16. The detailed interaction of compound 16 with binding site residues of 3CLpro were studied and are presented in Figure , 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 is presented with a green color, while the orange and blue bands indicate protein secondary structures, i.e., helices and β-strands, respectively (Figure ). These binding site residues possessed RMSF values of <2 Å.

Table 6

Number and Percentage of Residues Existing in Favored, Allowed, and Outlier Region for Simulated System

Entry	System	Number and percentage of residues in favored region	Number and percentage of residues in allowed region	Number and percentage of residues in the outlier
1	Compound 16-3CLpro	274 (90.1%)	30 (9.9%)	0.0
2	Indinavir-3CLpro	282 (92.8%	20 (6.6%)	2 (0.7%)
3	Unligated 3CLpro	280 (92.1%)	21 (6.9%)	3 (1.0%)

Figure 6

Simulation quality parameters for compound 16-3CLpro complex for 100 ns simulation in terms of total energy (E), potential energy (E_p), pressure (P), temperature (T), and volume (V).

Figure 7

RMSD plot for Cα of 3CLpro: (a) 16-3CLpro complex (blue color), (b) indinavir-3CLpro complex (orange color), and (c) unligated-3CLpro complex (gray color).

Figure 8

(A) Ligand RMSD plot for compound 16 and indinavir and (B) ligand-RMSF plot for compound 16-3CLpro complex, where brown line indicates the ligand fluctuations with respect to the binding site residues present on target protein and pink line shows the fluctuations where the ligand in each frame is aligned on the ligand in the first reference frame.

Figure 9

RMSF plot for Cα of 3CLpro residues in compound 16-3CLpro complex; residues are shown in three-letter code with green color belonging to binding site interacting with compound 16.

The distribution of the residue index of protein in secondary structure elements (SSE) was also studied (Figures S3–S5, Supporting Information). The SSE percentage for compound 16, indinavir, and unligated protein with values of 44.45%, 44.74%, and 44.45%, respectively, indicated negligible amounts with a 100 ns simulation period (Table S5, Supporting Information). Compound 16 interacted with Leu_27, His_41, Met_49, and Met_165 residues of proteins forming hydrophobic interactions and with Asn_142, Gly_143, Ser_144, Cys_145, His_163, His_164, and Glu_166 forming H-bond-water interactions. There was no detectable ionic interaction for the 16-3CLpro complex as depicted in Figure . The contacts formed by proteins with ligand 16 over the course of the trajectory are presented in the top panel and bottom panel in Figure , showing that the residues interacted with the ligand in each trajectory frame. Some protein residues showed more than one specific contact with the ligand, which are described with a darker shade of orange color. Overall, six parameters[35] were analyzed to elucidate the stability of ligand 16 in the 3CLpro receptor during the simulation of 100 ns as described in Figure . While the simulation process ligand RMSD remained constant, the overall RMSD of compound 16 was up to 2.771 Å. In the initial stage, fluctuation was observed from 0 to 18 ns, and afterward, a constant RMSD was observed during the entire simulation process. Fluctuation for the radius of gyration were noted until 40 ns, and subsequently, a stable conformation was obtained during the complete simulation process. The radius of gyration throughout the 100 ns simulation ranged from 3.275 to 4.524 Å for compound 16. The SASA plot revealed a fluctuating pattern for the initial 42 ns and then became stabilized until simulations were completed. MolSA and PSA plots also suggested stability of ligand 16 during the simulation process. Initially, compound 16 showed intramolecular H-bonds, which became stabilized after protein interaction, and no intramolecular H-bonds were observed at the end.

Figure 10

Histogram (stacked bar chart) showing compound 16-3CLpro forming H-bonds interactions (green color), hydrophobic interactions (gray violet color), and water bridges (blue color) during 100 ns simulation.

Figure 11

Timeline representation of interactions and contacts (H-bonds, hydrophobic, water bridges) in compound 16-3CLpro complex during 100 ns MD simulation.

Figure 12

Ligand properties during 100 ns simulations for compound 16: (A) ligand RMSD, root mean square deviation, (B) radius of gyration (rGyr), (C) NS34, (D) molecular surface area (MolSA), (E) solvent accessible surface area (SASA), and (F) polar surface area (PSA).

A 2D schematic representation for compound 16 is presented with color-coded rotatable bonds (Figure ). The rotatable torsional bond of compound 16 was supplemented by a dial (radial) plot and the same color bar plots. The radial plot explains the torsion angle conformation during simulation. The simulation starts from the center of the radial plot, and the time progress was plotted radially toward 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 conserving a protein-bound conformation.

Figure 13

Torsional analysis of ligand-16 NS34 conformations during 100 ns simulations.

Conclusion

COVID-19 is continuing to become a global burden on human health, and economic losses remain 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. The structure of a protein target is considerably helpful in identifying and designing the potential drug candidates. Virtual screening remains one of the incredible approaches for structure-based design of chemical scaffolds against a known target protein binding well within the active site.

In summary, we have screened approved antiviral drugs against 3CLpro. Our in silico screening results indicated that indinavir, a known inhibitor of HIV protease, has high potential against the main protease, 3CLpro of SARS-CoV-2, and thus was selected as a positive control for the design of new drug candidates based on the HEA scaffold. Molecular docking of newly designed compounds was carried out, and a total of 25 analogs were identified with improved or similar docking scores with respect to indinavir. Further, compound 16 was identified as a top hit from the designed library based upon ligand–protein interaction, SAR analysis, and ADMET profile. Both compound 16 and indinavir bind to domains I and II and the domain II–III linker of the targeted protein with H-bond, pi-pi, and hydrophobic interactions. Compound 16 and indinavir were further studied for MD simulation. RMSF values for indinavir and the 16 complex were below 2 Å. Interestingly, no outlier residue was noted for the compound 16-3CLpro complex as compared to two outlier residues found in the indinavir-3CLpro complex. Moreover, the presence of water molecules within the binding site of the protein suggested stability of the 16-3CLpro complex. RMSD and the ligand-RMSF percentage for Cα demonstrated that the stabilities of the 16-3CLpro complex and a protein-bound conformation was confirmed by torsional analysis, suggesting suitability as a potent candidate for further in vitro and in vivo explorations as an explorative drug for use against SARS-CoV-2.