Determination of potential inhibitors based on isatin derivatives against SARS-CoV-2 main protease (mpro): a molecular docking, molecular dynamics and structure-activity relationship studies.

SARS-COV-2, the novel coronavirus and root of global pandemic COVID-19 caused a severe health threat throughout the world. Lack of specific treatments raised an effort to find potential inhibitors for the viral proteins. The recently invented crystal structure of SARS-CoV-2 main protease (Mpro) and its key role in viral replication; non-resemblance to any human protease makes it a perfect target for inhibitor research. This article reports a computer-aided drug design (CADD) approach for the screening of 118 compounds with 16 distinct heterocyclic moieties in comparison with 5 natural products and 7 repurposed drugs. Molecular docking analysis against Mpro protein were performed finding isatin linked with a oxidiazoles (A2 and A4) derivatives to have the best docking scores of -11.22 kcal/mol and -11.15 kcal/mol respectively. Structure-activity relationship studies showed a good comparison with a known active Mpro inhibitor and repurposed drug ebselen with an IC50 value of -0.67 μM. Molecular Dynamics (MD) simulations for 50 ns were performed for A2 and A4 supporting the stability of the two compounds within the binding pocket, largely at the S1, S2 and S4 domains with high binding energy suggesting their suitability as potential inhibitors of Mpro for SARS-CoV-2.

Introduction

During December 2019, a novel coronavirus named “Severe Acute Respiratory Syndrome Corona Virus 2 (SARS-CoV-2), caused an outbreak of a respiratory disease in the city of Wuhan, capital of the Hubei province in China and since has been spreading globally (Wu et al., 2020b; Zhou et al., 2020) and the pulmonary disease caused by SARS-CoV-2 is known as COVID-19. The virus was named SARS-CoV-2 (Gorbalenya et al., 2020) due to the resemblance of the RNA genome with the previous coronavirus in 2003, i.e. the SARS-CoV. Both viruses belong to the clade b of the genus Betacoronavirus (Wu et al., 2020b; Zhou et al., 2020). Since the outbreak, the incredible capability of the virus to spread from human-to-human led to an exponential growth in the number of patients with world-wide spreading of the virus. According to the World Health Organization (WHO), the virus infects most of the patients with mild symptoms such as fever, cough and difficulty in breathing that are cured with proper care. However, in older patients or in patients with underlying health conditions, the disease can progress into fatal pneumonia and acute respiratory failure (Chen et al., 2020a; Zhou et al., 2020). The progress of the virus was due to its novel transfer from human to human through cough, sneeze, touch from an infected person i.e. through the common droplet infection (Chen et al., 2020c; Wu et al., 2020b; Zhou et al., 2020). As of July 2nd, 11,000,000 cases and over 525,000 deaths worldwide exerting an enormous social impact and severely effecting world economic status, as previous pandemics (Keogh-Brown & Smith, 2008).

Currently, many researchers have been working to develop specific drugs through repurposing of drugs for diseases like HIV, Ebola, malaria etc. (Fischer et al., 2020; Wang et al., 2016; Xu et al., 2020b; Zhang et al., 2020a). Still these significant developments have yet to achieve prevention of disease spreading. A key strategy for treatment viral infection has been inhibiting its main protease, such inhibitors are in clinical use as effective treatments for human immunodeficiency virus (HIV) and hepatitis C (Ghosh et al., 2016; Yang et al., 2006).

The SARS-CoV-2 virion (Scheme 1) has four structural proteins (Wu et al., 2020a) spike proteins (S), Envelope (E), Membrane protein (M) and nucleocapsid proteins (N) which holds the RNA genome, similar to the other coronaviruses. The Spike protein and Mpro (also called as 3CLpro) are the two main sites for drug development for SARS-CoV-2. Mpro plays a key role in mediating viral replication and transcription (Liu & Wang, 2020; Xu et al., 2020b; Jin et al., 2020b; Xu et al., 2020b). Mpro operates at 11 cleavage sites on the large polyprotein 1ab (replicase 1ab, ∼790 kDa), inhibiting the activity of this enzyme directly blocks viral replication (Jin et al., 2020b). As there are no human proteases with a homolog of Mpro, it is an ideal target for drug design as inhibitors are less likely to be toxic for humans (Zhang et al., 2020c). Thus, several studies describing Mpro inhibitors have been recently published (Chen et al., 2020b; Fan et al., 2020; Jin et al., 2020b; Liu & Wang, 2020; Raugi et al., 2016; Xu et al., 2020b).

Scheme 1.

Schematic representation of the structure of SARS-CoV-2.

The crystal structure of SARS-CoV-2 Mpro was recently solved and deposited to the protein data bank (PDB-6LU7) enabling rational design of specific inhibitors (Bzówka et al., 2020). There is a strong resemblance between SARS-CoV-2 and SARS-CoV-1, thus, many researchers were motivated to develop drugs based on their activities of the SARS- CoV-1 (Dai et al., 2020; Xu et al., 2020b). While some of the Mpro inhibitors are in preclinical or early clinical stages, there have been no approved treatment yet (Li et al., 2020; Pillaiyar et al., 2016). Hence, we were encouraged to perform a virtual screen of inhibitors that have a wide structural diversity and include a variety of structures to generate a large pool of potential Mpro inhibitors.

According to our literature survey, the presence of heterocyclic moieties was found in most of the potential drugs designed inhibitors for SARS-CoV-2 Mpro. Here, we have screened a wide variety of heterocyclic compounds including 118 previously studied anticancer, anti-bacterial and anti-microbial compounds in comparison to the published Mpro inhibitor N3 along, 5 natural products and 7 repurposed drugs.

Selection of candidate drugs that were screened

The 7 repurposed drugs are presently used by various researchers and medical institutes for the treatment of COVID-19, like Chloroquine (Cortegiani et al., 2020; Mallikarjuna et al., 2020), Hydroxychloroquine (Gautret et al., 2020), Remdesivir (Scavone et al., 2020), Favipiravir (Shiraki & Daikoku, 2020), Niclosamide (Xu et al., 2020a), Ebeselen (Jin et al., 2020a) and Eidd-2801 (Sheahan et al., 2020). Five natural compounds like Ursodeoxychloric acid, Quercetin, Kaempferol, Pinocembrin and Rutin are also included in this study. Chloroquine and Hydroxychloroquine are known drugs for treatment of malaria (Ben-Zvi et al., 2012; Bhattacharjee, 2016; Cortegiani et al., 2020). Remdesivir and Favipiravir are antiviral drugs potentially against several viruses (Scavone et al., 2020; Shiraki & Daikoku, 2020).

The broad spectrum antiviral drug Remdesivir, developed by the Gilead Sciences and the anti malaria drug Chloroquine were recently studied in an in vitro research work by G. Xiao et al., to be effective inhibitors for the 2019-nCoV (Wang et al., 2020). Chloroquine was also previously predicted to be an inhibitor for the SARS-CoV-1 (De Clercq, 2006). The other malaria drug Hydrochloroquine was also found to be effective in inhibiting SARS-CoV-2 infection in vitro and attenuate inflammatory response as given in a recent study (Liu et al., 2020). An experimental antiviral drug developed for the treatment of influenza Eidd2801 was also repurposed for SARS, MERS and SARS-CoV-2, still in pre-clinical trials though very fruitful results are not expected (Sheahan et al., 2020). Niclosamide, a drug designed to treat tapeworm infestations, is known for its broad antiviral activity, recent studies are performed to understand its activity against SARS-CoV-2 (Xu et al., 2020a). The natural compound Ursodeoxychloric acid is a mammalian bile acid found in bear and is used in treatment of cholesterol (Einarsson, 1994). Quercetin, Kaempferol, Pinocembrin are natural flavonoids found in several fruits and vegetables and used as antioxidants (Wang et al., 2018). Pinocembrin is also found in honey (Rasul et al., 2013). Rutin is a bioflavonoid, found in many citrus food and have powerful antioxidant properties (Enogieru et al., 2018).

Many of the screened drug molecules for our docking studies were taken from previously studied drugs, which were designed for anticancer activities. Some of them were also known for their antimicrobial and antibacterial activities (Afifi et al., 2017; Agarwal et al., 2016; Beedie et al., 2015; Bhatt et al., 2020; Blair et al., 2007; Fischmann et al., 2008; Gerova et al., 2016; Ivanova et al., 2007; Kok et al., 2008; Lane et al., 2001; Mendel et al., 2003; Pessoa et al., 2010; Schelman et al., 2011; Yadagiri et al., 2015; Zhang et al., 2011). Many of these drugs are already FDA approved drugs like Sunitinib (Bhatt et al., 2020) and many are in their pre-clinical trials like Zibotentan (Schelman et al., 2011).

Compounds based on oxidiazoles and isatin moieties, well known for vigorous research for anticancer, antiparasitic activities and such type of compounds were also screened on the basis of their structure similarities (Bhatt et al., 2020; Ivanova et al., 2007; Lakshmithendral et al., 2019; Mendel et al., 2003; Serafim et al., 2017; Venkatarao et al., 2019; Zhang et al., 2011). Other heterocyclic moieties with azo based structures and antibacterial, antitumor or anticancer activities were also screened (Bhatt et al., 2018; Jha & Ramarao, 2018; Tamokou et al., 2016). Antimicrobial compounds based on benzylpyrrolidine groups have also been screened for our cause (Sreekanth & Jha, 2020). As for comparison our in-silico screen included the very recent virtually screened 11 compounds that showed inhibitors of SARS-CoV2 (Fischer et al., 2020). In addition we selected compounds with a new heterocyclic ring system, piperdin-4-one semicarbazones, that exhibits moderate antibacterial and antifungal activities (Anand et al., 2019). A recent study by Zhang et al. (2020b) showed that α-keto amide inhibitors are very important inhibitors for MERS and predicted it to be better for SARS-CoV-2 as well. We have chosen the best structure (A118) for our docking analysis.

Our docking analyses showed high interaction of 12 novel heterocyclic drug moieties with SARS-CoV-2 main protease Mpro that can serve as potential novel drug candidates. Further Molecular dynamic (MD) analyses of the two best docked compounds A2 and A4, both derivatives of isatin and oxadiazole, showed that the interactions for these compounds are mostly governed by the hydrophobic and electrostatics interactions. The electrostatic charges and Molecular electrostatic potentials were also calculated for these compounds using Density Functional Theory (DFT) to understand their charge electrostatic charge, shape and the surface of the ligands. These calculations were done to predict the electronic distribution of the ligand which is one of the most important factor for analysing the binding between the protein and the ligand. Drug like properties of the compounds, ADME profile analysis, were also performed for these compounds.

This study suggests two compounds that bind Mpro extremely potently in-silico and may serve as the basis for potent inhibitors of COVID-19. Our calculations also shed light on the shape, electrostatic charge and surface of the binders showing adequate ADME properties.

Result and discussion

For our virtual screen of inhibitors of Mpro we selected structures of a library of newly designed compounds, A1-A118, in addition to natural products used as antiviral, anticancer, antibacterial drugs and COVID-19 drugs in clinical trials (Figure S1), as discussed in the introduction. All the structures were optimized with Gaussian 16 software at B3LYP/6-31 + G(d) level of theory (Gibbs free energy and the zero point energies are also given in Figure S1). Drugs used in clinical trials with specific crystal structures were taken directly without optimizing. All the structures were subjected to the docking analysis with the main protease Mpro of the SARS-CoV-2, PDB id: 6LU7. The free energy of binding (kcal/mol) and inhibitory constant (Ki) for these compounds are given in Table 1. The active site of Mpro is consists of five sub-sites (S1, S2, S1’, S2’, S4). S1 site consists of residues (Phe140, Leu141, Gln166, His163, Met165 and His172). The bulky hydrophobic S2 subsite (Met49 and Asp187, Glu189) (Wang et al., 2016), while S1’ (Thr25, Leu27, Cys38, Pro-39, Val42, and Cys145) is essential for catalysis (Xue et al., 2008), S2’ (Thr-26, Asn-28, Tyr-118, Asn-119, and Gly-143) and flexible S4 site (Leu167, Gln192) in nature (Mittal et al., 2020; Wang et al., 2016).

Table 1.

The calculated Ki, Free binding energy and bond length of hydrogen bonds and π-π interaction of selected FDA approved drugs, natural products containing heterocyclic moieties compared with the top 12 compounds with the best binding affinity to SARS-Co-V-2 Mpro.

S. No.	Compound(Pubchem ID)	InhibitoryConstant (Ki)	Free binding energy (kcal/mol)	H-bond and π-π interactions(bond length in Å)
1	N3 (native drug)	307.96 nm	−8.88	Cys145 (1.91 Å)
2	Ebselen (3194)	18.66 μm	−6.45	Glu166 (2.18 Å), Cys145
3	Favipiravir (492405)	378.76 μm	−4.67	–
4	Remdesivir (121304016)	676.24 nm	−8.42	Glu166-NH (2.09 Å), Glu166-OH (1.92 Å)
5	Chloroquine (2719)	7.27 μm	−7.01	Glu166 (2.04 Å), Arg188 (1.78 Å)
6	Hydroxychloroquine (3652 )	3.80 μm	−7.39	His164 (1.88 Å), Ser144 (2.21 Å), Leu141 (1.69 Å)
7	EIDD-2801 (145996610)	8.29 μm	−6.93	His41 (2.02 Å), Gly143 (2.21 Å), Ser144 (2.02 Å),Cys145 (2.11 Å), His 163 (2.14 Å),
8	Niclosamide (4477)	2.66 μm	−7.60	Gly143(1.76 Å), Cys145(2.16 Å),Glu166 (2.068 Å),
9	Ursodeoxycholic acid (31401)	123.64 nm	−9.42	Leu141 (2.22 Å), Met163 (2.16 Å),His164 (2.20 Å), Gln192 (2.09 Å)
10	Quercetin (5280343)	1.50 μm	−7.94	His164 (2.23 Å), Glu166 (2.12 Å), Asp187 (2.02 Å),
11	Kaempferol (5280863)	1.47 μm	−7.96	Glu166 (2.07 Å), Asp187 (2.05 Å), Thr190 (1.92 Å), Gln192 (2.026 Å)
12	Pinocembrin (68071)	2.07 μm	−7.76	Thr190(1.73 Å)
13	Rutin (6728944)	1.11 μm	−8.12	Asn142 (2.10 Å), His164(2.15 Å), Glu166 (1.93 Å), Thr190 (2.06 Å),Gln189-NH (1.79 Å), Gln189 C = O (2.19 Å), Gln189 C = O (2.08 Å).
14	A1	15.17 nm	−10.67	His163 (2.13 Å), Glu166 (2.03 Å)(π-π interactions with Phe140 and Tyr54)
15	A2	5.94 nm	−11.22	Glu166 (2.14 Å), Ser144, OH,(1.89 Å), Ser144, NH (1.99 Å)
16	A4	6.66 nm	−11.15	Cys145 (2.16 Å), His164 (2.21 Å),(π-π interaction with His172)
17	A5	30.51 nm	−10.25	His163 (1.73 Å), Thr190 (1.87 Å),Glu192 (1.94 Å)
18	A7	13.08 nm	−10.75	Glu166 (2.08 Å), His164 (2.13 Å),His163 (2.01 Å), (π-π interaction with Phe140)
19	A8	9.10 nm	−10.97	His163 (1.66 Å), Glu166 (2.02 Å),Glu192 (2.8 Å), Thr190 (1.87 Å)
20	A9	22.20 nm	−10.44	Glu166 (2.09 Å), Gln192 (1.85 Å)
21	A11	22.64 nm	−10.43	Gly13 (2.24 Å), Thr190 (1.86 Å),Glu192 (2.14 Å)
22	A12	10.84 nm	−10.87	Gln189 (2.27 Å), Glu166(2.15 Å), Glu192 (1.68 Å)
23	A20	11.33 nm	−10.84	Gly143 (1.74 Å)
24	A38	46.88 nm	−10.00	Glu166 (1.99 Å), His163 (1.93 Å)
25	A40	23.39 nm	−10.41	Glu166 (2.12 Å)
26	A68	11.67 nm	−10.82	Glu166 (2.23 Å)

Our obtained results were also compared with the recently published inhibitors for possible interaction with 6LU7-the main protease (Mpro) of the SARS-CoV-2. The docking analyses were also performed for the native drug N3 (−8.88 kcal/mol) and ebselen (−6.45 kcal/mol) as given in Table 1. Furthermore, we compared these to the FDA approved drugs Remdesivir (−8.42 kcal/mol), Chloroquine (−7.01 kcal/mol); Hydroxychloroquine (−7.39 kcal/mol); Niclossamide (−7.60 kcal/mol) and EIDD-2801 (−6.93 kcal/mol). The natural flavonoids Quercetin (−7.94 kcal/mol), Kaempferol (−7.96 kcal/mol), Pinocembrin (−7.76 kcal/mol) and Rutin (−8.12 kcal/mol) also showed very similar scores like the FDA approved drugs (Table 1). Ursodeoxychloric acid (−9.42 kcal/mol) showed a higher docking score than all the other natural products and the FDA approved drugs featured (Table 1). Our calculations (−8.42 kcal/mol) with AutoDock4.2 showed comparable interactions with the main protease of the SARS-CoV-2 with the Remdesivir drug by Ji et al., using Schrodinger docking suites (−7.215 kcal/mol) (Hall Jr & Ji, 2020).

Orientation and binding interaction of A2 and A4

The two best inhibitors found from our docking studies were A2 (−11.22 kcal/mol) and A4 (−11.15 kcal/mol) with calculated Ki of 5.94 nM and 6.66 nM, respectively, see structures in Table S1. The oxodiazole moiety of A2, oriented toward the S1 site and establishes a hydrogen bond with residue Glu166, and weakens the salt bridge between positively charged N-terminal Ser1 (one monomer) of the negatively charged Glu166 (another monomer), thus destablilizes the Mpro pocket (Huynh et al., 2020) this may be the reason for increasing the potency of compound A2. The Ebselen a Mpro inhibitor FDA approved drugs Remdesivir, Chloroquine and Niclosamide were also formed a similar type of hydrogen bonding with Glu166. The the Phenol segment of A2 has formed two H-bond with Ser144 by fully occupying S1' subsite. The chlorobenzene segment orients toward the S4 pocket, while the isatin segment is oriented towards the bulky hydrophobic S2 site (Figure 1). Compound A4 showed a different picture with its benzo[d]oxazol-2(3H)-1 segment, the benzyl moiety orienting towards the S2’ subsite while the oxazole moiety has accommodated well in S1’subsite by forming a strong hydrogen bond between the ring oxygen with the catalytic Cys145 residue; The chlorobenzene segment near oxazole moiety orients toward the S1 site and form a π-π interaction with the residue His172. The connecting 2° amine segment of the A4 is oriented towards the bulky hydrophobic S1 site and forms a hydrogen bond with the residue His164. The chlorobenzene segment near oxadiazole moiety orients towards the S4 pocket (Figure 1). In Figure 1, the other higher docked compounds like A2, A4, A8 and A20 are also represented. The best docked compounds A2 and A4 have been selected for further analysis with MD simulations using GROMACS.

Figure 1.

The binding mode of the four best docked compounds in the active site of the SARS-CoV-2 virus Mpro (PDB ID: 6LU7). The interacted amino acid residues and the distances in Å are given in yellow. Top, the ligands are shown together in the binding pocket. At bottom, 4 panels, individual compounds docked into the binding site of SARS CoV-2 virus Mpro. A2: Saffron; A4: Sky; A8: Ocean blue; A20: Dark brown.

Structural activity relationship

In this article, we performed screening of 16 distinct heterocyclic skeletons (Table 2) with a total of 118 compounds, 5 natural products and 8 repurposed drugs (including Ebselen, which has a similar structure to isatin (Jin et al., 2020b)) (Figure 2). The docking analysis of several isatin compounds A2, A4, A8, A20, showed the best interaction with the Mpro (above −10.84 kcal/mol) as presented in Table 1 and Figure 1. In addition, highly potent compounds with Ki under 10 nM, A2, A4 and A8 were comprised of three separate ring systems, branched from a central carbon enabling multiple binding interactions with Mpro.

Table 2.

Structural activity relationship of designed compounds.

1	A1 to A3		A2 (−11.22)
2	A4 to A6		A4 (−11.15)
3	A7 to A12		A8 (−10.97)
4	A15 to A29		A20 (−10.84)
5	A30 to A37 and A38 to A40		A39 (−10.41)
6	A42 to A45 and A88 to A92		A91 (−9.60)
7	A94 to A101		A95 (−9.58)
8	A75 to A87		A86 (−9.47)
9	A57 to A63		A58 (−9.44)
10	A102		A102 (−9.44)
11			Ursodeoxycholic acid (−9.42)
12	A50 to A56		A51 (−9.23)
13	A103 to A105		A105 (−8.96)
14	A106 to A112		A107 (−8.94)
15	A47		A47 (−8.07)
16	A48 and A49		A48 (−7.98)
17	A113 to A117		A117 (−7.30)

Figure 2.

IC50 and binding energies of potent compounds similar to Ebselen - a potent published SARS-CoV-2 Mpro inhibitor (Jin et al., 2020b). The top 3 were identified for SARS-CoV1 (Chen et al., 2005; Zhou et al., 2006).

According to a literature survey, isatin based compounds are well known for their anticancer (Bhatt et al., 2020; Popp, 1969) and antibacterial activities (Pandeya et al., 2000; Varma & Nobles, 1967). They have also been known for their antiviral activities for a various pathogen virus (Bauer & Sadler, 1960; Burger, 1979)89,90 including HIV (Pauwels et al., 1988; Teitz et al., 1994) and SARS-CoV-1 (Selvam et al., 2008). Oxadiazole, based derivatives have also been found very important for potent biological functions. Especially antiviral, anti-HIV, anti-HCV, anti-HBV, anti-HSV activities, etc (Dong et al., 2016; Li et al., 2011). Triazole derivatives also occupy a pivotal position in modern medicinal chemistry and several derivatives are well known for their applications in medicine especially anticancer and antiviral drugs (Kharb et al., 2011; Xu et al., 2019). The specific reason behind the isatin; oxadiazole and triazole scaffold acting as potential anticancer, antiviral drugs are their strong tendency to accept multiple hydrogen bonds from their target protein. The Oxidiazole based compounds A1, A2, A4, A5, A7-A9, A11 & A12 and triazole based compound A20 showed good binding energies and Inhibition constants. Previously designed Cyclin-dependent kinases inhibitors known for their anticancer activity, A38 and A40 also showed high binding interaction with the Mpro of the SARS CoV-2, though their inhibition constants were higher than the isatin oxidiazoles/triazoles based compounds (Table 1). Compound A68, also showed high scores in our calculation (−10.82 kcal/mol), also comprised of a three separate ring system. A68 inhibitor was already reported by Fischer et al.10. Other compounds with very efficient docked scores that were subjected to further analysis, from our docking calculations are A1 (−10.67 kcal/mol); A5 (−10.25 kcal/mol); A7 (−10.75 kcal/mol); A8 (−10.97 kcal/mol); A9 (−10.44 kcal/mol); A11 (−10.43 kcal/mol); A12 (−10.87 kcal/mol); A20 (−10.84 kcal/mol); A38 (−10.00 kcal/mol) and A40 (−10.41 kcal/mol). All other compounds with lower docked scores are given in the supporting information (Table S2). A2 (−11.22 kcal/mol) and A4 (−11.15 kcal/mol) compounds showed the highest docked scores and were subjected to further analysis.

Chelpg interaction

Molecular recognition between a protein and small molecules (ligand) occurs at their surfaces. The extent of binding forces in the protein-ligand complex largely depends on the electronic distribution of the ligand which can be predicted by the electrostatic charges of the ligand molecule. The electrostatic charges of the electronegative atoms of the A2 and A4 ligands are given in Figure 3. Such molecular details of the ligands helps predicting the binding probabilities of the ligands in the S1, S2, S1’, S2’ and S4 pockets as discussed above. Our charge analysis shown that more electronegative atoms are present in A2 predicting it to interact better with the protein than A4. The ChelPG charges of the electronegative atoms for the other potential ligands are given in the supporting information (Figure S1). As observed from the charge analysis, A1-an analogue of A2, showed less interaction to Mpro due to the presence of chlorine atom in place of the hydroxide. A5, an analogue of A4, showed low interactions despite of the presence of hydroxide ion. This may be due to the steric hindrance between the isatin and the chlorobenzene occurring, owing to the flexibility of bonds. A7, A8, A9, A11 and A12 having the similar structures showed very comparable interactions. The ChelPG charge analysis could not be performed for A7, A8 and A9 ligands due to restrictions in the Gaussian format with presence of bromine atom. The ligand A20, having a distinct structure with a hydrophobic end, showed higher interaction to the protein. The ChelPG charges of the other important ligands (A38, A40 and A68) with good interactions to the Mpro are also given in the Figure S1 in the supporting information.

Figure 3.

The ChelpG charge of the electronegative atoms for both the compounds A2 and A4 are given here. The circle in coloured lines and the Sj (j = 1,1',2,2',4) term denotes the pockets where these sections of the compounds may interact in the Mpro protein. Carbon: Ash; Oxygen: Red; Chlorine: Green; Nitrogen: Blue; Hydrogen: White.

Molecular electrostatic potential analysis

Molecular electrostatic potential calculations (MESP) were performed for A2 and A4 including the shape and surface of the compound since electronic distribution in the compound is one of the most important factor to predict the binding to a protein. In the MESP surface below, the electronegative atoms are shown in red and they act as regions for H-bond acceptor, while the electron poor atoms are designated in blue, which act as H-bond donors. The green colour portions are neutral in nature and are regions where π- and other types of π-staking interactions are important (Figure 4). The MESP analysis shows A2 can form at least three hydrogen bonds, while A4 can form two (red patches). The electron density of the A2 was found to be higher than the A4 (green to yellow). The incorporation of an OH group in A2, rather than a Cl as in A4, is the reason behind the increasing electron density of A2, which should increases its binding interactions. Such enhancement of the electron density is favourable for the π-staking interaction. The Cl atom for both compounds also occupies a large electron density surface, showed in green, with negative values which is also suitable for multiple π-staking interactions. The small red, yellow and blue patches on the large green surface of compounds are balancing the hydrophilic and hydrophobic parts which are essential for good binding to the protein. The MESP analysis of the other ligands as given in Table 1 are also calculated and are given in the Figure S2, in the supporting information. The red, yellow and blue patches on the large green surface of compounds in these ligands are observed and the interactions patterns are clear from the MESP.

Figure 4.

The calculated Molecular Electrostatic Potential surface for the A2 and A4 compounds showing the electron density surface. The red patches hydrogen bonding donor and blue hydrogen bond acceptors sites. π-interaction sites are given as green color and yellow patches showed higher electron density than the green sites. The surfaces shown correspond to an isosurface value of 0.03 electrons/a.u.3.

Molecular dynamics simulation

Further calculations using MD simulation tools were performed for A2 and A4 to check the interaction analysis and binding ability with the Mpro protein of the SARS-CoV-2 to evaluate the results of the docking and the DFT charge analyses.

System and ligand stability inside the Mpro active site

The stability of each simulated model was determined by the calculations of backbone atom RMSD for ligand − Mpro complex as shown in Figure 5. The RMSD of A2 and A4 with- Mpro rises till 20 ns and stabilizes after 10 ns at 2 Å. This clearly indicates that ligands are stabilized inside the pocket and don’t change their orientation in the active site of Mpro. Similarly protein backbone atoms for A2 showed stability around 2.5 Å. However, A4 showed higher RMSD and maintained at 3 Å, which directly resembels to the charge analyses. Furthermore, A2 and A4 hydrogen bond with Mpro are plotted along the 50-ns MD simulation, shown in Figure 6. These simulations indicated that on an average A2 and A4 show approximately two hydrogen bonds and are in accordance with our docking results and further suggested the stability of the ligand in active site.

Figure 5.

The root mean square deviation (RMSD) plots for the A2 and A4 compounds interacting with the Mpro protein for the entire 50 ns MD simulation. Top, RMSD of A2 and A4, and bottom of the RMSD of the backbone atoms of Mpro during binding of the A2 and A4. The green line represents A2 and blue line represents A4.

Figure 6.

Number of Hydrogen bond of A2 and A4 with Mpro plotted along the 50-ns MD simulation.

Other than the Hydrogen bond, hydrophobic interactions were established among non-polar amino acids and provide the stability of proteins in solution by shielding the non-polar amino acids in hydrophobic cores, away from the aqueous environment. Solvation plays an important role in ligand-protein association and has a strong impact on comparisons of binding energies for dissimilar molecules.

Change in Solvent Accessible Surface Area (SASA) directly indicates the change of hydrophobic core caused by changes in tertiary structure. The calculation of the SASA of A2 and A4 do not show any drastic change that indicate the ligand were stable during the simulation run, as shown in Figure 7, and showed very similar hydrophobic area. This indicates that A2 and A4 are less exposed to water molecules and buried in Mpro active site.

Figure 7.

Hydrophobic, Hydrophilic, and solvent accessible surface area for A2 and A4 with Mpro plotted along the 50-ns MD simulation. Blue colour showed ΔGsolv, the free energy of ligand desolvation, based on the solvent accessible surface plotted along the 50-ns. Red line showed the hydrophilic SASA; black line, the hydrophobic SASA and the green line showed the combination of hydrophilic and hydrophobic SASA.

Binding energy analysis

The calculated binding energy using the mmpbsa approach indicates that A2 has higher energy (−48.33 kcal/mol) in comparison to A4 (-33.10 kcal/mol), Table 3, similar to the finding in the docking analysis, Table 1. Individual components of energy for each system indicate that A2 and A4 both show higher van der Waals energies and very comparable SASA energies and polar salvation energies. However a significant difference was observed in the electrostatic energies. As explained in the ChelpG charge analysis (Figure 3), the electrostatic interactions for the A2 showed much higher interactions than the A4 ligand, which effects the binding energy difference of the two systems. Thus, the binding energies are governed by the van der Waals energy and electrostatic energy.

Table 3.

MM/PBSA binding free energy of the selected compounds compared with the known inhibitors of the SARSCoV-2 Mpro.

Complex	Van der Waal energykcal/mol	Electrostatic energykcal/mol	Polar solvation energykcal/mol	SASA energykcal/mol	Bindingenergykcal/mol
A2	−58.38+/− 0.26	−10.08+/−0.23	24.74+/−0.15	−4.61+/−0.02	−48.33+/−0.36
A4	−49.89+/− 0.27	−2.43+/−0.19	23.61+/−0.24	−4.39+/−0.02	−33.10+/−0.20

Energy contribution of each residue to binding

To identify important interacting residues with A2 and A4 the contribution of each residues to the binding energy was calculated and plotted (Figure 8) indicating three important regions: amino acids 26 to 56, amino acids 136 to 148, and amino acids 162 to 194. These regions showed significant energy contribution to the binding energy. Leu27 Asn51 and Pro52 showed negative energy in both systems while His41 showed positive energy in A2. In the second region (amino acid 136 to 148) A4 showed higher energy values for residues Leu141, Asn142, Gly143, Ser144, Cys145. In the third region (amino acid residue 162 to 194) A2 had higher binding energy values for interactions with residues Met165, Leu167, Pro168, Phe185, Asp187, Thr190, Ala191 and Gln192. While A2 interacted with amino acid from region 162 to 192 and residue 44, A4 showed good interaction with region 140 to 146.

Figure 8.

The contribution energy of each residue to Mpro. Top box is a general energy plot of all residues. Bottom three boxes: the Y axis presents Energy [kcal/mol] and the X axis indicated the energy contribution of each residue by number in Mpro, green triangles represent A2 and blue triangle A4. Clear perturbations are detected at amino acids 26 to 56, amino acids 136 to 148 and amino acid 162 to 194. The respective structures of the ligand-protein complexes are given above.

The Mpro complex interacts with the ligands A2 and A4 with the different amino acid residues between energy −2.4 kcal/mol to −0.24 kcal/mol. A2 formed more interaction with high binding energy (up to 2.4 kcal/mol) in the active site, Figure 9. Cys44, Met165, Gln189 and Ala191 are major residues that highly contributing to binding with Mpro, this indicates the ligand fits well into the active site of Mpro. However, residues participating in A4 binding are somewhat different, very few residues bind with more than −0.95 kcal/mol, the major binding residue were found to be Gln189, Met165, Cys145 and Pro168. Several amino acids were common to binding both A2 and A4: Leu27, Cys44, Met49, His164, Met165, Leu167, Asp187, Gln189, Thr190, Ala191 and Gln192 these also contribute to the binding energy. Overall this analysis supports the better binding interactions of the A2 on the active site of Mpro than A4, Figure 10.

Figure 9.

Major residues contributing to the binding of A2 and A4 compounds are plotted. The Y axis represents the energy in kcal/mol for the energy contribution and X axis showed amino acid position for each residue.

Figure 10.

The different interactions the compounds A2 and A4 at the active site of the Mpro protein of the SARS-CoV-2 calculated using the Molecular dynamics simulations. Carbon: Ash; Oxygen: Red; Chlorine: Green; Nitrogen: Blue; Hydrogen: White. The different types of interactions are also given in the figure for both the compounds.

Flexible region and catalytic dyad of Mpro

The RMSF captures, for each residue, the fluctuation from its average position, gives insight into the flexibility of regions of the Mpro. The fluctuations per residue of COVID-19 Mpro indicated that after the binding, three regions, amino acid residues from 40 to 54, 124 to 144 and 180 to 200, and in addition 272 to 306 showed higher fluctuations, Figure 11. Moreover, these residues are also participating in ligand binding. These might suggest that these flexible regions are important in substrate or ligand binding.

Figure 11.

The RMSF plot for each residue out of the four highly fluctuating regions, is show over the 50 ns MD simulation. A2 (green) and A4 (blue); (A): showing fluctuations in region 40 to 54; (B): Showing fluctuation in region 124 to 144 & 180 to 200; (C): showing fluctuation in region 272 to 306 during Mpro-Ligands complex.

ADME parameters prediction

In order to screen the pharmacokinetic behaviour, the in-silico ADME parameters of the designed ligands with free energy of binding higher than −10 kcal/mol were calculated using the SwissADME software (http://www.swissadme.ch). Lipinski’s rule of five is a prerequisite to ensure drug-like properties when using rational drug design. The designed derivatives follow the four properties of Lipinski’s rule of five (mol. wt. ≤ 500 Da; log P o/w ≤ 5; HBD ≤ 5; HBA ≤10). All the ligands have showed ADME parameters in an acceptable range and possess significant drug-like characteristics based on Lipinski’s rule of 5, Table 4.

Table 4.

In-silico predicted ADME properties of the designed active compounds.

Compound	Mol. Wt.	Number of HBD	Number of HBA	MR	Log P o/w	Log S(ESOL)	Solubility(mg/ml)	Lipinski Rule
A1	450.27	0	5	119.1	4.53	−5.98	4.66e-04	Yes; 0 violation
A2	431.83	1	6	116.02	3.52	−5.25	2.41e-03	Yes; 0 violation
A4	453.28	1	5	118.36	4.80	−6.54	1.31e-04	Yes; 0 violation
A5	434.9	2	6	115.37	3.66	−5.80	6.89e-04	Yes; 0 violation
A7	544.18	1	5	131.38	4.16	−7.17	3.70e-05	Yes; 1 violation: MW > 500
A8	525.74	2	6	128.40	3.17	−6.43	1.96e-04	Yes; 1 violation: MW > 500
A9	554.74	1	7	135.20	3.92	−6.63	1.29e-04	Yes; 1 violation: MW > 500
A11	446.84	2	6	120.70	2.86	−5.52	1.35e-03	Yes; 0 violation
A12	475.84	1	7	127.50	2.26	−5.72	8.99e-04	Yes; 0 violation
A20	432.47	0	4	131.48	3.18	−5.78	7.14e-04	Yes; 0 violation
A38	341.41	1	2	102.48	3.18	−4.52	1.04e-02	Yes; 0 violation
A40	427.31	1	4	111.50	3.89	−5.58	1.12e-03	Yes; 0 violation

Mol. Wt.: Molecular weight, HBD: Hydrogen bond donor, HBA: Hydrogen bond acceptor, MR: Molar Refractivity. Log Po/w: Octanol/water partition coefficient, Log S: Aqueous solubility.

Conclusion

In conclusion, we have studied the binding of a series of heterocyclic moieties to the active site of Mpro of the SARS-CoV2. Molecular docking analysis show that the isatin based derivatives possess excessive interactions to Mpro, with binding energies of over −10 kcal/mol. The best docked ligands A2 and A4 consisted of an Isatin moiety along with an oxidiazole ring, and both moieties are highly used in anticancer, antiviral, antibacterial drug development. Further, the binding interactions were more vividly discussed with a comprehensive MD simulation analysis using GROMACS, and binding energies were calculated using the MMPBSA method on the active site of the protein. These analyses show hydrogen bonds and π-π interactions with the residues His41, Cys145, Val42, Cys44, Thr45, Ser46, Glu47, Asp48, Met49 and Glu166, Met165, Leu167 and Pro168. Our analysis also shows that hydrophobic and electrostatic interactions are most important to predict the binding energies of the ligands to the protein and that A2 interacts more strongly than the A4. This finding was reinforced by the electrostatic charge analysis and MESP analysis through DFT calculations. All the ligands have showed ADME parameters in an acceptable range and possess significant drug-like characteristics based on Lipinski’s rule of 5. Our findings indicate that both isatin and oxidiazoles based drugs are indeed very good candidates as potential inhibitor for Mpro of SARS-CoV-2 and A2 and A4 compounds are the best candidates with extremely high potency in silico and further in vitro analysis will be performed as future studies.

Computational section

Quantum chemical calculations

All the ligands were initially optimized using Density functional theory with Gaussian 16 suite of programs to obtain the most stable structure of each system (Ogliaro et al., 2016) except the repurposed drugs with crystal structures available. All the compounds were optimized at B3LYP/6-31 + G(d) (Becke, 1993; Lee et al., 1988) level of theory in gas phase. Positive vibrational frequencies confirmed the structures to be minimum. PDB structures were saved using Gaussview and used for the Molecular docking studies. Electrostatic charges were calculated for the specific compounds using the ChelpG charge analysis (Martin & Zipse, 2005; Rozas, 1997; Tsuzuki et al., 1996) implemented in Gaussian software. The Molecular Electrostatic Potential (MESP) was generated for both the ligands at the same level of theory.

Molecular docking studies

The crystal structure of SARS-CoV-2 has been retrieved from the rcsb.org (PDB ID: 6LU7) (Jin et al., 2020b) and used to generate initial 3D coordinates. Co-crystallized water molecules were deleted along with addition of polar Hydrogen and compute gasteiger charges. The structures of selected inhibitors were superimposed against pre-docked ligand in the PDB, and the latter was then removed to generate initial conformation of inhibitor at the active site of SARS-CoV-2. Grid box was then determined by the native ligand position on the binding site (Thr26, Tyr54, Phe140, Asn142, Gly143, Ser144, Cys145, His163, His164, Glu166, and His172). Finally both the Auto-grid and Auto-dock were run with the default parameters (Badavath et al., 2016a). All molecular docking analysis were performed by AutoDock4.2 using earlier reported protocol (Badavath et al., 2016b; Badavath et al., 2015; Badavath et al., 2017; Badavath et al., 2016c; Gangireddy et al., 2019; Munusamy et al., 2019; Nath et al., 2018; Nayak et al., 2015; Vishnu Nayak et al., 2013) software version 4.2 (Forli et al., 2012). After completing the molecular docking process, the generated DLG (.dlg) file was analyzed using AutoDock v4.2 tools for all the possible interaction and binding free energies following the methodology (Huey et al., 2007). Top scoring molecules from largest cluster were considered for the interactions.

Gromacs 4.6.2 (Hess et al., 2008; Van Der Spoel et al., 2005) with GROMOS96 54a7 force field (Schmid et al., 2011) was used for MD simulation studies of two system 50 ns each. The PRODRG2 Server (Schuttelkopf & Van Aalten, 2004) was used to generate the topology of A2 and A4. Each system was placed in the centre of cubic box having distance of 1.0 nm between protein and edge of the simulation box and solvated with explicit water molecules, with the SPC solvation model (Mark & Nilsson, 2001). Before proceeding towards energy minimization, all the systems were neutralized by adding Na+/Cl− ions accordingly. Steepest descent method was used for energy minimization of each system. MD simulations with NVT (isochoric-isothermal) and NPT (isobaric-isothermal) ensembles (N = constant particle number, V = Volume, P = Pressure, T = Temperature) were performed for 1 ns to equilibrate the protein-ligand system for constant volume, pressure (1 atm) and temperature (300 K). To calculate electrostatic interaction, Particle Mesh Ewald (PME) algorithm (Darden et al., 1993; Essmann et al., 1995) was used with grid spacing of 1.6 Å and a cutoff of 10 Å and LINCS method (Hess et al., 1997) was used to restrain the bond length. Finally, 50 ns of production MD was performed for each system and trajectories were saved at every 2 ps.

Binding energy calculation for A2 and A4

G_mmpbsa (Kumari et al., 2014) was developed using two widely used open source software i.e. GROMACS and APBS and it has similar user interface like other GROMACS tools. The tool calculates components of binding energy using MM-PBSA method except the entropic term and energetic contribution of each residue to the binding using energy decomposition scheme. Binding free energy estimation of all the docked complexes was done by MM/PBSA method through g_mmpbsa module. For each simulated system, to calculate average binding energy by using bootstrap analysis, a total of 200 snapshots were taken from the last 10 ns of the trajectory on the intervals of 50 ps. The module calculates electrostatic interactions, van der Waals interactions, polar solvation energy and non-polar solvation energy.

ADME calculation

ADME parameters of the designed ligands were calculated using the SwissADME (http://www.swissadme.ch) (Daina et al., 2017) accessed on- 5th May 2020. All the compounds were drawn in Chem Draw software and smiles were taken as input, before being used in SwissADME. Lipinski’s rule of five (Lipinski et al., 1997).

Click here for additional data file.