Potential of NO donor furoxan as SARS-CoV-2 main protease (Mpro) inhibitors: in silico analysis.

The sharp spurt in positive cases of novel coronavirus-19 (SARS-CoV-2) worldwide has created a big threat to human. In view to expedite new drug leads for COVID-19, Main Proteases (Mpro) of novel Coronavirus (SARS-CoV-2) has emerged as a crucial target for this virus. Nitric oxide (NO) inhibits the replication cycle of SARS-CoV. Inhalation of nitric oxide is used in the treatment of severe acute respiratory syndrome. Herein, we evaluated the phenyl furoxan, a well-known exogenous NO donor to identify the possible potent inhibitors through in silico studies such as molecular docking as per target analysis for candidates bound to substrate binding pocket of SARS-COV-2 Mpro. Molecular dynamics (MD) simulations of most stable docked complexes (Mpro-22 and Mpro-26) helped to confirm the notable conformational stability of these docked complexes under dynamic state. Furthermore, Molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) calculations revealed energetic contributions of key residues of Mpro in binding with potent furoxan derivatives 22, 26. In the present study to validate the molecular docking, MD simulation and MM-PBSA results, crystal structure of Mpro bound to experimentally known inhibitor X77 was used as control and the obtained results are presented herein. We envisaged that spiro-isoquinolino-piperidine-furoxan moieties can be used as effective ligand for SARS-CoV-2 Mpro inhibition due to the presence of key isoquinolino-piperidine skeleton with additional NO effect. Communicated by Ramaswamy H. Sarma.

Introduction

Novel Coronavirus (SARS-CoV-2) causing COVID-19 disease was originated from Hubei province of China, and has been exponentially increasing in whole world with United States of America, China, Europe, India and Middle East countries becoming the greatest concern. It is one of three Coronaviruses which have been associated with pneumonia, including severe acute respiratory syndrome coronavirus (SARS-CoV) and Middle-East respiratory syndrome coronavirus (MERS-CoV). The World Health Organization announced that the outbreak of this deadly and fast-spreading infectious pneumonia constitutes a global health emergency (WHO, 2019-nCoV Situation Report, 2020). As on 22nd June 2020, the total numbers of cases around the world were recorded to be 9,072,803 with more than 471,178 deaths (worldometers.info/coronavirus).

WHO has been declared COVID-19 as a global pandemic disease but currently there is no effective treatment available and we only depend on preventive and supportive therapies (Cheng, 2019). Researchers around the globe are desperately working to discover either effective vaccines or antiviral drugs, however these experimental methods for drug discovery are expensive and time-consuming. In silico based screening methods offer easy way, fast and low-cost techniques for novel testable hypotheses to discover potential lead drug candidate (Cheng, 2019). There are many clinical trials in progress by pharmaceutical companies to afford the potential agents but very limited data with respect to in vitro and in vivo activities of the lead drugs is available in literatures which are in clinical studies for treatment of COVID-19 (Lu, 2020).

Nitric oxide (NO) is involved in a broad range of processes and act as an important signaling molecule. Thus, the designing of NO-donors/hybrid molecules with NO releasing group may be beneficial for a variety of disorders, such as: cardiovascular, inflammatory, bacterial, fungal, parasitic, ocular diseases and cancer, arthritis, asthma, cerebral ischemia, Parkinson’s disease, neurodegenerative diseases, seizures and viral (Serafim et al., 2012). The well know NO donor, furoxan is 1,2,5-oxadiazole-2-oxide that could release NO in the presence of thiols as well as by enzymatic action in tissue with complex process that probably involves generation more than one redox form of NO (Gasco et al., 2004).

Nitric oxide (NO) is an antiviral effector of the immune system and inhibits the replication of a variety of viruses (Saura et al., 1999). NO can inhibit the replication of viruses that encode cysteine proteases, such as members of the Picornavirus family and the Coronavirus family (Mannick, 1995; Reiss & Komatsu, 1998). Endogenously generated NO can inhibit the replication cycle of SARS-CoV and hence exhibit the desired antiviral effect (Åkerström et al., 2009). NO inhibits the replication of SARS-CoV by dual mechanisms. First, NO or its derivatives cause a reduction in the palmitoylation of nascently expressed spike (S) protein which affects the fusion between the S protein and its cognate receptor i.e. angiotensin converting enzyme 2 (ACE-2). Second, NO or its derivatives cause a reduction in viral RNA production in the early steps of viral replication, and this could possibly be due to an effect on one or both of the cysteine proteases encoded in Orf1a of SARS-CoV (Åkerström et al., 2009).

In-vitro studies have identified significant inhibition activity of NO against SARS-coronavirus infection by S-Nitroso-N-acetylpenicillamine, a nitric oxide donor (Keyaerts et al., 2004). In patients with severe acute respiratory syndrome (SARS), inhalation of nitric oxide (NO) has improved arterial oxygenation and enabled the reduction of inspired oxygen therapy and airway pressure (Chen et al., 2004). Currently, clinical trial testing of inducible nitric oxide (iNO) is underway in COVID-19 infected patients complicated with Acute respiratory distress syndrome (ARDS) [NCT04306393 and NCT04305457] (clinicaltrials.gov).

Theophylline derivative and pyrimidone derivative were identified as inhibitor of RNA binding agent against 2019-nCoV N protein (N terminal domain) by in-silico studies (Sarma et al., 2020). Many natural compounds were studied for their properties to block the activity of the angiotensin converting enzyme 2 (ACE-2) as a receptor for SARS-CoV-2 and as a potential therapeutic target of the COVID-19 virus using in silico study (Abdelli et al., 2020). The extensive research is required by using computer simulations such as (A) Exploration of the energetic binding affinity of SARS-CoV-2 Mpro with new inhibitors based on free energy calculations. (B) Investigation of the structural properties, flexibility, conformational changes of the SARS-CoV-2 Mpro, and study interaction pattern between virus and membrane, virus and inhibitor. C) Monitoring of the thermodynamics (virus-water interaction) properties of SARS-CoV-2 Mpro in the presence/absence of the antiviral inhibitor (Boopathi et al., 2020). Thus computational study could be an effective and swift approach on discovering novel uses for existing medications or approved antiviral drugs to inhibit SARS-CoV-2 Mpro (Al-Khafaji et al., 2020; Khan et al., 2020a, 2020b; Lobo-Galo et al., 2020; Muralidharan et al., 2020; Pant et al., 2020) or SARS-CoV-2 RdRp (Elfiky, 2020a). The researchers have also studied the potential application of natural origin phytochemicals and bioactive derivatives with an anti-viral properties as inhibitors of SARS-CoV-2 Mpro using molecular docking (Aanouz et al., 2020; Das et al., 2020; Elfiky, 2020b; Gyebi et al., 2020; Islam et al., 2020; Joshi et al., 2020; Sinha et al., 2020). The advent of high resolution experimental structure of the target, the predictions of virtual screening studies and binding energy calculations are generally more accurate. The high resolution experimental structure of the target, the main protease of SARS-CoV-2 (Jin et al., 2020; Zhang et al., 2020), was utilized in the current study as the target for molecular docking study.

Keeping in view the importance of early screening of NO donor hybrids for effective anti-viral attributes, the present study focus on evaluation of phenyl furoxan as potential drug candidates by using molecular docking, MD simulation and MM-PBSA approach. The predictions of this study will provide information that can be utilized further for development of efficient drugs by testing it in vitro, in vivo and through clinical trials to control SARS-CoV-2.

Experimental

Molecular docking

Preparation and selection of target

The crystal structure of the SARS-CoV-2 Mpro bound to potent broad-spectrum non-covalent inhibitor X77 (PDB ID: 6W63) was selected for the docking analysis (www.rcsb.org; www.pdb.org). The protein structure was selected on the basis of resolution and method of determination. The selected protein structure was analyzed for the Ramachandran plot using V life MDS 4.6 (www.vlifesciences.com). Selected protein structure of SARS-CoV-2 Mpro was refined for removal of native ligand structure X77 keeping its native conformation intact (Pettersen et al., 2004). The refined protein structure was utilized further for docking analysis.

Preparation of ligand structures

Molecular structures of the selected furoxan derivatives was developed using 2D molecular builder in V life MDS and further converted in to 3D geometry and optimized via application of merck molecular force field (MMFF) to get optimized ligand structures (vlifesciences.com). These optimized ligand structures were utilized for docking analysis.

Redocking analysis

The redocking protocol is utilized for the validation of the docking protocol. The redocking was carried out via docking of co crystallized ligand X77 in SARS-CoV-2 Mpro. The RMSD was determined between docked ligand X77 and co crystallized ligand X77. The RMSD value less than 1.0 Å was considered to be satisfactory and which will ascertain that the binding mode of the selected furoxan derivative will be similar to that of co crystallized ligand X77.

Docking analysis

Docking simulations was performed to predict the possible binding mode of the selected furoxan derivatives. The Grip based docking analysis followed by redocking using AutoDock Vina (Trott & Olson, 2010) which uses Lamarckian Genetic Algorithm (LGA) was performed keeping ligand structures (Table 1) in the flexible conformations. The rotational angle between the two conformations was kept at 10° and total number of rotation to 30 in grip based docking. The best docking pose was selected on the basis of the docking score, binding energy and types of interactions.

Table 1.

Docking result of selected NO donor furoxan.

Sr. No	Structure & Name	Binding Affinity (kcal/mol)	Docking Score	Interactions
				H bond	Aromatic	Charge	Hydrophobic
1	3-({4-[bis(4-fluorophenyl)methyl] piperazin-1- yl}methyl)-4-(4-chloro phenyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.1	−58.33		HIS41	GLU166	LEU141 ASN142 MET165 GLU166
2	3-({4-[bis(4-fluorophenyl)methyl] piperazin-1- yl}methyl)-4-phenyl-1,2,5-oxadiazol-2-ium-2-olate	−6.7	−66.94	SER144	HIS41	HIS41	ASN142 MET165
3	3-({4-[bis(4-fluorophenyl)methyl] piperazin-1-yl}methyl)-4- (4-methoxyphenyl)-1,2,5-oxadiazol-2-ium-2-olate	−8.1	−67.36		HIS41	GLU166	ASN142 MET165 GLU166
4	3-({4-[bis(4-fluorophenyl)methyl] piperazin-1-yl}methyl)-4- (2-nitro phenyl)-1,2,5-oxadiazol-2-ium-2-olate	−8.4	−46.14	GLN189			MET49 GLN189
5	3-({4-[bis(4-fluorophenyl)methyl] piperazin-1-yl}methyl)-4- (2-methyl phenyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.6	−67.09	ASN142	HIS41	GLU166	ASN142 MET165 GLU166
6	3-({4-[bis(4-fluorophenyl)methyl] piperazin-1-yl}methyl)-4- (4-nitro phenyl)-1,2,5-oxadiazol-2-ium-2-olate amine	−7.9	−63.99		HIS41		MET49 ASN142 GLU166
7	3-({4-[bis(4-fluorophenyl)methyl] piperazin-1-yl}methyl)-4- (4-fluoro phenyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.2	−70.55	SER144	HIS41 HIS163		MET165 GLU166 LEU167 PRO168
8	3-{[4-(diphenylmethyl)piperazin-1-yl] methyl}-4-phenyl-1,2,5-oxa diazol-2-ium-2-olate	−7.5	−76.22		HIS41	GLU166	MET165 GLU166
9	3-{[4-(diphenylmethyl)piperazin-1-yl]methyl}- 4-(4-chlorophenyl)-1,2,5 -oxadiazol-2-ium-2-olate	−8.5	−66.26	CYS145	HIS41 HIS163	GLU166	LEU141 ASN142 GLU166
10	3-{[4-(diphenylmethyl)piperazin-1-yl]methyl}- 4-(4-methoxyphenyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.1	−70.45	GLU166	HIS41	HIS41	HIS41 ASN142 HIS164 MET165 GLU166
11	3-{[4-(diphenylmethyl)piperazin-1-yl]methyl}- 4-(2-nitrophenyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.7	−64.58	GLU166 CYS145		GLU166	LEU141 ASN142 GLU166
12	3-{[4-(diphenylmethyl)piperazin-1-yl]methyl}- 4-(4-nitrophenyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.9	−69.80	CYS145	HIS163	GLU166	LEU141 ASN142 GLU166
13	3-{[4-(diphenylmethyl)piperazin-1-yl]methyl}- 4-(2-methoxyphenyl)-1,2,5-oxadiazol-2-ium-2-olate	−8.4	−47.40	HIS41	HIS41		THR25 LEU27 HIS41 VAL42 CYS44
14	3-{[4-(diphenylmethyl)piperazin-1-yl]methyl}- 4-(4-fluorophenyl)-1,2,5-oxadiazol-2-ium-2-olate	−8.4	−57.37	CYS145	HIS41 HIS163	GLU166	ASN142 MET165 GLU166
15	3-({3-oxo-2,3-dihydro-1H-spiro [isoquinoline-4,4′-piperidine]-1′-yl} methyl)-4-phenyl-1,2,5-oxadiazol-2-ium-2-olate	−8.0	−40.26	GLU166	HIS41		GLU166 PRO168
16	4-(2-nitrophenyl)-3-({3-oxo-2,3-di hydro-1H- spiro[isoquinoline-4,4′-piperidine]-1′-yl} methyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.5	−57.77	SER144	HIS41		MET165 GLU166
17	4-(2-methoxyphenyl)-3-({3-oxo-2,3-dihydro-1H- spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−7.6	−75.44	GLY143 CYS145 GLN189	HIS41		MET49 ASN142 CYC145 HIS165 ARG188 GLN189
18	4-(2-methoxyphenyl)-3-({3-oxo-2,3-dihydro-1H- spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−8.3	−62.94		HIS41		MET49 MET165 GLU166 ARG188 GLN189
19	4-(4-nitrophenyl)-3-({3-oxo-2,3-dihydro-1H-spiro [isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1,2, 5-oxadiazol-2-ium-2-olate	−8.8	−43.83	GLN192	HIS41		MET49
20	4-(4-methoxyphenyl)-3-({3-oxo-2,3-dihydro- 1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−8.4	−88.64	GLY143, SER144	HIS41		HIS41, ASN142, MET165, GLU166, PRO168. ASP187, ARG188
21	4-(4-bromophenyl)-3-({3-oxo-2,3-dihydro-1H- spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−8.5	−66.78	GLN189	HIS41		CYS44 MET149 MET165 GLU166 ARG188 GLN189
22	4-(3-nitrophenyl)-3-({3-oxo-2,3-di hydro-1H- spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−9.6	−90.93	CYS145 SER144	HIS41 HIS163		HIS41, CYS44, MET49, ASN142, GLY143, CYS145 MET165 GLU166
23	4-(4-flurophenyl)-3-({3-oxo-2,3-dihydro- 1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−8.3	−89.44	GLN189, GLY143.	HIS41 HIS163		MET49, ASN142, MET165, ARG188 GLN189
24	4-(2-chlorophenyl)-3-({3-oxo-2,3-dihydro- 1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl} methyl)-1,2,5-oxadiazol-2-ium-2-olate	−8.1	−84.66	SER144	HIS163		CYS44, MET49, ASN142, CYS145, MET165.
25	4-(3-methoxyphenyl)-3-({3-oxo-2,3-dihydro-1H- spiro[isoquinoline-4,4′-piperidine]-1′-yl} methyl)-1,2,5-oxadiazol-2-ium-2-olate	−7.5	−80.77	GLY143	HIS41 HIS163		LEU141 ASN142 MET165 GLU166 ASP187 ARG188
26	4-(3-bromophenyl)-3-({3-oxo-2,3-dihydro- 1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl} methyl)-1,2,5-oxadiazol-2-ium-2-olate	−9.8	−90.91	CYS145 SER144.	HIS41 HIS163.		HIS41, CYS44, MET49, ASN142, GLY143, CYS145 MET165 GLU166
27	4-(2,4-dichlorophenyl)-3-({3-oxo-2,3-dihydro-1H- spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−8.0	−61.69		HIS41		MST19 MET165 GLU166 ARG188 GLN189
28	4-(2,3-dimethoxyphenyl)-3-({3-oxo-2,3-dihydro- 1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1, 2,5-oxadiazol-2-ium-2-olate	−8.0	−78.90	GLY143	HIS41		HIS41 CYS44 MET49 ASN142 MET165 GLU166

Molecular dynamics (MD) simulation study

The MD simulations study of Mpro-X77 (Control) and most stable docked complexes obtained as per binding energy and docking scores i.e. M and M was performed by using GROMACS 5.0.5 (www.gromacs.org) program (Abraham et al., 2015) on Linux operating system to examine the structural behavior in terms of stability of M and M complexes in comparison to control under dynamic state. The receptor and ligand molecules from control and docked complexes were separated, and their topology files were built separately. The topology file of receptor Mpro was built using optimized potentials for liquid simulations of all atoms (OPLS-AA) force field (Kaminski et al., 2001). Whereas, topology files of ligand X77 (control), 22 and 26 were generated using an automated online PRODRG server (Van Aalten et al., 1996). For dynamics study the simulation parameters were kept uniform and constant for all three complexes. Initially the solvation of all the complexes was carried out in the system of a cubic box with a size of 1.0 nm by using SPC (simple point-charge) water molecules to provide an aqueous environment. The systems were then neutralized by adding three Na+ ions. The periodic boundary conditions (PBC) were applied in all directions. Energy minimization of the solvated complexes was performed with the steepest descent (SD) method for 50,000 steps at 300 K to remove bad contacts until a tolerance of 1000 kJ/mol was achieved. The linear constraint solver (LINCS) algorithm was used to constrain all bond lengths (Hess et al., 1997). A cutoff of 1.4 nm for the van der Waals interaction and 1.2 nm for electrostatic interactions were used for the simulation of all three complexes. The energy minimization of complexes was then followed by equilibration in two phases. In the first phase of equilibration, an NVT (constant number of particles, volume, and temperature) ensemble with a constant temperature of 300 K was used with a coupling constant of 0.1 picosecond for 500 picosecond, whereas in second phase an NPT (constant number of particles, pressure, and temperature) ensemble with a constant pressure of 1 bar was used with a coupling constant of 2 picosecond for 1 ns equilibration period. For both ensembles of equilibration, the coupling scheme of V-rescale was used and Particle Mesh Ewald (PME) algorithm (Essmann et al., 1995) was used to calculate long-range electrostatic interactions. After equilibration, the production MD run was performed for 100 ns to carry out structural analysis on all complexes. All resulting trajectories obtained from 100 ns MD simulations were analyzed using GROMACS utilities (Abraham et al., 2015). The root mean square deviation (RMSD) was calculated over the 100 ns production simulation for all back-bone atoms of each complex. Similarly, root mean square fluctuations (RMSF), radius of gyration (Rg) relative to the initial structure and H bonds were also calculated for both the complexes and control to assess the stability under dynamic state. MD simulations for all complexes were performed in replicates to check the consistency and reproducibility of the results.

Binding free energy calculation by MM-PBSA

The binding free-energy (ΔGbind) calculations for trajectories obtained from 100 ns MD simulation of complexes viz; Control, M and M were performed using g_mmpbsa tool which implements the Molecular Mechanics Poisson-Boltzmann Surface Area (MM-PBSA) to estimate interaction free energies (Kumari et al., 2014). g_mmpbsa tool represents the ΔGbind of receptor with ligand (inhibitor) in the solvent by using following equations: where, Gcomplex corresponds to total free energy of the complex, Greceptor and Gligand are the total energies of the only receptor and ligand respectively. where, ΔGMM indicates the average molecular mechanics potential energy in vacuum and ΔGsolv represents the free energy of solvation. whereas, ΔGvdw corresponds to Van der Waals energy and ΔGelec represents electrostatic energy.

To calculate binding free energies, the VDW energy, electrostatic energy, molecular mechanics energy, polar solvation energy, and solvent-accessible surface area (SASA) energy were calculated for the trajectories obtained from 100 ns MD simulation of all three complexes. The contribution of each residue in binding free energy was calculated by MmPbSaDecomp.py python script (Kumari et al., 2014). This has been helpful in determining the crucial residues and their thermodynamic state involved in interaction with the SARS-CoV-2 Mpro – inhibitor complexes. The binding energy calculated using g_mmpbsa shows an apparent correlation of 0.80-0.85 with the experimental binding free energy (Kumari et al., 2014).

Results and discussion

Multiple small molecules have been reported to exhibit potent SARS-CoV Mpro inhibition activity (Ghosh et al., 2020). Moreover, Andrographolide (Enmozhi et al., 2020) and resveratrol (Wahedi et al., 2020) were found to be potential inhibitors of SARS-CoV-2 main protease. The 28 different nitric oxide donor furoxan derivatives coupled with benzhydrylpiperazine and spiro-isoquinolino-piperidine from our previous research work (Prabhuling et al., 2020; Pudukulatham et al., 2016) were selected for the current in silico analysis.

Development of small molecules with selective action on the Mpro of the COVID-19 is challenging task as protease inhibition is often correlated with the peptide like structures. The development of imidazole carboxamide derivatives as potent SARS-CoV-2 Mpro inhibitors open the avenue for discovery of small peptide or the peptide like structures as an attractive approach for investigation of potent inhibitors against SARS-CoV-2 Mpro. As selected furoxan derivatives have peptide like structures and due to their structural resemblance with reported imidazole carboxamide (X77) they are virtually analyzed against SARS-CoV-2 Mpro to explore their binding potential.

The results of the predicted pharmacokinetics properties of the analyzed furoxan derivatives are displayed in Figures 1–7 and Table 1. The molecular docking revealed that the spiro-isoquinolino-piperidine class of molecules have more strong binding potentials than benzhydrylpiperazine class of molecules. The docking analysis indicated that the spiro-isoquinolino-piperidine ring is responsible for the anchoring of the molecule in the binding pocket in such a way that other parts of molecule will be adopting the perpendicular orientation in the bind pocket which is responsible for the key interactions in the binding pocket of SARS-CoV-2 Mpro. In case of the benzhydrylpiperazine class of molecules the planner ring orientation of the piperazine pushing interactive phenyl ring in the lower plane and two phenyl ring in the co planner conformation which might be responsible for steric repulsion and subsequently misfitting in the binding pocket of SARS-CoV-2 Mpro.

The crystal structure of the SARS-CoV-2 Mpro bound to potent broad-spectrum non-covalent inhibitor X77 (PDB ID: 6W63) was selected for the docking analysis. The resolution of the selected SARS-CoV-2 Mpro was found to be 2.10 Å and as it is X-ray determined structure, which was found to be suitable for docking process. Ramachandran plot of the selected protein structure was determined to ascertain the structural geometry and to validate the selected protein structure. The Core region of the protein was found to be contain 84.16%, allowed region was found to contain 11.22% of amino acids, which indicated selected protein structure is in good shape and having excellent structural geometry (Figure 1).

Figure 1.

Ramachandran plot of SARS-CoV-2 Mpro (PDB ID: 6W63).

Redocking analysis was performed to validate the docking protocol. The RMSD between docked non-covalent inhibitor X77 and co-crystalized X77 was found to be 0.11 Å which validated the docking protocol. Results of the redocking process indicated that using applied docking protocol to the selected furoxan derivatives will give similar type of binding potential to that of co-crystalized X77 (Figure 2).

Figure 2.

Binding pose of co-crystallized ligand X77 and 22 with crystal structure of SARS-CoV-2 Mpro (PDB ID: 6W63) (UCSF chimera).

Docking analysis of the studied furoxan derivatives were found to have similar binding ability to SARS-CoV-2 Mpro as compared to reported inhibitors. Derivative 4-(3-nitrophenyl)-3-({3-oxo-2,3-dihydro-1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1,2,5-oxadiazol-2-ium-2-olate (22) have binding affinity of −9.6 kcal/mol, and docking score of −90.93, with important hydrogen binding interactions with CYS145 (2.208 A0) and SER144 (2.164 A0) and aromatic interactions with HIS41 (5.281 A0) and HIS 163 (4.628 A0) as shown in Figure 3.

Figure 3.

Binding pose for 22 with crystal structure of SARS-CoV-2 Mpro (PDB ID: 6W63) (Generated using Biopredicta Module of Vlife MDS 4.6).

Derivative 4-(3-bromophenyl)-3-({3-oxo-2,3-dihydro-1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1,2,5-oxadiazol-2-ium-2-olate (26) shown binding affinity of −9.8 kcal/mol and docking score of −90.91 as compared to 22. Hydrogen binding interactions with CYS145 (2.153 A0) and SER144 (1.823 A0) and aromatic interactions with HIS41 (5.175 A0) and HIS 163(4.760 A0) as shown in Figure 4.

Figure 4.

Binding pose for 26 with crystal structure of SARS-CoV-2 Mpro (PDB ID: 6W63) (Generated using Biopredicta Module of Vlife MDS 4.6).

Derivative 4-(4-fluorophenyl)-3-({3-oxo-2,3-dihydro-1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1,2,5-oxadiazol-2-ium-2-olate (23) was shown binding affinity of −8.3 kcal/mol and docking score of −89.44. Important hydrogen binding interactions was observed with GLN189 (1.926 A0), GLY143 (1.565 A0) and aromatic interactions with HIS41(4.521 A0) and HIS 163(4.994 A0) as shown in Figure 5.

Figure 5.

Binding pose for 23 with crystal structure of SARS-CoV-2 Mpro (PDB ID: 6W63) (Generated using Biopredicta Module of Vlife MDS 4.6).

Derivative 4-(4-methoxyphenyl)-3-({3-oxo-2,3-dihydro-1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1,2,5-oxadiazol-2-ium-2-olate (20) displays binding affinity of −8.4 kcal/mol and docking score of −88.64. Hydrogen bond interactions with GLY 143 (1.507 A0), SER144 (2.533 A0) and aromatic interactions with HIS41 (4.323 A0) as shown in Figure 6.

Figure 6.

Binding pose for 20 with crystal structure of SARS-CoV-2 Mpro (PDB ID: 6W63) (Generated using Biopredicta Module of Vlife MDS 4.6).

4-(2-chlorophenyl)-3-({3-oxo-2,3-dihydro-1H-spiro[isoquinoline-4,4′-piperidine]-1′-yl}methyl)-1,2,5-oxadiazol-2-ium-2-olate (24) showed binding affinity of −8.1 kcal/mol and docking score −84.66, along with interaction with SER144 (2.364 A0) via formation of hydrogen bond and with HIS163 (4.797 A0) via formation of aromatic interaction as shown in Figure 7.

Figure 7.

Binding pose for 24 with crystal structure of SARS-CoV-2 Mpro (PDB ID: 6W63) (Generated using Biopredicta Module of Vlife MDS 4.6).

Molecular dynamics simulation study

The MD simulation study of most stable docked complexes viz: M and M helped us to confirm the stable binding ability of furoxan derivatives 22 and 26 towards Mpro receptor structure and compare it with experimentally known inhibitor of Mpro i.e. X77. Thus, MD simulation study provided a clear picture of the overall stability of M and M systems in comparison with Mpro-X77 (Control) (Parulekar & Sonawane, 2018a, 2018b). All systems are subjected to 100 ns MD simulations and their stability in dynamic state was investigated by using GROMACS utilities for calculating RMSD, RMSF, Rg and H bonds (Abraham et al., 2015). The RMSD is a significant parameter used to analyze the stability as per output obtained from MD trajectories which is used to verify the convergence of simulations. RMSD of the protein backbone atoms are plotted as a function of time to check the stability of each complex throughout the simulation. From Figure 8A it is clear that control system reached stable equilibrium after 10 ns whereas, both M and M systems reached equilibrium after 15 ns. The average RMSDs for all simulated systems were calculated and shown in Table 2. These average RMSDs vary between 0.18 nm to 0.28 nm for all systems (Table 2). The RMSD values of Control and M complex were found to be in the range of 0.1 nm to 0.2 nm and remained stable over entire simulation time after reaching equilibrium (Figure 8A). The RMSD values of M complexes are in the range of 0.2 nm to 0.3 nm with slight higher deviation and remained stable throughout 100 ns simulation after attaining equilibrium (Figure 8A). The RMSF parameter in MD simulation is used as a mean describing flexibility differences among residues. The higher value of RMSF generally indicates more flexibility, whereas the low RMSF value indicates limited movements of residues during simulation in relation to its average position. Hence, we analyzed fluctuations of Cα atoms from both docked complexes along with control to study the flexibility of the residues from respective structures in presence of furoxan derivatives 22, 26 and X77 (control). The RMSF values of both the docked complexes along with control are shown in Figure-8B. These RMSF results clearly depicts that residues of Mpro in the presence of 22 and 26 showed similar type of behavior as in presence of X77 (control) [Figure 8B]. Precisely, the active site residues of Mpro i.e. CYS145, SER144, HIS41, HIS163 and GLU166 involved in binding with 22 and 26 showed significantly less fluctuations indicating the stable hydrogen and aromatics interactions of these residues with both the inhibitors as observed for X77 inhibitor used as control (Table 1, Figure 8B). We also determined Rg and calculated average values as reported in Table 2 to understand the level of compactness in the structure of Mpro in presence of furoxan derivatives 22, 26 and X77 (control). Interestingly, the Rg values of both the complexes and control were found to be in the range of 2.2 nm to 2.3 nm and were consistent over entire 100 ns simulation time scale (Figure 8C). Thus, this suggests the similarity in compactness of receptor protein Mpro in the presence of 22, 26 and X77 (control) [Figure 8C]. Additionally, the effective binding ability of furoxan derivatives 22, 26 with SARS-CoV-2 Mpro protein was illustrated by H bonds analysis (Figure 9). Both inhibitors 22, 26 and more prominently 26 showed strong and consistent hydrogen bonding with Mpro in comparison to X77 (control) over entire 100 ns simulation run (Figure 9). Altogether, MD simulation results for both the complexes were in accordance with control system which significantly suggests the effective binding ability of inhibitors 22, 26 as like X77 (control) towards SARS-CoV-2 Mpro receptor structure (Figure 8A–C and Figure 9).

Figure 8.

Analysis of MD simulation trajectories for control and docked complexes in terms of root mean square deviations (RMSD), root mean square fluctuations (RMSF), and radius of gyration (Rg). (A) Backbone RMSDs of Mpro-X77 (control), M and M systems for 100 ns simulation time. (B) RMSF plot of Cα atoms from Mpro receptor structure in presence of X77 (control), 22 and 26 inhibitor molecules (C) Rg of Mpro receptor protein in presence of X77 (control), 22 and 26 inhibitor molecules for 100 ns showing similarity in compactness of SARS-CoV-2 Mpro protein.

Table 2.

The average RMSD and Rg for all the 100 ns simulated systems.

System	Backbone RMSD (nm) with S.D.	Radius of gyration Rg (nm) with S.D.
Control	0.18 nm (0.02)	2.20 nm (0.01)
Mpro-22	0.20 nm (0.02)	2.24 nm (0.01)
Mpro-26	0.28 nm (0.03)	2.25 nm (0.02)

*S.D.= standard deviation.

Figure 9.

The comparative hydrogen bonds analysis of control and both the docked complexes (M and M) over 100 ns simulation run.

Binding free energy calculations and residue contributions study by MM-PBSA

The binding free energy calculations on the basis of different energetic terms for control, M and M complexes confirmed the binding efficiency of X77 (control) and furoxan derivatives 22, 26 at the known inhibitor binding site region of SARS-CoV-2 Mpro structure (Table 3). As seen from Table 3 control showed involvement of higher VDW energies for stabilizing the Mpro-X77 complex whereas, M and M complexes showed significant involvement of higher electrostatic energies for binding thus suggesting prominent role of hydrogen bonds in stabilizing both the docked complexes (Table 3, Figure 9). The control system seems to be more stabilized by non-bonded contacts in comparison to hydrogen bonds which could be the reason for higher VDW energies and lower electrostatic energies observed for control (Table 3, Figure 9). Similarly, polar salvation energy values were also obtained higher for M and M complexes in comparison to control (Table 3). The estimated binding energies (ΔGbind) for control was obtained as −130.014 ± 9.977 whereas, ΔG for M and M was observed to be −183.860 ± 9.106 and −171.972 ± 8.757 respectively (Table 3). Hence, on comparative account 22, 26 showed more effective binding than X77 (control) towards SARS-CoV-2 Mpro receptor on different energetic terms (Table 3).

Table 3.

Binding energy components comparison for SARS-CoV-2 Mpro protein in complex with X77 (control) and most stable inhibitor molecules 22, 26 as obtained from in silico investigation.

Complex	ΔEvdwa	ΔEelec	ΔEMM	ΔGpolar	ΔGnon-polar	ΔGbinding
Control	−168.312 ± 13.357	−7.801 ± 4.448	−176.113 ± 11.457	61.034 ± 4.082	−14.934 ± 0.054	−130.014 ± 9.977
Mpro-22	−123.792 ± 7.327	−151.609 ± 9.399	−275.402 ± 15.326	102.911 ± 7.036	−11.370 ± 0.993	−183.860 ± 9.106
Mpro-26	−113.942 ± 7.049	−245.111 ± 13.287	−359.053 ± 19.781	219.180 ± 13.912	−12.099 ± 0.602	−171.972 ± 8.757

aΔEvdw, ΔEelec, ΔEMM, ΔGpolar and ΔGnon-polar are binding energy components of van der Waals, electrostatic, molecular mechanics, polar and non-polar solvation (SASA) energies, respectively. ΔGbinding is the total binding energy. The unit of energy is kJ/mol.

Furthermore, individual energy decomposition of all residues from SARS-CoV-2 Mpro receptor structure was calculated to investigate the key residue of Mpro involved in binding with X77 (control) and 22, 26 (Figure 10A–C). The results obtained from the calculation of energetic contribution of individual residues clearly confirm the effective binding ability of 22, 26 in comparison to X77 (control) towards Mpro (Figure 10A–C). These results are totally in accordance with the results obtained for the binding free energy calculations of all three complexes on the basis of different energetic terms (Table 3, Figure 10A–C). As seen from Figure 10A–C less residue wise energy contribution is observed for control whereas, significant high values (especially electrostatics and polar solvation) are observed for complexes M and M thus reflecting the strong binding nature of 22, 26 towards SARS-CoV-2 Mpro structure. Residue wise decomposition results confirms that the interacting residues from binding pocket of SARS-CoV-2 Mpro as obtained from docking analysis (Table 1) also shows significantly higher energetic contribution in binding with potent furoxan derivatives 22, 26 in comparison to control (Figure 10A–C).

Figure 10.

The energetic contribution of individual residues from all three complexes to binding energy in kJ/mol for 100 ns MD simulation, (A) all residue contribution of the Mpro-X77 (control) complex in the binding energy (kJ/mol). (B) all residue contribution of the M complex in the binding energy (kJ/mol) and (C) all residue contribution of the M complex in the binding energy (kJ/mol).

Conclusion

Thus, overall in silico analysis performed in the current study reveal the effective binding affinity of investigated spiro-isoquinolino-piperidine-furoxan derivatives in comparison to the experimentally known inhibitor of Mpro i.e. X77 (control) and also they were found to be superior to the corresponding benzhydrylpiperazine derivatives. Further modification of these spiro-isoquinolino-piperidine-furoxan derivatives can lead to discovery of potent SARS-CoV-2 Mpro inhibitors having the advantage of exogenous NO release feature. Although the molecules have shown prominent results in the in silico analysis, the wet lab experiment are required to validate the results. Present results can act as starting point for development of potent selective SARS-CoV-2 Mpro inhibitors.