Identification of a potential SARS-CoV2 inhibitor via molecular dynamics simulations and amino acid decomposition analysis.

Considering lack of validated therapeutic drugs or vaccines against contagious SARS-CoV2, various efforts have been focused on repurposing of existing drugs or identifying new agents. In an attempt to identify new and potential SARS-CoV2 inhibitors targeting specific enzyme of the pathogen, a few induced fit models of SARS-CoV2 main protease (Mpro) including N-aryl amide and aryl sulfonamide based fragments were subjected to a multi-step in silico strategy. Sub-structure query of co-crystallographic fragments provided numerous ZINC15 driven commercially available compounds that entered molecular docking stage to find binding interactions/modes inside Mpro active site. Docking results were reevaluated through time dependent stability of top-ranked ligand-protease complexes by molecular dynamics (MD) simulations within 50 ns. Relative contribution of interacted residues in binding to the most probable binding pose was estimated through amino acid decomposition analysis in B3LYP level of theory with Def2-TZVPP split basis set. In confirmation of docking results, MD simulations revealed less perceptible torsional distortions (more stable binding mode) in binding of ZINC_252512772 (ΔGb -9.18 kcal/mol) into Mpro active site. H-bond interactions and hydrophobic contacts were determinant forces in binding interactions of in silico hit. Quantum chemical calculations confirmed MD results and proved the pivotal role of a conserved residue (Glu166) in making permanent hydrogen bond (98% of MD simulations time) with ZINC_252512772. Drug-like physicochemical properties as well as desirable target binding interactions nominated ZINC_252512772 as a desirable in silico hit for further development toward SARS-CoV2 inhibitors. Highlights A few N-aryl amide/aryl sulfonamide based fragments were subjected to a multi-step in silico strategy to afford potential SARS-CoV2 Mpro inhibitors. MD simulations revealed less perceptible torsional distortions (more stable binding mode) in binding of ZINC_252512772 (ΔGb -9.18 kcal/mol) into Mpro active site. H-bond interactions and hydrophobic contacts were determinant forces in binding interactions of in silico hit. Quantum chemical calculations confirmed MD results and proved pivotal role of a conserved residue (Glu166) in making permanent hydrogen bond (98% of MD simulations time) with ZINC_252512772. Communicated by Ramaswamy H. Sarma.

Introduction

Coronaviruses (CoVs) are enveloped viruses with positive single-stranded RNA genomes infecting both animals and human (Fung & Liu, 2019). CoVs are caused of several diseases from common cold to more severe highly pathogenic ones such as SARS (Severe acute respiratory syndrome), MERS (Middle-east respiratory syndrome) and Coronavirus disease 2019 (Covid-19) with higher mortalities (Cui et al., 2019). A novel coronavirus SARS-CoV-2 (severe acute respiratory syndrome coronavirus 2) is responsible for Covid-19 that first appeared in Wuhan region of China within December 2019 (Zhu et al., 2020). Due to high transmission rate, SARS-CoV-2 spread quickly in multiple countries and the WHO declared an outbreak of COVID-19 a pandemic on March 2020 (https://www.who.int/emergencies/diseases/novel-coronavirus-2019; Gorbalenya et al., 2020). At the time of writing present manuscript, more than 200 countries were affected with more than 2600000 confirmed cases and about 180000 deaths making it a serious threat to global public health (https://www.who.int/emergencies/diseases/novel-coronavirus-2019).

Similar to other members of beta-coronaviruses (SARS-CoV and MERS-CoV), SARS-CoV-2 cause viral pneumonia by attacking lower respiratory system, entry into alveolar epithelial cells, rapid replication and triggering a strong immune response leading to pulmonary tissue damage (Villar et al., 2019). More vulnerable population that must be majorly managed on prevention of SARS-CoV2 infection, are the people with decreased immune function such as the elderly and individuals suffering from basic diseases (Chen et al., 2020).

Given the urgency of highly contagious SARS-CoV2 and considering the lack of validated therapeutic drugs or vaccines against this pathogen, several research groups rely on repurposing FDA approved drugs with proven effectiveness against similar infectious diseases (Boopathi et al., 2020; Liu et al., 2020). Numerous drugs, such as ribavirin, lopinavir and interferon have been tried against SARS or MERS, although the efficacies of some drugs remained controversial (Zumla et al., 2016). A few reports indicated that broad-spectrum antiviral drug arbidol that has been used against influenza outbreaks, could enter the clinical trials (Li & De Clercq, 2020). Optimistic results on the application of chloroquine and remdesivir against clinical isolate of 2019-nCoV have been attained (Wang et al., 2020) and remdesivir is also currently under clinical investigations against Ebola virus infection (Mulangu et al., 2019). Barcitinib was proposed as a potential anti-Covid-19 agents on the basis of its anti-inflammatory effect and possibly reducing viral entry (Richardson et al., 2020). In an open-label control study, Cai et al. demonstrated that favipiravir could control Covid-19 progression and viral clearance (Cai et al., 2020). Another open-label non-randomized clinical trial survey indicated that viral load reduction/disappearance in hydroxychloroquine administered Covid-19 patients were reinforced by azithromycin (Gautret et al., 2020). In addition to synthetic drugs (Adeoye et al., 2020; Lobo-Galo et al., 2020), traditional Chinese medicines (Ren et al., 2020; Tahir et al., 2019) and other medicinal plants (Aanouz et al., 2020; Enmozhi et al., 2020; Islam et al., 2020; Umesh et al., 2020) also draw attention in treatment of Covid-19.

Elucidation of SARS-CoV2 pathogenic mechanisms with the aim of identifying potential drug targets, contribute significantly to successful drug repurposing purposes and also design and development of novel therapeutic agents with probably less side effects (Al-Khafaji et al., 2020). Human CoVs are characterized by structural proteins such as spike protein (Protein S) and also nonstructural proteins including 3-chymotrypsin-like cysteine protease (3CLpro), papain-like protease (PLpro) and RNA dependent RNA polymerase (RdRp) (Ibrahim et al., 2019). During viral infection progression, both structural and nonstructural proteins play a key part, proposing that targeting them could be a promising strategy for controlling COVID-19 infection (Duan et al., 2020). It should be noted that some of the considered proteins were previously proven to be drug discovery targets in the case of SARS-CoV and MERS-CoV (ul Qamar et al., 2020) and may be looked upon as druggable targets in the case of SARS-CoV2 due to the sequence similarity with SARS-CoV (ul Qamar et al., 2020) and also similar binding sites (Joshi et al., 2020; Liu et al., 2020). Hoffmann et al. conducted an interesting research on TMPRSS2 (host-cell produced transmembrane protease serine 2) that facilitates binding of spike protein to angiotensin-converting enzyme 2 (ACE2) on the host cell surface membrane. It was demonstrated that SARS-CoV-2 uses the SARS-CoV ACE2 and TMPRSS2 for host cell entry and clinically proven TMPRSS2 inhibitor camostat mesylate could block SARS-CoV-2 entry (Hoffmann et al., 2020).

Along with conducting in vitro and clinical assessments with the aim of successful drug repurposing, few evaluations were also dedicated to in silico development or computational molecular modeling of new phytochemicals or synthetic compounds (Elfiky, 2020d; Elfiky & Azzam, 2020; Gyebi et al., 2020; Wahedi et al., 2020) or FDA approved drugs against different SARS-CoV2 targets such as RdRp (Elfiky, 2020a, 2020b, 2020c), 3CLpro (Khan, Jha, et al. 2020a; Khan, Zia, et al. 2020b; Muralidharan et al., 2020; Pant et al., 2020), helicase (Beck et al., 2020), hydrolase (Elmezayen et al., 2020), nucleocapsid protein (Sarma et al., 2020), envelope protein (Gupta et al., 2020), cell-surface Heat Shock Protein A5, ACE2 and spike protein (Abdelli et al., 2020; Ibrahim et al., 2020; Sinha et al., 2020). Moreover, some studies focused on the design of a vaccine and a preventative peptidomimetic SARS-CoV2 antagonist (Enayatkhani et al., 2020; Hasan et al., 2020; Robson, 2020). SARS-CoV2 protease enzymes are of significant interest in targeting and modeling studies since they have determinant role in cleaving viral polyproteins into effector proteins, deubiquitinating certain host cell proteins and immune suppression (Baez-Santos et al., 2015; Liu et al., 2020).

Given the vital role of SARS-CoV2 protease in viral infection, current study focused on the identification of in silico hit compounds as potential SARS-CoV2 inhibitors via multi-step modeling strategy comprising docking, molecular dynamics (MD) simulations and ab initio based amino acid decomposition analysis of a ZINC retrieved chemical entities (Sterling & Irwin, 2015). Our structure based simulations run on high-resolution 3D structures of SARS-CoV2 main protease (Mpro) in complex with several fragment ligands (PDB IDs 5R7Y, 5R7Z, 5R80, 5R81, 5R82, 5R83 and 5R84) which were recently unveiled by Fearon et al. (Fearon et al., 2020). To get more information, hierarchical view of the multi-step simulations is depicted in Figure 1.

Figure 1.

Hierarchical view of the multi-step simulations to identify hit compounds as potential SARS-CoV2 inhibitors; Covid-19 related 3 D Mpro structures were screened and subjected to AutoDock4.2 validation step to afford three holo structures with PDB IDs 5R80, 5R81 & 5R84 (Fearon et al., 2020). Corresponding co-crystallographic fragment ligands were subjected to substructure query in Zinc15 database with the criteria set on all subsets including purchasable compounds (Sterling & Irwin, 2015). 295 compounds were selected and entered into molecular docking study to attain three top-ranked analogues with highest binding affinities to related protease targets. Time dependent stability of top-ranked ligand-protease complexes were checked by MD simulations within 50 ns and the cycle might be repeated for next top-ranked candidates if non-desirable results were achieved. Subsequent amino acid decomposition analysis was performed on in silico hit compounds to acquire ligand-residue intermolecular binding energies by functional B3LYP in association with split valence basis set using polarization functions (Def2-TZVPP). Schematic 3 D representation of fragment-enzyme complexes were generated by NGL which is a WebGL based 3 D viewer incorporated into PDB (Rose et al., 2018).

Materials and methods

Target dataset

3D structures of SARS-CoV2 Mpro in complex with several fragment ligands (5R7Y, 5R7Z, 5R80, 5R81, 5R82, 5R83 and 5R84) were retrieved from the Protein Data Bank (www.rcsb.org) (Fearon et al., 2020). Co-crystallographic dimethyl sulfoxide and water molecules were all removed from the original PDB files and non-polar hydrogens merged. Acquired files were converted to pdbqt format after assigning Kollman charges and solvation parameters via AutoDock Tools program (ADT) (Morris et al., 2009; Sanner, 1999).

Ligand dataset

Internal validation of AutoDock4.2 software was done by seven co-crystallographic ligands (5R7Y, 5R7Z, 5R80, 5R81, 5R82, 5R83 and 5R84) using their interacted conformations as starting point. Detailed characteristics of cognate ligands are summarized in internal validation section. All archived ligands that were retrieved in MOL format from Zinc15 database (Sterling & Irwin, 2015) were finally converted to pdbqt using ADT subsequent to the assignment of Gasteiger charges and solvation parameters (Morris et al., 2009; Sanner, 1999). Physicochemical descriptors and ADME parameters were computed by web-based SwissADME program (Daina et al., 2017).

Molecular docking

Lamarckian genetic algorithm (LGA) incorporated into AutoDock4.2 package was used to run ligand flexible docking simulations on Zinc driven small molecules (Morris et al., 2009; Sanner, 1999). All the processing steps and docking parameters were arranged according to previous protocols (Alikhani et al., 2018). Post-docking representation and analysis of ligand-enzyme interactions were performed by the aid of the protein-ligand interaction profiler (PLIP) web server (Salentin et al., 2015). Besides protein-ligand interaction profiler (PLIP) package, discovery studio visualizer (ViewerLite) 3.5 was also used to graphically represent the molecules.

MD Simulation study

The all-atom MD simulations were performed using GROMACS 4.6.5 simulation package (Hess et al., 2008) with standard AMBER99SB-ILDN force fields (Lindorff‐Larsen et al., 2010). The compatible topologies of the ligands were prepared using the General Amber Force Field (GAFF) (Wang et al., 2004) and ANTECHAMBER suite of programs (Case et al, 2005), applying the semiempirical AM1BCC charges (Jakalian et al., 2002). MD simulations were performed according to previous studies (Ebadi et al., 2013; Razzaghi-Asl et al., 2018). Briefly, two 10000 and 5000 steps steepest descent and conjugated gradient minimizations were performed to remove potential clashes between water, ions, and ligand-enzyme in the simulation box. Following initial optimization, the system was equilibrated at 300 K during 500 ps NVT and 1000 ps NPT ensembles. 50 ns MD production was performed in NPT ensemble without any restraint. The temperature and pressure were controlled using the modified V-rescale thermostat from Berendsen (Berendsen et al., 1984) with 0.1 ps time constant and Parrinello-Rahman barostat (Parrinello & Rahman, 1981) with 2 ps time constant for coupling respectively.

Amino acid decomposition analysis

Analyzing of intermolecular binding energy components was performed via constructing a chemical environment including key interacted SARS-CoV2 residues around a docked molecule. N-terminals of interacted amino acids were acetylated while C-terminals were methyl amidated without disturbing original conformational and configurational features. Polar hydrogen bonds were optimized using B3LYP/Def2-TZVPP method through heavy atom fixing (HAF) approximation (constrained optimizations) (Fogarasi et al., 1992). All the ligand-residue binding energies were estimated by the same method and basis set. The whole calculations were performed with the ORCA ab initio quantum chemistry package (Neese, 2012).

Results

To check the validity of LGA in predicting binding modes within SARS-CoV2 Mpro binding site, candidate co-crystalized ligands were extracted from original complex structures and subjected to re-docking into the corresponding induced fit models (5R7Y, 5R7Z, 5R80, 5R81, 5R82, 5R83 and 5R84). For this purpose, root-mean-square deviation (RMSD) of the Cartesian coordinates of re-docked co-crystallographic ligand atoms were used as validation criteria (Vyas et al., 2008). On the basis of re-docking results, adaptable predictability levels (≤ 2 Å) could be achieved with 100 independent GA runs and 2.5 × 107 maximum number of evaluations for complexes 5R80, 5R81 and 5R84 (Table 1).

Table 1.

AutoDock4.2 validation results for different 3 D holo SARS-CoV2 protease-ligand complexes retrieved from PDB.

No.	PDB code	Resolution (Ǻ)	Co-crystallographic ligand	Top-ranked population (out of 50)	ΔGb (kcal/mol)	RMSD from reference (Å)
1	5R7Y	1.65		23	−5.82	6.42
2	5R7Z	1.59		33	−6.05	4.17
3	5R80	1.93		49	−5.02	1.99
4	5R81	1.95		49	−6.82	1.58
5	5R82	1.31		45	−5.16	7.25
6	5R83a	1.58		50	−6.94	6.93
7	5R84	1.83		32	−7.21	0.40

Drug-enzyme inhibition constants (Ki) were estimated via calculation of free binding energies (equation 1) which was incorporated into the force-field of AutoDock4 as a dependent variable (equation 2):

In equation 1, ΔGb is indicative of free binding energy (kcal/mol), R is the gas constant (Cal.K−1.mol−1) and T represents temperature in kelvin (298.15 K). A number 2.71828 is indicates Napier's constant. In equation 2, ΔGTorsional is a term which shows loss of torsional free energy upon binding to enzyme binding site. EvdW and EH-bond and EElectrostatic define hydrophobic interaction energy (van der Waals energy), hydrogen bond interaction energy and electrostatic interaction energies. EDesolvation represents desolvation energy or ‘hydrophobicity’ for drug target interaction. To explain more, EDesolvation represents an energy required to replace ligand-water bonds with that of ligand target ones.

Corresponding co-crystallographic fragments (PDB IDs 5R80, 5R81 & 5R84) were subjected to substructure query in Zinc15 database with a criteria set on ‘all subsets’ to afford 295 purchasable compounds for docking simulations (Sterling & Irwin, 2015). Achieved docked poses were ranked to attain best and worst conformations associated with their binding affinities (Table 2).

Table 2.

AutoDock4.2 results for top-ranked Zinc15 driven compounds with highest binding affinities to SARS-CoV2 protease; within each structure, primary fragment used as a criteria for substructure query is depicted by dashed rectangle/ellipsoid frames.

No.	PDB code	Top-ranked compound (Best conformation)	Top-ranked population (out of 50)	ΔGb (kcal/mol)
				Best conformation	Co-crystallographic conformation	Worst conformation
1	5R80	ZINC_170619726	33	−7.14	−5.02	−6.52 ZINC_74941908
2	5R81	ZINC_108126011	16	−9.39	−6.82	−6.52 ZINC_74941908
3	5R84	ZINC_252512772	20	−9.18	−7.21	−6.77 ZINC_70069

To characterize binding modes of screened compounds inside the active site of SARS-CoV2 Mpro, all intermolecular interactions of ZINC_170619726, ZINC_108126011 and ZINC_252512772 with key residues, schematic representation of ligand binding poses are depicted in Figure 2. It should be noted that for salt bridge interactions, the distance is reported between charge centers.

Figure 2.

Schematic 3 D representation of binding poses of ZINC_170619726, ZINC_108126011 and ZINC_252512772 inside binding site of SARS-CoV2 protease with PDB accession codes a) 5R80, b) 5R81 & c) 5R84, respectively; Yellow spheres within ZINC_170619726 conformation indicates charge centers of salt bridge between SARS-CoV2 and ligand.

MD Simulations

All systems reached steady-state after a few ps from the beginning of MD simulation. The convergence of energy and temperature was evaluated by calculating mean and RSD during the 50 ns MD simulation. The temperature of the holo 5R84, ZINC_170619726, ZINC_108126011 and ZINC_252512772 systems were stable at 300 K (300.0 (RSD: 0.38) in four systems). The energy of holo 5R84, Zinc_170619726, ZINC_108126011, and ZINC_252512772 systems reached −173007 (RSD: 0.16%), −168959 (RSD: 0.16), −170846 (0.16) and −173011 (0.16) kcal/mol, respectively and were stable during 50 ns MD simulation.

The stability of SARS-CoV2 Mpro during MD simulation was evaluated through monitoring the RMSD matrix, the radius of gyration, and intra-enzyme hydrogen bonds. The 3 D structure of SARS-CoV2 Mpro was stable during the 50 ns MD simulation. The RMSD of protein in crystallographic structure (5R84) increased to 0.12 nm in 20 ps at the beginning of production run. The RMSD was stable during 50 ns MD simulation (0.19 ± 0.014 nm (RSD: 7.7%)). The radius of gyration (Rg) was constant during 50 ns MD simulation (2.21 ± 0.009 nm). The mean number of intra-hydrogen bonds in SARS-CoV2 Mpro for crystallographic structure (5R84) was calculated 231.7 ± 7.3.

Figure 3 shows the RMSD matrix and radius of gyration of three potential inhibitors of SARS-CoV2 Mpro found via virtual screening. The maximum RMSD changes in the RMSD matrix of three systems were less than 0.3 nm indicating the convergence of SARS-CoV2 Mpro to equilibrium structure.

Figure 3.

RMSD matrix and radius of gyration plots of Zinc_170619726, ZINC_108126011 (PDB 5R80), and ZINC_252512772 (PDB 5R84) in complex with SARS-CoV2 Mpro.

The mean number of intra-hydrogen bonds in SARS-CoV2 Mpro for Zinc_170619726, ZINC_108126011, and ZINC_252512772 systems were calculated as 227.3 ± 6.9, 226.0 ± 7.2 and 225.6 ± 6.9 in each frame of trajectories, respectively.

MD simulation of holo structure of SARS-CoV2 Mpro revealed that the Mpro formed stable complex with Z31792168 (Figure 4). Z31792168 was stable in the active site of SARS-CoV2 Mpro. Main conformational changes about 9 ns (red line in the RMSD matrix) was due to the rotation of cyclohexane ring and did not change the orientation of Z31792168 in the active site.

Figure 4.

The orientation of Z31792168 in holo structure of SARS-CoV2 Mpro (PDB 5R84) after 6, 10 and 40 ns of MD simulations. Cys145, Glu166, and Thr190 were depicted as blue sticks to represent the dimensions of the active site for better comparison.

The stability of Zinc_170619726 in complex with SARS-CoV2 Mpro was evaluated in terms of RMSD matrix and clustering analysis. The conformational changes in Zinc_170619726 during the 50 ns MD simulation led to 14 clusters by considering 2 nm cutoff and using gromos method. The size of the top three clusters was obtained as 83.8, 6.5, and 4.8% of the total population, respectively. Conformational changes of Zinc_170619726 in the active site of SARS-CoV2 Mpro are depicted in Figure 5. The red color in the margin of the RMSD matrix before 7 ns indicated the major conformational change in comparison to the initial structure. In predicted docking pose, the hydroxyl moiety of Zinc_170619726 formed a hydrogen bond with residue Thr190. But during the first 7 ns of MD simulations, it began to turn and orient toward Glu166. Zinc_170619726 formed an average hydrogen bond of 0.6 ± 0.7 per ns indicating that the hydrogen bonds formed between Zinc_170619726 and SARS-CoV2 Mpro were not relatively stable. No hydrogen bond with Glu166 was detected by analyzing the MD trajectories. The sulfone group of Zinc_170619726 formed a hydrogen bond with Gln192 in 37% of MD simulation time.

Figure 5.

The orientation of Zinc_170619726 (PDB 5R80) in the active site of SARS-CoV2 Mpro for the most populated cluster (83.8%), after 4, and 40 ns of MD simulations. Cys145, Glu166, and Thr190 were depicted as blue sticks to represent the dimensions of the active site for better comparison.

Figure 6 shows that ZINC_108126011 experienced a high level of conformational variation. The clustering of trajectories during the 50 ns MD simulation resulted in 109 clusters. The size of the most populated cluster was 31.9% of the total population. At the time of 8 ns of MD simulation, ZINC_108126011 leaped out of the binding site of SARS-CoV2 Mpro. After a short period (< 4 ns), ZINC_108126011 bound almost in the same pose as the initial docking structure. But in the newly formed complex, quaternary amine was more exposed to the bulk water. The amine group of ZINC_108126011 formed a salt bridge with Glu166 in 31.9% of trajectories. During this time, ZINC_108126011 formed a stable hydrogen bond with Asn142 indicating a different induced fit in the structure of SARS-CoV2 Mpro. In the last 10 ns of MD simulation, the amino moiety was only exposed to the bulk water.

Figure 6.

The orientation of ZINC_108126011 (PDB 5R81) in the active site of SARS-CoV2 Mpro for the most populated cluster (31.9%), after 8, and 16 ns of MD simulations. Cys145, Glu166, and Thr190 were depicted as blue sticks to represent the dimensions of the active site for better comparison.

Conformational changes of ZINC_252512772 in the initial predicted binding mode were less than Zinc_170619726 and ZINC_108126011 (Figure 7). The maximum change of the RMSD matrix was 0.34 nm. ZINC_252512772 formed a permanent hydrogen bond with Glu166 during 98% of MD simulations time.

Figure 7.

The orientation of ZINC_252512772 (PDB 5R84) in the active site of SARS-CoV2 Mpro for the most populated cluster (93.3%), after 8, and 22 ns of MD simulations. Cys145, Glu166, and Thr190 were depicted as blue sticks to represent the dimensions of the active site for better comparison.

The short range electrostatic and vdW interactions between studied ligands and the active site of SARS-CoV2 Mpro is reported in Table 3.

Table 3.

MD estimated van der Waals and electrostatic energies of 5r80 (Zinc_170619726), 5r81 (ZINC_108126011) and 5r84 (ZINC_252512772).

Compound	Electrostatic (kcal/mol)	vdW (kcal/mol)
5r84 (holo structure)	−11.3 (±3.7)	−32.2 (±3.2)
5r80	−11.7 (±5.0)	−32.0 (±2.4)
5r81	−6.6 (±3.6)	32.6 (±3.6)
5r84	−16.9 (±3.8)	−41.4 (±2.4)

Intermolecular binding energy components

On the basis of molecular docking and MD simulations results, ZINC_252512772 exhibited superior binding features and less conformational variations during 50 ns inside active site of SARS-CoV2 Mpro. In order to determine the relative contribution of each interacted Mpro residue in binding to ZINC_252512772, we were convinced to estimate the binding energies of the most probable binding pose of ZINC_252512772 with key residues surrounding active site in B3LYP level of theory with Def2-TZVPP split basis set (Table 4). It should be noted that Cys145 and Thr190 were also taken into consideration as the interacted environment due to MD outputs. All the ligand-residue binding interaction energies (ΔEb) were calculated through equation 3:

Table 4.

Estimated binding energies of ZINC_252512772 inside the binding site of SARS-CoV2 Mpro (PDB code 5R84) in B3LYP level of theory with the Def2-TZVPP as split basis set.

Residue	Binding energy (kcal/mol)
His41	2.51
Pro52	0.13
Phe140	−0.08
Asn142	−3.16
Cys145	−0.35
Met165	5.27
Glu166	−31.66
Arg188	1.98
Gln189	0.90
Thr190	−0.09

ELR is indicative of residue ligand interaction energy while ER and EL stand for electronic energies for residues and ligand, respectively.

Pharmacokinetics prediction

In silico prediction of pharmacokinetics profile was performed through SwissADME calculator (Table 5).

Table 5.

Estimated physicochemical and ADME properties of selected in silico hits (ZINC_170619726, ZINC_108126011 and ZINC_252512772) predicted by swiss ADME.

Compound	Physicochemical properties
	MW	RTBs	HBA	HBD	TPSA (Å2)	logS / Class	Consensus Log Po/w
ZINC_170619726	299.34	4	6	1	92.29	−2.05	1.10
ZINC_108126011	325.47	4	5	2	83.81	−2.55	1.26
ZINC_252512772	355.46	5	3	2	88.32	−3.43	2.12
	Pharmacokinetics
Compound	GIA	BBB permeant	P-gp substrate	CYP1A2 inhibitor	CYP2C9 inhibitor	CYP2D6 inhibitor	CYP3A4 inhibitor
ZINC_170619726	High	No	No	No	No	No	No
ZINC_108126011	High	No	Yes	No	No	No	No
ZINC_252512772	High	No	Yes	No	No	Yes	No

Discussion

Analysis of binding affinity/mode

Growing fragments into drug-like structures via adding layers of rings, substituents and linkers has been regarded as a key approach for drug design (Jacquemard & Kellenberger, 2019). Desirable fragment libraries are appropriate starting points for developing privileged medicinal structures provided that further interactions to hotspots inside a target binding region can be reached upon molecular extension. Analysis of docked poses exhibited that ZINC_170619726, ZINC_108126011 and ZINC_252512772 were located in the binding site of SARS-CoV2 protease (Figure 8).

Figure 8.

Cluster of Zinc driven docked poses (ZINC_170619726, ZINC_108126011 and ZINC_252512772), co-crystallographic fragments (thick stick) and top-ranked poses (thick stick) within binding site residues of SARS-CoV2 protease binding site; PDB IDs a) 5R80, b) 5R81 & c) 5R84.

Assessed compounds participated in H-bond interaction and hydrophobic contacts with 13 hydrophilic and hydrophobic residues of SARS-CoV2 protease. In addition to H-bond interactions and hydrophobic contacts, Zinc_170619726 built salt bridge via carboxylate moiety with His41 imidazole nitrogens (depicted with yellow spheres as charge centers in Figure 2 and the reported distance is between charge centers). His41 is one of the catalytic dyad residues building an active site of the enzyme (Das et al., 2020; Yang et al., 2003). Comparison of interaction models indicated that Zinc_170619726 tended to made more hydrophilic bonds (Table 6) while hydrophobic contacts were dominant for ZINC_252512772 through involvement of cyclic carbon atoms (Table 7). Although more diverse residues were interacted, Zinc_170619726 was associated with lower free binding energy with regard to other three molecules.

Table 6.

Estimated binding interactions of Zinc_170619726 inside the binding site of SARS-CoV2 protease (Chain A; PDB 5R80).

No.	Residue	Interaction type	Distancea (nm)	Hydrogen bond donor angle (°)	Ligand atom/group	Protein atom
1	Met165	Hydrophobic	3.47	–	C19	CB (Side chain)
2	Glu166	Hydrophobic	3.59	–	C7	CB (Side chain)
3	Leu167	Hydrophobic	3.66	–	C13	CD2 (Side chain)
4	Pro168	Hydrophobic	3.97	–	C13	CG (Side chain)
5	Glu166	Hydrogen bond	2.88	105.98	O4	NH (Backbone)
6	Arg188	Hydrogen bond	1.74	122.07	O16H	O (Backbone)
7	Thr190	Hydrogen bond	2.17	144.22	O16	NH (Backbone)
8	Gln192	Hydrogen bond	2.01	131.95	O16	NE2H (Side chain)
9	His163	Salt bridge	5.10	–	Carboxylate	Imidazole Ns

aFor hydrophobic interactions distance is reported between interacting carbon atoms; For H-bonds distance is reported between hydrogen and acceptor atom; For salt bridges distance is reported between charge centers.

bFor Hydrogen bonds donor angle is reported as the angle between donor, acceptor and hydrogen atoms.

Table 7.

Estimated binding interactions of ZINC_252512772 inside the binding site of SARS-CoV2 protease (Chain A; PDB 5R84).

No.	Residue	Interaction type	Distance a (nm)	Hydrogen bond donor angle (°)	Ligand atom/group	Protein atom
1	His41	Hydrophobic	3.10	–	C5	CB (Side chain)
2	Pro52	Hydrophobic	3.01	–	C4	CG (Side chain)
3	Phe140	Hydrophobic	3.75	–	C19	CB (Side chain)
4	Met165	Hydrophobic	3.85	–	C16	CB (Side chain)
5	Glu166	Hydrophobic	3.28	–	C19	CB (Side chain)
6	Asn142	Hydrogen bond	1.95	135.79	O25	NH (Backbone)
7	Glu166	Hydrogen bond	2.16	114.39	N24H	OE1 (Side chain)
8	Glu166	Hydrogen bond	1.96	165.50	N13	NH (Backbone)
9	Arg188	Hydrogen bond	3.60	104.22	N10H	O (Backbone)
10	Gln189	Hydrogen bond	3.56	106.91	N10	NE2H (Side chain)

aFor hydrophobic interactions distance is reported between interacting carbon atoms; For H-bonds distance is reported between hydrogen and acceptor atom; For salt bridges distance is reported between charge centers.

bFor Hydrogen bonds donor angle is reported as the angle between donor, acceptor and hydrogen atoms.

Met165 and Glu166 were found to form H-bond and hydrophobic interactions with all of the three molecules (Tables 6–8) and hence may be regarded as key residues in maintaining relevant complexes (Kumar et al., 2020). Perpendicular orientation of pyridyl ring and piperidinyl ring of ZINC_252512772 in the active site of SARS-CoV2 protease provided two hydrogen bonds with side chain carboxylate oxygen and backbone NH of Glu166 (Figure 2c). Similar H-bond pattern could be detected for ZINC_170619726 between ester carbonyl oxygen and Glu166 backbone NH (Figure 2a). Comparable hydrogen bond interaction was recorded between Glu166 carboxylate and amino methyl substituent of ZINC_108126011.

Docking results proposed that N-aryl amide derivative (ZINC_252512772) could bind to SARS-CoV2 protease with high binding energy (ΔGb −9.18 kcal/mol) and hence might be a good anti-coronavirus candidate for further development strategies. Most of the established interactions were hydrophobic in nature (Table 7) and probably Glu166 was the most significant residue that participated in three key interactions with C19, pyridyl nitrogen and terminal amide nitrogen of ZINC_252512772. Similar results could be found on drug repurposing studies of an anti-asthmatic drug montelukast which showed highly binding affinity to SARS-CoV2 Mpro through hydrophobic contacts with Thr24, Leu27, His41, Phe140, Cys145, His163, Met165, Pro168 and His172 (Wu et al., 2020).

In the case of sulfonamide derivatives (ZINC_170619726 and ZINC_108126011), it seemed that non-rigidified sulfonamide (ZINC_108126011) could bind more tightly to the SARS-CoV2 binding site. Attachment of a primary amine into cyclopentyl ring of ZINC_108126011 provided a good HBD (hydrogen bond donor) and at the same time HBA (hydrogen bond acceptor) site to interact through H-bonding with Glu166 and Leu167. Similar H-bond interaction could not be recorded for ZINC_170619726 due to the lack of hydrogen on sulfonamide nitrogen and also its non-protruding orientation. Moreover, none of the sulfonamide oxygens formed H-bonds with SARS-CoV2 residues.

MD Analysis

Enzymes are flexible macromolecules that can adopt a limited induced-fit models in order to achieve probably best steric and electronic complementary with their ligands. To evaluate the stability of predicted SARS-CoV2 Mpro complexes with regard to the dynamic characteristics of the enzyme, 50 ns MD simulations of three systems (ZINC_170619726, ZINC_108126011 and ZINC_252512772) were conducted in explicit water. The obtained results of three studied systems were compared to the result of 50 ns MD simulation of crystallographic holo structure (PDB ID: 5R84).

The evaluation of temperature and energy during MD simulations indicated that the systems reached to the steady-state and conservation of energy law was satisfied within four systems.

A comparison of MD results indicated that density of polar surface area in the structure of SARS-CoV2 Mpro inhibitor is the dominant factor in stability of ligand-enzyme complex. Distribution of polar surfaces (O and N atoms) in the structure of ZINC_252512772 is more than ZINC_170619726 and ZINC_108126011. The active site of SARS-CoV2 Mpro is exposed to the surface of enzyme and polar hydroxyl and quaternary amine in the structures of ZINC_170619726 and ZINC_108126011 tended to be more exposed to the bulk water than interacting with enzyme active site. This structural feature has been seen in the structure of other known protease inhibitors especially in -navir and -previr classes (Lv et al., 2015). Appropriate but not the most potent interaction with ionized carboxylate of Glu166 eased the binding and also increased the stability of the ligand-enzyme complex.

Comparing the RMSF of three systems during the 50 ns MD simulations showed that ZINC_252512772 formed tight binding with SARS-CoV2 Mpro. The RMSF in residue Glu166 indicated a significant change among three systems. The RMSF of Glu166 in ZINC_252512772, ZINC_108126011 and Zinc_170619726 was 0.08, 0.13 and 0.12 nm respectively. ZINC_108126011 with quaternary amine formed less favorable interactions with Glu166. In comparison, ZINC_252512772, which interacted with Glu166 via its primary amid moiety, had the most binding energy. Interestingly, ZINC_108126011 and Zinc_170619726 had more constructive interactions with hydrophobic residues (His41, Pro52, and Phe140). In this regard, the most noticeable change in the binding pose of Zinc_170619726 was the protruding of methyl benzoate moiety into the SARS-CoV2 Mpro active site (Figure 4). Similar shift in binding pose could be detected for ZINC_108126011 in a way that tetrahydroquinoline ring permeated into the hydrophobic pocket and made the sulfonamide and quaternary amine more water exposable (Figure 5).

The electrostatic interaction decreased in the following order: 5r84 > 5r80 ∼ 5r84 (holo structure) > 5r81. The vdW interaction of 5r84 with the active site of SARS-CoV2 Mpro was more than 5r81, 5r84 (holo structure) and 5r80.

According to obtained results, ZINC_252512772 showed better steric and electronic complementary fitness with the SARS-CoV2 Mpro active site in comparison to ZINC_108126011 and Zinc_170619726.

MD results confirmed that within 50 ns simulations, less conformational changes were met in binding of ZINC_252512772 with regard to ZINC_170619726 and ZINC_108126011. More stable ZINC_252512772-Mpro complex convinced us to determine relative contribution of each interacted Mpro residue in binding to the most probable conformation of the corresponding ligand. On the basis of achieved results (Table 4), five residues (Phe140, Asn142, Cys145, Glu166 & Thr190) were associated with attractive binding forces while His41, Pro52, Met165, Arg188 and Gln189 were associated with repulsive binding forces.

Probably the most favorable result could be achieved in the case of Glu166 which is one of the catalytic dyad residues. Quantum chemical calculations confirmed MD results and indicated pivotal role of Glu166 in making permanent hydrogen bond (98% of MD simulations time) with ZINC_252512772 (ΔEb −31.66 kcal/mol). Our calculations showed that Glu166 was responsible for tight binding of ZINC_252512772 into the Mpro active site. Interaction of Glu166 with a few synthetic and phytochemical agents were previously reported (ul Qamar et al., 2020) and moreover functional role of conserved residue Glu166 within substrate-binding site of SARS-CoV was pointed before (Lin et al., 2004).

Some hydrophobic contacts were associated with lower binding energies (Phe140 − 0.08 kcal/mol) in B3LYP/Def2-TZVPP level of calculation. Second-ranked interacting residues were found to be Asn142 with −3.16 kcal/mol supported binding energy. Modeling studies showed a hydrogen bond between Asn142 backbone NH and terminal amide oxygen of ZINC_252512772 and estimated binding energy for MD simulated complex might be attributed to the stability of this key interaction. Despite prediction of H-bond interaction between Arg188 backbone oxygen and internal amide NH of ZINC_252512772 by molecular docking, QM calculations did not dedicate attractive binding force to the corresponding interaction (1.98 kcal/mol). According to this, we were prompted to estimate the binding energy to Arg188 in docked pose of ZINC_252512772 which was found to be −0.98 kcal/mol. Different binding energies confirmed MD results and indicated that predicted hydrogen bond between Arg188 backbone oxygen and internal amide NH of ZINC_252512772 might not be a stable interaction since considerable torsional deviation occurred within the most probable binding pose of the ligand during MD time.

In silico ADMET prediction

Physiochemical properties have substantial effects on the behavior of compounds within a living system and therefore it is significant to predict ADMET properties through hit identification. Results of estimated parameters indicated that ZINC_170619726, ZINC_108126011 and ZINC_252512772 met the criteria for drug-likeness and followed Lipinski’s rule of five.

Conclusion

A new pathogen (SARS-CoV2) that first emerged in Wuhan region of China in December 2019 spread quickly into all countries of the world. Given an urgent need toward therapeutic drugs or vaccines against the pathogen, seeking for new privileged anti-SARS-CoV2 molecules with good synthetic accessibility or commercial availability is now an important attitude for structure-based hit/lead generation. In the current study potential SARS-CoV2 inhibitors were proposed through in silico analysis of ZINC15 driven commercially available compounds against Mpro as the specific enzyme of pathogen with determinant role in maturation of viral nonstructural proteins (Nsps) and hence viral life cycle. For this purpose, a few co-crystallographic ligands of SARS-CoV2 Mpro with N-aryl amide and aryl sulfonamide based fragments were used to identify new and potential enzyme inhibitors via screening ZINC driven commercially available developed structures. Molecular docking simulations proposed ZINC_108126011 and ZINC_252512772 as superior Mpro binders with free binding energies comparable to a few re-purposed drugs. H-bond interactions and hydrophobic contacts were dominant attractive forces in binding to target. MD simulations on top-ranked docked poses demonstrated less conformational changes (more stable binding mode) for ZINC_252512772 within Mpro active site during 50 ns trajectory. Per-residue binding energy analysis of ZINC_252512772-Mpro complex determined relative contribution of individual interacted residues in binding to the selected hit. Quantum chemical calculations confirmed MD results and proved the pivotal role of Glu166 in making permanent hydrogen bond (98% of MD simulations time) with ZINC_252512772. Obtained results showed that ZINC_252512772 might be a desirable in silico hit for further development of potent SARS-CoV2 inhibitors. Although additional effort remains to be made for in vitro, in vivo and clinical validations, it is hoped that proposed compound could serve as appropriate hit molecules for further development of anti-Covid-19 agents.

Table 8.

Estimated binding interactions of ZINC_108126011 inside the binding site of SARS-CoV2 protease (Chain A; PDB 5R81).

No.	Residue	Interaction type	Distancea (nm)	Hydrogen bond donor angle (°)	Ligand atom/group	Protein atom
1	Met165	Hydrophobic	3.60	–	C6	CB (Side chain)
2	Pro168	Hydrophobic	3.63	–	C17	CB (Side chain)
3	Pro168	Hydrophobic	3.45	–	C16	CD2 (Side chain)
4	Gln189	Hydrophobic	3.75	–	C26	CB (Side chain)
5	Glu166	Hydrogen bond	1.68	154.57	N22	OE1H (Side chain)
6	Glu166	Hydrogen bond	1.91	162.59	N13H	O (Backbone)
7	Leu167	Hydrogen bond	3.60	106.70	N22H	O (Backbone)

aFor hydrophobic interactions distance is reported between interacting carbon atoms; For H-bonds distance is reported between hydrogen and acceptor atom; For salt bridges distance is reported between charge centers.

bFor Hydrogen bonds donor angle is reported as the angle between donor, acceptor and hydrogen atoms.