In search of SARS CoV-2 replication inhibitors: Virtual screening, molecular dynamics simulations and ADMET analysis.

Severe acute respiratory syndrome has relapsed recently as novel coronavirus causing a life threat to the entire world in the absence of an effective therapy. To hamper the replication of the deadly SARS CoV-2 inside the host cells, systematic in silico virtual screening of total 267,324 ligands from Asinex EliteSynergy and BioDesign libraries has been performed using AutoDock Vina against RdRp. The molecular modeling studies revealed the identification of twenty-one macrocyclic hits (2-22) with better binding energy than remdesivir (1), marketed SARS CoV-2 inhibitor. Further, the analysis using rules for drug-likeness and their ADMET profile revealed the candidature of these hits due to superior oral bioavailability and druggability. Further, the MD simulation studies of top two hits (2 and 3) performed using GROMACS 2020.1 for 10 ns revealed their stability into the docked complexes. These results provide an important breakthrough in the design of macrocyclic hits as SARS CoV-2 RNA replicase inhibitor.

Introduction

Novel coronavirus (CoV) has been recognized as a deadly respiratory pathogen and causative agent for severe acute respiratory syndrome coronavirus-2 (SARS CoV-2) or coronavirus disease 2019 (COVID-19) [1] since its outbreak in the Hubei province of Wuhan, China in 2019 [2,3]. It has caused an unprecedented pandemic worldwide with the higher fatality rate than the previously known SARS CoV. The background of the virus traces back nearly two decades to a closely related previously emerged pathogen having genomic similarity of 75–80% [4]. The epidemiological study on the virus brings in to notice the first case to be reported on December 2019 in China, which was then observed with the aggressive or violent form by January 2019 with increasing number of cases transmitting in the other parts of the world [5,6]. By January 2019, approximately eighteen countries were reported to be infected with the virus emerged in China led to the declaration as “Public Health Emergency of International Concern” by World Health Organization (WHO) on Mar 11, 2020 [7]. As on Jul 9, 2021 there have been reported 185,291,530 confirmed cases with 4,010,834 fatali cases worldwide [8]. The entire world has to suffer from the huge numbers of challenges due to this global pandemic as all the vital sectors like health, education, economic, social, pharmaceutical, etc. have been widely affected. Unemployment, events cancellation, prohibited public places, declined foreign and domestic trades during the pandemic caused a huge fall in the global economy [9,10].

The novel SARS CoV-2 viral pathogens responsible for the global pandemic have been derived from the previously emerged strain reported in 2002–03 which are suspected to be transmitted into humans through several animal hosts such as bats, civets, pangolins and camels [11]. The clinical symptoms in the patients infected with coronaviruses may vary from moderate symptomatic to severe asymptomatic characteristics including headache, cough, cold, tiredness, aches, fever, loss of taste or smell, diarrhoea, discoloration of finger or toes, breathing difficulties to dysfunction of organs. The viruses such as SARS and MERS are transmitted into humans through the direct or indirect exposure of the respiratory secretions [12,13].

The large genome of the coronaviruses including 30 kb long nucleotide sequence has been comprised within the single stranded positive sense enveloped RNAs [14]. The crown like structural morphology of these viruses observed under electron microscope are due to the club-like glycoproteins [13]. The CoV belongs to the β-coronaviruses among the four major genuses of the large family Nidovirales [15]. The critical insights into the viral pathogenesis and replication may shine a light towards the comprehensive understanding regarding the virus as presented in Fig. 1 . The pathogenesis of corona virus begins with the fusion and attachment to the host receptor angiotensin converting enzyme-2 (ACE2) through its spike glycoprotein. Upon invasion and after releasing its genomic material inside the host cell, the genome gets translated into several polypeptides followed by the translation into sixteen non-structural polyproteins (NSPs). The 3’ end of the sub-genomic RNA gets transcribed into several structural proteins which are assembled on the endoplasmic reticulum (ER) and golgi body to release the viral progenies outside the plasma membrane [16].

Fig. 1

Life cycle of corona virus inside the host cell. ACE2:Angiotensin Converting Enzyme 2; 3’UTR:Untranslated region; pp1a/pp1ab: polyproteins; S: spike glycoprotein; M: membrane protein; E: envelope protein; N: nucleocapsid protein.

Despite of current antiviral strategies including the anti-viral agents and vaccines to combat the deadly coronavirus, these could not prove themselves effective against the virulent disease. In order to generate the scopes for effective inhibitors of SARS CoV-2, several drugs have been repurposed and implemented to combat this ailment, though they have not shown any reliable results necessitating the demand for some potential inhibitors [17,18]. In view of this fact, several vaccines containing live attenuated, RNA-DNA derived, viral vectored, protein contained and inactivated vaccines have been developed with positive results against the catastrophic virus, SARS CoV-2 [19]. In contrary, the vaccines under phase III trials have also been noted with several side effects after the vaccination such as headache, fatigue, chills, diarrhoea, arthralgia, myalgia, redness, itching, swelling, etc. [20] To exemplify, a nucleotide analog, remdesivir was declared under emergency use authorization by United States food and drug administration (USFDA) and was evaluated under in vivo and in vitro studies for the inhibition of viral replication inside the host body [21]. Despite of its efficacy, it has not been found effective solution to tackle the race for the inhibition of the deadly virus due to its impotency in patients with immunodeficiency and heart failure, hypotension, respiratory dysfunction, as well as renal impairment [22].

The conventional approach to the drug design and development may take plenty of resources such as cost, time and manpower with no guarantee to afford the desired candidates with efficacy and potency. Therefore, the computer assisted drug design (CADD) have been preferred over several years above the traditional approaches reducing the cost-burden by benefiting the feasibility and improved results after removing the gap between chemical and biological science [23]. Among these CADD tools, the use of molecular modeling should not be undermined as they provide significant interactions of identified hits against their biological targets to understand their mode of action [24], [25], [26], [27], [28]. Recently, we have performed in silico screening of significant phytochemicals against RNA-dependent RNA Polymerase (RdRp) and main protease (Mpro) using molecular docking followed by their druggability using ADMET assay and stability of their complexes using molecular dynamics (MD) simulations in search of effective SARS CoV-2 inhibitors [29]. Motivated by the potential of computational chemistry and in search of potent inhibitors of SARS CoV-2, herein we have performed the virtual screening of total 267,324 ligands from 2020–01 Asinex EiteSynergy (91,473) and BioDesign (175,851) libraries using AutoDock Vina against RdRp using molecular docking followed by their further assessment using rules for druggability, ADMET assay and MD simulations.

Results and discussion

Analysis of protein structure

Xu et al. have reported the cryo-EM structure of the RdRp from SARS-CoV-2, co-crystallized with remdesivir having the resolution of 2.8 Å (PDB ID: 7BV2) [30,31]. In this complex, partial double-stranded RNA template was implanted into the central channel of the RdRp where remdesivir was covalently fused into the primer strand at the first replicated base pair and terminated the chain elongation. The core component of the viral replication complex is the non-structural proteins (nsp12) of the RdRp [32,33].

Molecular docking

The molecular docking has been performed against the selected protein (PDB ID: 7BV2) [30] using AutoDock Vina [34] to evaluate the binding mode of ligand and interactions in the active site. The required 267,324 ligands have been obtained from 2020–01 Asinex EiteSynergy (91,473) and BioDesign (175,851) library from the Asinex online database, a part of PubChem database consists of small molecules, which can be easily synthesized in organic chemistry labs [35]. The preprocessed protein and ligands in PDBQT format have been used for the molecular docking process. The binding energy of the hit molecules (2–22) have been compared with the co-crystallized ligand remdesivir (1, Table 1 , Entry 1). Further, these hits (2–22) have been found with better binding energy as compared to other promising RdRp inhibitors such as ribavirin (23), penciclovir (24), favipravir (25), [36] molnupiravir (26) and sofosbuvir (27) [37]. Total twenty-one identified macrocyclic hits have been found to possess binding energy lesser than -12.5 Kcal/mol and presented in Table 1 and their 2D-possess with the formed interactions have been presented in Figs. S1–S4 (see supporting informartion).

Table 1

The identified hits (2–22) having more binding energy than clinically used remdesivir (1) and other RdRp inhibitors (23–27).

Compound No.	Dataset ID	2D-Structure	Binding energy (Kcal/mol)
1	Remdesivir	Image, table 1	-9.2
2	LAS 51620435	Image, table 1	-13.1
3	LAS 51620429	Image, table 1	-13
4	LAS 51609341	Image, table 1	-13
5	LAS 51620411	Image, table 1	-12.9
6	LAS 51609377	Image, table 1	-12.8
7	LAS 51620382	Image, table 1	-12.8
8	LAS 51624260	Image, table 1	-12.8
9	LAS 51624264	Image, table 1	-12.8
10	LAS 51624295	Image, table 1	-12.8
11	LAS 52131039	Image, table 1	-12.8
12	LAS 51609335	Image, table 1	-12.7
13	LAS 51614073	Image, table 1	-12.7
14	LAS 51609336	Image, table 1	-12.6
15	LAS 51620402	Image, table 1	-12.6
16	LAS 51620430	Image, table 1	-12.6
17	LAS 51620433	Image, table 1	-12.6
18	LAS 51624245	Image, table 1	-12.6
19	LAS 51609330	Image, table 1	-12.5
20	LAS 51620424	Image, table 1	-12.5
21	LAS 51624284	Image, table 1	-12.5
22	LAS 51624301	Image, table 1	-12.5
23	Ribavirin	Image, table 1	-7.4
24	Penciclovir	Image, table 1	-7.2
25	Favipravir	Image, table 1	-6.7
26	Molnupiravir	Image, table 1	-10
27	Sofosbuvir	Image, table 1	-8.7

The 3D-poses of identified ligands 2–7 bound at the active site of RdRp have been generated using Biovia Discovery Studio 2020 [38] and presented in Fig. 2 a–f. Compound 2 formed conventional hydrogen bonds (HB) with the amino acid Asn496 from chain A and nucleotides Ade11 and Ade13 of chain P, Ade14, Ade15 and Ade19 of chain T. Further, it formed van der Waals interactions with Ile548 and Ade14 of chain P, π-anion interaction with Urd18 (chain P), π-alkyl interaction with Ala547, Ade13 and Ade14 (chain T, Fig. 2a). Additionally, the presence of sulfonyl group in the structure of 2 could lead the formation of π-sulfur interaction with Ade15 of chain P. Pyrimidine derivative 3 made conventional HB with Ade14 of chain P and Ade11 of chain T, C-H bond with Lys500 and Urd12 of chain T, π-π bond with Phe441, π-alkyl bond with Ala547 and Ile548 along with π-anion bond with Urd17 of chain P (Fig. 2b). Compound 4 formed conventional HB with the amino acid, Asn496 and nucleotides, Ade13 and Ade14 of chain P and Ade14 and 15 of chain T. Additionally, it formed π-π bond with Urd12 of chain T (Fig. 2c). By the replacement of 2-methylpyrimidine ring of compound 3 with pyridine ring in compound 5 retained almost similar interactions including van der Waals bond (Lys500), π-π bond (Phe441, π-alkyl bond with Ala547 and Ile548), π-anion bond (Urd17 and Urd18 of chain P) and conventional HB (Ade14 of chain P and Ade11 of chain T, Fig. 2d). Compound 6 found to form conventional HB with Ade14, Gua16 and Ade19 of chain P, Ade13 of chain T, Lys500 and Asn497. It also formed C-H bond with Asn497, π-anion bond with Ade19 of chain P and π-alkyl bond with Ala547 and Urd18 of chain P (Fig. 2e). N-Methylsulfonyl piperidine substituted compound 7 interacted at the active site of RdRp by forming π-alkyl bond (Ala547, Ile548 and Ade13 of chain T), π-sulfur bond (Ade14 of chain P), π-anion bond (Urd17 and Urd18 of chain P), C-H bond (Lys500) and conventional HB (Ade11 of chain T and Ade14 and 15 of chain P, Fig. 2f). In general, these interactions with the ligands have supported our hypothesis on the tight binding of identified ligands into the active site of protein.

Fig. 2

The 3D-interactions of ligands 2 (a), 3 (b), 4 (c), 5 (d), 6 (e) and 7 (f) complexed with RdRp. Poses have been generated using Biovia Discovery Studio 2020 [38].

2,4-Disubstituted pyridine derivative 8 has been found with binding energy of -12.8 Kcal/mol (Table 1) due to interaction with RdRp including π-alkyl bond with Ade19 of chain P and conventional HB with Ade11 of chain T, Ade14 and Ade19 of chain P (Fig. 3 a). 1,3-Dimethyl-1H-pyrazol-4-yl sulfonyl derivative compound 9 (Fig. 3b) found to bind tightly through the formation of π-anion bond with Urd13 of chain P, non-polar C-H bond with Urd18 of chain P, π-sulfur bond with Ade19 of chain P and HB with Ade11 (chain T) Ade14 and Ade19 (chain P). Interactions including C-H bond (Lys500), π-alkyl bond (Ala547), π-sulfur bond (Ade14 of chain T), π-anion bonds (Asp845 and Urd18 of chain P) and conventional HB (Ade13, Ade14 of chain T and Gua16, Ade15 of chain P) have been found with the docked complex of 1-sulfonyl pyrrolidine-3-carboxamide 10 against RdRp (Fig. 3c) with similar binding energy. Apart from the interactions including conventional HB (Ade14 of P chain, Ade13 and Ade14 of T chain), two π-π bonds (Urd12) and π-alkyl bond (Urd10 from chain T), N-cyclopenyl indole derivative 11 formed halogen interactions between chlorine atom present in it with Cyt15 and Ade14 of chain T (Fig. 3d). The replacement of pyrimidine-2,4(1H,3H)-dione-6-methyl in ligand 6 with pyridin-2-yl-acetamide in 12 has been resulted in almost similar binding energy (-12.8 and -12.7 kcal/mol, respectively) along the with formation of π-anion bond with Ade19 of chain P, π-alkyl bond with Ala547 and Urd18 of chain P, conventional HB with Ade14, Ade15 and Ade19 belonging to chain P (Fig. 3e). Several π-alkyl interactions are found in docked complex of compound 13 with amino acids Phe441, Ile548, Lys500, nucleotides Urd17, Gua16 (chain P) and Ade11, Urd12 (chain T). Further, ligand 13 formed conventional HB with Lys545, Arg836 and Ade14 of chain T and π-anion bond with Urd18 of T chain and π-σ bond with Ala547 (Fig. 3f).

Fig. 3

3-D interactions of ligands 8 (a), 9 (b), 10 (c), 11 (d), 12 (e) and 13 (f) with RdRp (PDB ID: 7BV2). Interactive poses of hits have been generated using Biovia Discovery Studio 2020 [38].

Compound 14, with binding energy of -12.6 kcal/mol, formed conventional HB with Ade19 of chain P and Ade11, Ade14 and Cyt15 of chain T, π-anion bond with Urd13 of chain P. Additionally, nucleotide Urd18 has found with π-π bond, π-sulfur bond and π-alkyl bond with the ligand 14 (Fig. 4 a). Compound 15, a metamer of compound 9, formed completely dissimilar interactions in comparison with 9 such as C-H bond with Lys500, π-sulfur bond with Ade13 of T chain, π-alkyl bonds with Ala547 and Ile548, π-anion bonds with Urd17 and Urd18 from chain P and conventional HB with Asn496, Ade11 and Ade13 of chain T (Fig. 4b). Next, 2,3-dihydrobenzo[b][1,4]dioxin-6-sulfonyl derivative 16 found to form π-sulfur bond with Ade14 of chain T, π-alkyl bond with Ile548, π-anion bond with Urd17 and Urd18 of chain P and HBs with the amino acid Asn496 and nucleotides Ade15 and Ade19 from P chain and Ade13 and 14 from T chain (Fig. 4c). Compound 17 being N-(4-flurophenylsulfonyl)pyrrolidine derivative formed two π-alkyl bonds (Ala547), π-sulfur bond (Ade17 of chain T), C-H bonds (Lys500 and Urd12 of chain T), π-anion bonds (Asp845, Urd17 and 18 from P chain), conventional HB (Ade13, Ade14 and Cyt15 of chain T). Further, 4-fluoro on phenyl ring interacted with Asn496 through halogen interactions (Fig. 4d). N-Ethyl 2-substitued benzo[d]imidazole derivative 18 interacted with the protein through C-H bonds with Gua16 and Urd17 of chain P and π-alkyl bond with Ade15, Gua16, Ade19 and Urd20 of chain P (Fig. 4e). 2-Morpholinyl-piperdin-4-yl compound 19 has been found to interact with π-anion bond with Ade19 of P chain, π-alkyl bond with Urd18 of P chain and conventional HB with Ade13 of chain T, Ade19 of chain P and Asn497 (Fig. 4f) with binding energy of -12.5 kcal/mol.

Fig. 4

3D interactions of ligands compound 14 (a), compound 15 (b), compound 16 (c), compound 17 (d), compound 18 (e) and compound 19 (f) with 7BV2. Poses have been generated using Biovia Discovery Studio 2020 [38].

Ligands 20–22 have been observed with similar binding energy of -12.5 kcal/mol (Table 1). Ligand 20 formed conventional HB and π-sulfur bond with Ade15 of chain P, C-H bonds with Lys500 and Ade14 of chain P, π-alkyl bonds with Ala547 and Ile548, π-anion bonds with Urd17 and 18 of chain P. One of the fluorine atom of phenyl ring interacted with Urd13 of chain P (Fig. 5 a). Compound 21, being the positional isomer of 8, made π-anion bond with Urd17 of chain P and conventional HB with Lys846, Arg858, and Ade11 of chain T and Ade14 of chain P. Additionally, two fluorine atoms of trifluoromethyl group in 21 were found to interact with Asp845 (Fig. 5b). α-Naphthyl-1H-pyrrolidonone derivative 22 formed conventional HB with Ade15 of chain P, C-H bond with Asn496, π-alkyl bonds with Ala547 and Urd18 of chain P and π-anion bonds with Arg836 and Ade19 of chain P (Fig. 5c) to accommodate on the active site of RdRp.

Fig. 5

3D interactions of ligands 20 (a), 21 (b) and 22 (c) with 7BV2 (RdRp). Poses have been generated using Biovia Discovery Studio 2020 [38].

Analysis of physicochemical properties and ADMET parameters

The identified macrocyclic hits have been analyzed for the several physicochemical properties (Table 2 ) such as molecular weight (MW), hydrogen bond donors (HBD), hydrogen bond acceptors (HBA), number of rotatable bonds (RB) and partition co-efficient (Log P) using SwissADME [39] and pkCSM [40,41]. All the hits have better physicochemical properties than clinically used remdesivir (1). Almost all the hits have molecular weight beyond 500 including remdesivir. No compound from the hits has found violating the criteria of hydrogen bond acceptors, hydrogen bond donors and rotatable bonds. Only six compounds (8, 11, 13, 18, 21 and 22) exceeded the maximum limit of Log P value. These exceeded properties may be tackled through novel sophisticated techniques during development of formulation with a view of improving the drug-likeness. Thus, most of the hits may have potential to be a better drug molecule in comparison with remdesivir (1, Table 2).

Table 2

Physicochemical properties of the macrocyclic hits (2–22) and marketed drug, remdesivir (1).a

Comp No.	MWb	HBDc	HBAd	RBe	Log Pf	No. of violationsg
1	602.58	4	12	14	2.31	2
2	625.82	1	8	6	3.78	1
3	589.77	1	6	6	4.36	1
4	528.64	0	5	3	3.80	1
5	574.75	1	5	6	4.66	1
6	494.58	2	5	3	2.01	0
7	625.82	1	8	6	3.78	1
8	572.66	1	7	4	5.24	2
9	585.76	1	7	4	3.05	1
10	571.70	1	7	4	3.76	1
11	560.08	1	5	3	7.05	2
12	504.62	1	5	4	4.06	1
13	513.07	1	4	2	6.12	2
14	562.72	0	5	3	4.48	1
15	585.76	1	7	4	3.05	1
16	570.70	1	7	5	3.54	1
17	571.70	1	7	4	3.76	1
18	571.75	1	4	4	5.19	2
19	532.67	0	5	3	3.94	1
20	585.73	1	7	4	4.15	1
21	572.66	1	7	4	5.24	2
22	553.69	1	4	3	5.11	2

Parameters calculated using pkCSM [40,41],

Molecular weight (≤500 Da),

Number of hydrogen bond donors (≤ 5),

Number of hydrogen bond acceptors (≤ 10),

Number of rotational bonds (≤ 10),

Partition coefficient in oil to water (≤ 5).

Table 3.

Table 3

Predicted absorption and distribution properties of the hits (2–22) and remdesivir (1).a

Comp No.	Ali log Sb	MRc	tPSA (Å2)d	CaCO2 permeabilitye	Human intestinal absorption (% absorbed)f	VDss (Human)g	Fraction unbound (Human)h	P-gp inhibition (yes/no)i
1	-6.01	150.43	213.36	0.51	69.44	-0.38	0.03	Yes
2	-5.89	182.34	130.18	0.99	98.37	-0.40	0.00	Yes
3	-5.57	178.39	104.73	1.00	97.33	-0.09	0.02	Yes
4	-5.37	157.77	84.74	1.49	94.72	0.28	0.00	Yes
5	-5.56	175.62	91.84	0.98	89.36	0.04	0.01	Yes
6	-3.95	143.15	115.57	1.04	78.89	-0.13	0.21	Yes
7	-5.89	182.34	130.18	0.90	92.22	-0.47	0.00	Yes
8	-6.56	161.6	74.77	0.79	90.74	-0.02	0.00	Yes
9	-5.05	168.38	122.22	0.95	86.35	-0.45	0.00	Yes
10	-5.75	161.52	104.4	0.95	95.03	-0.31	0.00	Yes
11	-6.79	165.08	64.96	1.03	91.72	0.55	0.13	Yes
12	-4.89	149.43	91.84	0.88	95.13	0.06	0.06	Yes
13	-7.47	149.63	71.78	0.65	89.23	0.46	0.00	Yes
14	-6.77	165.58	112.98	0.92	93.39	0.14	0.02	Yes
15	-5.24	168.38	122.22	1.02	86.64	-0.51	0.00	Yes
16	-5.78	158.12	119.62	1.13	87.27	-0.31	0.00	Yes
17	-5.94	161.52	104.4	1.02	95.32	-0.39	0.00	Yes
18	-6.66	178.16	79.7	0.74	90.82	0.61	0.02	Yes
19	-5.15	160.85	75.21	1.45	96.05	0.02	0.00	Yes
20	-6.31	166.33	104.4	1.00	100	-0.39	0.00	Yes
21	-6.56	161.6	74.77	0.77	90.62	-0.19	0.00	Yes
22	-6.19	171.7	78.95	0.76	93.19	0.48	0.00	Yes

Parameters calculated using pkCSM [40].

Aqueous solubility descriptor (≤ 0),

molar refractivity (≤ 155),

topological polar surface area (≤ 150 Å2),

Caco-2 cell permeability (log Papp in 10−6 cm/s > 0.09),

absorption (human, % > 30),

volume of distribution (human, log L/kg) (low if < − 0.15 and high if > 0.45),

fraction unbound, and

ability to inhibit the P-glycoprotein.

Rules for the drug-likeliness such as Lipinski's rule of 5 (Ro5), Ghose rule, Veber's rule, Egan's rule and Muegge's rule have been also studied to ensure the drug-likeliness of 1–22 and number of violations have presented graphically in Fig. 6 . According to the rule of 5, the molecule can serve as better drug candidate with the suitable properties including MW (≤ 500 Da), HBD (≤ 5), HBA (≤ 10), LogP (≤ 5) and RB (≤ 5) [42]. Not more than one violation of these parameters will reduce the druggability of the molecule. Most of the hits (except 8, 11, 13, 18, 21 and 22) have found to meet the desired properties suggesting their candidature for the anti-viral drug. As per the Veber's rule, compound should have the polar surface area less than 140 Å2 to have the better oral bioavailability [43]. None of hits violated Veber's rule suggesting their good oral bioavailability. All hits except 11 qualified the criteria for oral bioavailability laid down by the Egan rule with desired tPSA (0–132 Å2) and log P (-1–6) [44]. Muegge's filter includes several parameters including MW (200–600), lipophilicity (XLOGP3, -2–5), tPSA (≤ 150), cyclic rings (≤ 7), carbon atoms (> 4), heteroatoms (> 1), RB (≤ 15), HBA (≤ 10), and HBD (≤ 5) [45] and only few compounds (2, 7, 8, 11, 13, 18 and 21) violated with note more than one parameter to qualify it. Ghose rule indicated the desired requirements of MW (160–480 Da), Log P (-0.4–5.6), MR (40–130) and atoms (20–70) to be the good drug candidate [46] and all the hits meet the suggested properties. However, as compared to Remdesivir, the identified hits may have potential to act as a good drug candidate in the future as they obey the rules for druggability and drug-likeness.

Fig. 6

Result of the hits against different violations for the rules of druggability and drug-likeness.

Next, we have studied the various ADMET (absorption, distribution, metabolism, excretion and toxicity) parameters like water solubility (Ali log S), molar refractivity (MR), topological polar surface area (tPSA), CaCO2 cell permeability, intestinal absorption, volume of distribution (Vd), unbound fraction of drug and ability to inhibit the P-glycoprotein substrate for the identified hits. All the molecules have satisfied the criteria for the Ali log S, tPSA and CaCO2 cell permeability. Most of the hits with sufficient human intestinal absorption indicated their better oral bioavailability. All the hits also qualified the criteria for Vd along with the inhibition of the P-glycoprotein substrates. Further, we analyzed the additional parameters to predict the metabolism, excretion and toxicity profile of remdesivir and these hits (Table 4 ). Most of the molecules inhibited CYP2D6 but not CYP3A4 along with medium renal clearance except compound 16 and 22 (with low clearance). All the hits except hits 4, 5, 18 were found as the substrates for renal uptake transporter in proximal convoluted tubule (OCT2). All the hits were found non-cytotoxic (hERG cell line), non-mutagenic (AMES toxicity), non-toxic to dermis (skin sensitization) and found safe through prediction of oral rat acute and (LD50) chronic toxicity (LOAEL). Thus, these macrocyclic hits (2–22) might be good and successful drug candidates in the future.

Table 4

Metabolism, excretion and safety parameters of remdesivir (1) and identified hits (2–22).a

Comp No.	CYP2D6 inhibitorb	CYP3A4 inhibitorc	CLTd	Renal OCT2 substratee	AMES toxicityf	hERG I toxicityg	LD50h	LOAELi	Skin sensitizationj
1	No	No	0.16	No	No	No	2.25	2.27	No
2	No	Yes	0.56	No	No	No	3.09	1.54	No
3	No	Yes	0.50	No	No	No	3.20	1.33	No
4	No	Yes	0.27	Yes	No	No	2.88	-0.08	No
5	No	Yes	0.49	Yes	No	No	3.34	-0.45	No
6	No	Yes	0.42	No	No	No	2.40	0.92	No
7	No	Yes	0.54	No	No	No	2.95	1.51	No
8	No	Yes	0.11	No	No	No	2.70	0.80	No
9	No	Yes	0.36	No	No	No	2.97	1.78	No
10	No	Yes	0.37	No	No	No	2.80	1.44	No
11	No	Yes	0.55	No	No	No	3.24	-0.19	No
12	No	Yes	0.37	No	No	No	2.65	1.60	No
13	No	Yes	-0.09	No	No	No	2.49	0.71	No
14	No	Yes	0.05	No	No	No	2.92	-0.42	No
15	No	Yes	0.40	No	No	No	3.11	1.63	No
16	No	Yes	0.43	No	No	No	3.05	1.69	No
17	No	Yes	0.40	No	No	No	2.98	0.94	No
18	No	Yes	0.44	Yes	Yes	No	2.37	1.47	No
19	No	Yes	0.40	No	No	No	3.18	-0.12	No
20	No	Yes	0.41	No	No	No	3.03	0.97	No
21	No	Yes	0.17	No	Yes	No	3.00	-0.08	No
22	No	Yes	-0.03	No	No	No	2.46	1.79	No

Parameters calculated using pkCSM [40],

ability to inhibit CYP2D6 enzyme,

ability to inhibit CYP3A4 enzyme,

total renal clearance; high (>1 mL/min/kg), medium (> 0.1 to < 1 mL/min/kg) or low (≤ 0.1 mL/min/kg),

ability to inhibit renal OCT2 substrate;

AMES toxicity;

hERG I toxicity;

oral rat acute toxicity (LD50);

oral rat chronic toxicity (LOAEL);

skin sensitisation.

The Brain Or IntestinaL EstimateD permeation method (BOILED-Egg) provide the in silico estimation of accessibility of molecules to gastrointestinal (GI) tract and blood-brain barrier [47]. Hence, we also become interested to predict the permeability of the identified macrocyclic hits using SwissADME [39]. The boiled egg model revealed all the hit molecules possessed satisfactory GI absorption along with inhibition of the P-glycoprotein, a protein responsible for efflux of drugs from cells (Fig. 7 ). All the compounds (except 19) were found with probability of permeation to BBB and henceforth, the least chances of neurotoxicity of CNS with the proposed compounds 2–18 and 20–22 may be approximated.

Fig. 7

BOILED egg model of hit molecules (2–22) generated with the help of SwissADME. The yellow and colorless regions indicate BBB and GI permeability, respectively. The blue circles indicate inhibition of P-glycoprotein.

MD simulation

To study the stability of ligand into the binding site of protein, molecular dynamics (MD) simulations provide the better understanding of ligand through several statistical parameters [48,49]. Several scientific literatures support the reliability and satisfactory stability for the MD run (≤ 10 ns) through in silico endeavor [50], [51], [52], [53], [54], [55]. In a similar line of approach, the hit complexes obtained from the evaluation of ADMET and drug-likeliness properties, were subjected for MD simulations using GROMACS 2020.1 to assess the stability of the ligands (2 and 3) in the active site of docked complex at various time points up to 10 ns [56,57].

Next, we subjected the complex of compound 2 to MD simulation studies and the graphical representation of plots of statistical parameters has been presented in Fig. 8 . The ligand-receptor complex of 2 with RdRp was found with root mean square deviation (RMSD) value with an average of 0.757 nm (Fig. 8a) for the ligand and 0.267 nm (Fig. 8b) for the protein which indicate the stability of the ligand with least deviation in the active site of remdesivir. The radius of gyration (RoG) for the same complex ranging from 2.94 to 3.01 nm with an average of 2.972 nm (Fig. 8c) indicated the compactness of the complex. The plots of solvent assessable surface area (SASA) and hydrogen bonds has been presented in Fig. 8d and Fig. 8e, respectively. The surface area accessed by the solvent molecules was found within the satisfactory range from 447–461 nm2 with an average of 454.72 nm2 (Fig. 8d). From the plot of number of HB vs simulation time, maximum four HB were found between the ligand and receptor within the time period of 10 ns (Fig. 8e).

Fig. 8

The schematic plots of RMSD-L (a), RMSD-P (b), RoG (c), SASA (d) and HB (e) for the complex of compound 2 with RdRp.

The satisfactory stability of the ligand with least deviation in the active site of complex has been evident from RMSD value with an average of 0.942 nm (ligand, Fig. 9 a) and 0.270 nm (protein, Fig. 9b) for the ligand-receptor complex of compound 3 with RdRp. The RoG for the same complex was found in the range of 2.94–3.00 nm with an average of 2.966 nm (Fig. 9c) indicating their compactness. The surface area accessed by the solvent/water molecules was found within the range from 446–460 nm2 with an average of 454.215 nm2 (Fig. 9d). From the plot of number of HB vs time, maximum two HBs were found between the ligand and receptor within the time period of 10 ns (Fig. 9e). The electrostatic (coulombic short-range, Coul-SR) and van der Waals/hydrophobic (LJ-SR) energies have been represented in Table 5 for the complexes of identified hits 2 and 3 with RdRp. These complexes have been observed to be stabilized by the van der Waals or hydrophobic interactions significantly over the electrostatic or coulombic contributions.

Fig. 9

The schematic plots of RMSD-L (a), RMSD-P (b), RoG (c), SASA (d) and HB (e) for the complex of compound 3 with RdRp.

Table 5

The electrostatic steric interactions (KJ/mol) for the complex of respective compounds (2 and 3) with RdRp protein.

Compd. No.	Energy (KJ/mol)
	Electrostatic interaction (Coul-SR)	van der Waals/Hydrophobic interactions (LJ-SR)
2	-48.7869 ± 4.6	-115.263 ± 9.1
3	-32.4944 ± 8.6	-108.078 ± 8.1

Conclusion

In conclusion, the in silico based virtual screening of total 267,324 ligands from 2020-01 Asinex EiteSynergy (91,473) and BioDesign (175,851) libraries against RdRp using AutoDock Vina lead to discovery of twenty-one macrocyclic hits (2–22) with better binding energy than remdesivir (1), clinically used SARS CoV-2 inhibitor. The several interactions between ligand and protein helped to understand their binding mode into the active site of protein. The drug-likeliness and druggability studied using several rules for druggability and analysis of ADMET profiles of these hits (2–22) revealed the better bioavailability and safety profile than 1. MD simulation studies up to 10 ns using GROMACS 2020.1 of the top two hits (2 and 3) indicated their stability into the active site of complex. Further, the structural modifications of these molecules may have the scope with significant inhibition of SARS CoV-2 in future.

Experimental section

Collection and preparation of data

The 3-dimentional structure of the SARS-CoV-2 RdRp evaluated through cryo electron microscopy (resolution: 2.8 Å) having antiviral remdesivir as a co-crystallised inhibitor (PDB ID: 7BV2) was obtained from RCSB protein data bank [30]. The protein was prepared for further molecular modelling using AutoDock Vina [58]. All the water molecules were deleted along with addition of polar hydrogen atoms and Kollman charges to complete the protein in terms of polarity and charges, respectively. The generation of the receptor grid box around the macromolecule of co-crystallised ligand remdesivir was achieved with the size of 90.768 × 99.896 × 99.788 Å (x, y and z) and the co-ordinates of centers 70 × 70 × 70 Å (x, y, and z). The prepared protein has been kept in PDBQT format which has been further used for the molecular docking. The present docking protocol has been validated by docking of remdesivir using the same and observed with RMSD of less than 1.

The Asinex BioDesign library 2020.1 having 175,851 molecules and Asinex EliteSynergy having 91,473 molecules was downloaded from Asinex database [35]. All the ligand molecules have been optimized and converted into PDBQT format using OpenBabel and further used for molecular docking process [59].

The multi-ligand molecular docking was performed using the prepared protein using AutoDock Vina to reveal the binding affinity and interactions of the ligands with RdRp protein (PDB ID: 7BV2) [30]. All the ligands (267,324) were docked on the site using the optimized grid box and ten poses per ligand were generated. The docked pose of remdesivir was compared with the co-crystallized inhibitor remdesivir to validate the present docking protocol. The docking interactions of hits were visualized and analyzed using Biovia Discovery studio 2020 [38]. The 2D interactions of docked compounds (2–22) with RdRp have been presented in supporting information (Figs S1–S4).

Analysis of physicochemical properties, drug-likeliness and ADMET parameters

Various physicochemical properties like molecular weight, number of hydrogen bond donors, number of hydrogen bond acceptors, oil to water partition co-efficient and number of rotational bonds of the selected hits (2–22) and remdesivir (1) were studied using Swiss ADME [39] and pkCSM [40,41]. Several rules such as Lipinski's rule of 5, Ghose rule, Veber's rule, Egan's rule and Muegge's rule were applied to the selected hits using the SwissADME to ensure their druggability and to identify numbers of violated parameters. The ADMET properties like water solubility, molar refractivity, topological polar surface area, CaCO2 cell permeability, intestinal absorption, volume of distribution, fraction unbound, total renal clearance, hepatotoxicity and ability to inhibit the P-glycoprotein were also studied using pkCSM.

Molecular dynamics (MD) simulation

Two molecules (2 and 3), obtained from the manual analysis of the hits found from the results of molecular docking and ADMET analysis, were selected for MD simulation using GROningen MAchine for Chemical Simulations (GROMACS) 2020.1 [56,57] software. CHARMM36 (Chemistry at Harvard Macromolecular Mechanics) as an all atom force field [60] and CHARMM General Force Field (CGenFF) server to retrieve the topology of (2 and 3) were used [61,62]. After solvation (TIP3P water model), neutralization (Na+ and Cl− ions), equilibration [canonical (NVT) and isobaric-isothermic (NPT) ensemble for 100 ps], the complexes were subjected to molecular dynamics run for 10 ns. By analyzing the plots of temperature, pressure and energy equilibria of the trajectories (Figs S5 and S6, see supporting information), we found the system satisfactorily equilibrated.

CRediT authorship contribution statement

Prinsa R. Nagar: Methodology, Data curation, Writing – original draft, Writing – review & editing. Normi D. Gajjar: Methodology, Data curation, Writing – original draft, Writing – review & editing. Tejas M. Dhameliya: Data curation, Conceptualization, Supervision, Writing – original draft, Writing – review & editing.

Declaration of Competing Interest

The authors declare no competing financial interest.