Structure-based identification of potential SARS-CoV-2 main protease inhibitors.

To address coronavirus disease (COVID-19), currently, no effective drug or vaccine is available. In this regard, molecular modeling approaches are highly useful to discover potential inhibitors of the main protease (Mpro) enzyme of SARS-CoV-2. Since, the Mpro enzyme plays key roles in mediating viral replication and transcription; therefore, it is considered as an attractive drug target to control SARS-CoV-2 infection. By using structure-based drug design, pharmacophore modeling, and virtual high throughput drug screening combined with docking and all-atom molecular dynamics simulation approach, we have identified five potential inhibitors of SARS-CoV-2 Mpro. MD simulation studies revealed that compound 54035018 binds to the Mpro with high affinity (ΔGbind -37.40 kcal/mol), and the complex is more stable in comparison with other protein-ligand complexes. We have identified promising leads to fight COVID-19 infection as these compounds fulfill all drug-likeness properties. However, experimental and clinical validations are required for COVID-19 therapy.Communicated by Ramaswamy H. Sarma.

Introduction

The severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has emerged as a global pandemic as it has affected the entire population (Pang et al., 2020). The genome sequence of SARS-CoV-2 revealed a close similarity to the coronavirus (CoV) that caused an outbreak of severe acute respiratory syndrome (SARS) in 2003 (Asrani et al., 2020; Lai et al., 2020). World Health Organization (WHO) declared COVID-19 as pandemic because more than 70,000 people died and about 2.5 million infected (Chang et al., 2020). After H1N1 (2009), Polio (2014), Ebola in West Africa (2014), Zika (2016) and again Ebola in the Democratic Republic of Congo (2019); WHO has declared COVID-19 as the sixth public health emergency of global distress (Lai et al., 2020). The current situation is rapidly advancing; thus, the ultimate extent and severity of this pandemic remain to be critical. The mortality rate of SARS-CoV-2 is about 3.8% which has been reported as the lowest than other SARS-CoV (10%) and Middle East respiratory syndrome coronavirus (MERS-CoV) (37.1%). However, its infection rate is more than 10 times higher than these CoVs (Ahn et al., 2020). The most characteristic symptom of COVID-19 patients is a high level of respiratory distress, which needs immediate intensive care facility (Lai et al., 2020; Zou et al., 2020). So far, there has been no effective and accurate treatment to cure COVID-19 (Fatima et al., 2020). Hence, there is an urgent need to make significant efforts to develop therapeutic interventions and diagnostic methods to control coronavirus infections (Asrani et al., 2020).

CoVs belong to the family of Coronaviridae containing a single-strand of positive-sense RNA viruses. These viruses can be classified into four genera: alpha, beta, gamma, and delta. The current SARS-CoV-2 belongs to the beta genus and is commonly known to infect humans (Menachery et al., 2015; Ahmed et al., 2020). The genome length of this virus is about 27-32 kb encoding both structural and non-structural proteins (Zhang & Holmes, 2020). Among the structural proteins, membrane, envelope, nucleocapsid, and spike proteins contribute significantly to virus transmission and its replication in the host cells (Naqvi et al., 2020; Shanmugaraj et al., 2020). PLpro and 3CLpro proteases are vital in virus replication and hence considered a promising drug target. Structural and mechanistic information will help to design potent and selective inhibitors of 3CLpro and PLpro that will eventually be implicated to address COVID-19 as these proteases are indispensable for virus assembly and replication (Zhang et al. 2020).

The unexpected emergence of COVID-19, have underscored the urgent need for effective preventive and therapeutic measures for the development of antiviral therapy (Kumari et al., 2020). To address the critical health challenge of COVID-19, many therapeutic targets for the development of vaccines and drugs are under the explorative stage (Asrani & Hassan, 2020; Wrapp et al., 2020). COVID-19 is difficult to control because of heterogeneous member's periodical cycle and significant overlap with human and wild animal ecologies. Currently, there are no approved specific antiviral therapies available for SARS-CoV-2. Several attempts have been made using approved antivirals drugs (ribavirin and lopinavir-ritonavir) and immunomodulators (corticosteroids, interferons, etc.) (Grein et al. 2020; Li et al., 2020; Morse et al., 2020; Wang et al., 2020). However, these approaches are not much effective in the case of COVID-19 because of significant differences in the surface, structural, and enzymatic proteins of SARS-COV-2 in comparison to the SARS-CoV and MERS-CoV (Grein et al. 2020; Kumari et al., 2020). Design and development of SARS-COV-2-specific direct-acting antiviral drugs could be potentially obtained by targeting conserved enzymes such as 3 C-like protease, papain-like protease, non-structural protein 12 (nsp12), and RNA-dependent RNA polymerase (Zumla et al., 2016; Tang et al., 2020).

The main protease (Mpro)is a key enzyme that plays a critical role in the SARS-CoV-2 life cycle as this enzyme influence the process of viral replication and transcription (Jin et al.; Naqvi et al., 2020), making it an attractive target for antiviral drugs discovery (Ahn et al., 2020). Here, we have used SARS-CoV-2 Mpro as a target enzyme to facilitate the rapid search of antiviral compounds with clinical potential in the therapeutics of COVID-19. Recently, the three-dimensional structure of SARS-CoV-2 Mpro has been determined which is highly similar to that of the SARS-CoV Mpro, with a 96% sequence identity with the RMS deviation of 0.53 Å (Zhang et al. 2020). Although, a minor variation in the sequence and structure of Mpro from both the viruses, a significant difference was noticed in the binding pattern and inhibition of Mpro by inhibitors designed for SARS-CoV-2. Hence, approved protease inhibitors such as disulfiram, lopinavir, and ritonavir have been reported to be active against SARS and MERS CoVs, but are ineffective in the case of SARS‐CoV-2 (Kandeel & Al-Nazawi, 2020).

The bound inhibitor showed the formation of several interactions with the binding pocket residues. The N atom of lactam moiety, N atom of amide functional group of the neighboring His164, and O atom of benzyloxy group formed hydrogen bonds with His163, His164, and Gly143, respectively. The alkyl moiety of leucinamide subunit of N3-ILP is surrounded by the hydrophobic side chain of His41, Met49, Tyr54, Met165, and Asp187 of Mpro. The solvent-exposed area of the Val side chain of the inhibitor indicated that this site could be substituted by a wide range of functional groups. The Ala side chain N3-ILP was surrounded by the side chains of Met165, Leu167, Phe185, Asp192, and Asp189, which formed a hydrophobic pocket. P5 makes van der Waals contacts with P168-A and the backbone of residues 190–191, Figure 1.

Figure 1.

Structural representation of COVID-19 Mpro enzyme in complex with N3-ILP inhibitor; A. Close up view of the binding pocket of Mpro accommodated N3-ILP, B. Close up view of active residues interacting with the inhibitor.

Here, the structural features of already known peptide-like N3-ILP inhibitor of SARS-COV-2 Mpro enzyme have been employed to generate a small library of similar compounds to gain a short‐term and specific solution to treat COVID‐19 patients (Shamsi et al., 2020; Ton et al., 2020). To address this general challenge, we developed an integrated approach of drug discovery using pharmacophore modelling, screening of PubChem library, molecular docking, and molecular dynamics (MD) simulation to find potential preclinical leads for the therapeutic management of COVID-19 (Mohammad et al., 2020). This strategy will provide a route-map for tailoring the antiviral inhibitors that might help in assisting and developing novel SARS-CoV-2 Mpro lead compounds to be implicated in the treatment of COVID-19 patients.

Methods

System preparation

The crystal structure of SARS-CoV-2 Mpro in complex with a peptide-like inhibitor was retrieved from the Protein Data Bank (PDB ID: 6LU7) (Jin et al., 2020). The structure of the enzyme was pre-processed, minimized, and refined using the Protein Preparation Wizard implemented (Madhavi Sastry et al., 2013) in Schrödinger suite (Schrödinger Release 2020-1: Protein Preparation Wizard; Epik). This involved eliminating crystallographic waters, missing hydrogens/side-chain atoms were added to assign appropriate charge and protonation state to relieve the steric clashes among the residues. The OPLS-2005 force-field was used for the energy minimization using a root mean square deviation (RMSD) cut-off value of 0.30 Å. The preparation of the ligand, N3 peptide-like inhibitor (N3-IPL) and the studied compounds were performed using LigPrep (Schrödinger Release 2020-1: LigPrep) module of Schrodinger Suite which performs addition of hydrogens, adjusting realistic bond lengths and angles, ionization states, correct chiralities, stereo chemistries, tautomers, and ring conformations (Fakhar et al., 2020). Partial charges were assigned to the structures using the OPLS-2005 (Harder et al., 2016) force-field, and the resulting structures were subjected to energy minimization until their average RMSD reached 0.001 Å. Epik ionization tool (Schrödinger Release 2020-1: Epik) was used to set the ionization state at physiological pH = 7.4.

Preparation of inhibitor-like ligand library

We have filtered and retrieved 409 compounds from the PubChem database considering Lipinski's rule of five (Lipinski, 2004). All the compounds were selected based on 80% structural similarity and the main scaffold of N3-inhibitor and subsequently considered for further virtual screening analysis.

Identification of pharmacophore hypotheses

For the structure-based pharmacophore modeling, PHASE module (Dixon et al., 2006) implemented in Maestro 11.6 was used with the default set of six chemical features: hydrogen bond acceptor (A), hydrogen bond donor (D), hydrophobic contacts (H), negative ionizable (N), positive ionizable (P), and aromatic ring (R) to construct the most representative features of the Mpro active sites. The five 3 D-features were generated using Hypothesis Generation for Energy-Optimized structure-based pharmacophores considering the excluded volumes within 5 Å of refined ligand for the enzyme (Loving et al., 2009; Salam et al., 2009). Pharmacophore features were selected based on essential interaction contacts with the key residues of the enzyme accommodated the inhibitor. The resulted pharmacophore features contain the functional groups involved in their bioactivity of a targeted enzyme. The Excluded volumes include all atoms within 5 Å of the refined ligand for the target.

Screening of Mpro inhibitors

The obtained five 3 D-pharmacophore features were exported and set as a reference for PHASE-based virtual screening to screen the library of 409 compounds from the PubChem database which was retrieved and filtered with 80% structural similarity to the known N3-Inhibitor of Mpro enzyme (Dixon et al., 2006). Out of 409 compounds, 171 were generated based on the highest PHASE screen score and matched ligand sites. Both the quantity and quality of a similar feature was taken into account in the Phase-Screen-Score factor.

Docking-based virtual screening

Molecular-docking-based virtual screening was performed using the Glide-Based virtual screening workflow of Maestro 11.6 to identify suitable compounds that strongly bind to Mpro enzyme (Halgren et al., 2004). The receptor grid was generated as center coordinates (X = −10.81 Y = 12.41 Z = 68.93) using two cubical boxes having a common centroid to organize the calculations: a larger enclosing and a smaller binding box with dimensions of 24 × 24 × 24 Å and 32 × 32 × 32 Å, respectively. The grid box was centered on the centroid of the ligand in the complex, which was sufficiently large to explore a larger region of the enzyme structure. The ligands were docked using three docking protocols which starts with “High throughput Virtual Screening” (HTVS) followed by “Standard Precision” (SP) and then by “Extra-Precision” mode (XP). Finally, 171 input compounds were evaluated using Docking-Based Virtual Screening and filtered to the final 20 compounds based on the docking scores and XP-GScores.

Admet properties assessment

The QikProp 5.6 module (QikProp) of Schrodinger was used to predict absorption, distribution, metabolism, excretion, and toxicity (ADMET) properties of the considered compounds to generate the ADMET related descriptors. This protocol predicts significant physicochemical and pharmacokinetic-based descriptors of the compounds based on Lipinski’s rule of five (Lipinski, 2004). ADMET properties of the top five compounds and N3-Inhibitor were assessed and analyzed using QikProp 5.6 module and these compounds were designated for the final molecular dynamic (MD) simulations analysis.

MD Simulations

To understand the physical basis of the structure and function of biological macromolecules, MD simulation substantiates to be an essential approach (Amir et al., 2020; Beg et al., 2019; Dahiya et al., 2019; Gulzar et al., 2019; Gupta et al., 2019; Gupta et al., 2019). This technique assists in discovering the structural dynamics and how it is coupled to the biomolecular function of the enzyme (McCammon et al., 1977; Karplus & McCammon, 2002; Arnittali et al., 2019). AMBER 18 package (Case et al., 2018) was used to execute 100 ns MD simulations on all the prepared complexes using (Graphics Processing Unit) GPU accelerated version of Partial Mesh Ewald Molecular Dynamics (PMEMD) simulations (Lee et al., 2018). The ff99SB (Hornak et al., 2006) and the general AMBER force fields (GAFF) (Wang et al., 2006) were employed to parametrize the enzyme and the considered ligands using LEaP module implemented in Amber 18 (D.A. Case). The ANTECHAMBER (Wang et al., 2001) module was used to assign atomic partial charges for the ligands employed in General Amber Force-Field (GAFF) (Wang et al., 2004). The system was solvated using the TIP3P (Harrach & Drossel, 2014) explicit water in a cubic box with an 8 Å box edge. The Na + counter ions were added randomly to neutralize the complex. The partial Mesh Ewald (PME) (Harvey & De Fabritiis, 2009) method was used to account for the long-range electrostatic forces using a cutoff of 12 Å, and the SHAKE algorithm (Ryckaert et al., 1977) was used to constrain all the hydrogen atom bonds. Energy minimizations were performed in two stages with 2500 steps of steepest descent minimization followed by 2500 of the conjugated gradient to remove the bad contacts. The first stage was followed with a harmonic restraint of 500 kcal mol-1 A-2 on the solute molecule whereas, ions and water molecules were relaxed. In the second stage of minimization, the restraints were removed and the whole system was relaxed. Each minimized complex was then gradually heated up from 0 K to 300 K for 200 ps to keep the solute using a weak harmonic restraint of 10 kcalmol-1 A-2. The 50 ps density equilibration with weak restraints followed by the 500 ps constant pressure equilibration at 300 K were performed at constant pressure using Berendsen barostat (Lin et al., 2017). Ultimately, the production phase of 100 ns MD simulation was performed on all the complexes at a constant temperature of 300 K and constant pressure at 1 atm.

Post-dynamic analyses

The atomic coordinates of enzymes bound with the inhibitors were further saved after every 1 ps and the trajectory curves were calculated using the CPPTRAJ module integrated into AMBER 18 package (Roe & Cheatham, 2013). The root means square deviation (RMSD) of Cα atoms, root means square fluctuation (RMSF) of each residue in the complex, a radius of gyration (Rg), solvent accessible surface area (SASA), intramolecular and intermolecular hydrogen bond formation and thermodynamic calculations of all the systems were calculated. Origin software was used for MD trajectories analysis (Janert, 2009).

Binding free energy calculations

The relative binding free energies were calculated using Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) binding free energy (Wang et al., 2019). All water molecules and counterions were stripped using the CPPTRAJ module. The binding free energies (ΔGbind) were calculated with the MM/GBSA method for each system as below:

The free energy term, ΔGbind is computed using the following equations:

Where,

The gas-phase energy (ΔEgas) is the sum of the internal (Eint), van der Waals (EvdW) and Coulombic (Eelec) energies, (Eq. 4). The solvation free energy is the combination of polar (GGB) and nonpolar (Gnonpolar, solvation) contributions (Eq. 5). The polar solvation GGB contribution was calculated using the Generalized Born (GB) solvation model with the dielectric constant 1 for solute and 80.0 for the solvent (Onufriev & Case, 2019). However, the nonpolar free energy contribution was estimated using Eq. 6, where the surface tension proportionality constant, γ, and the free energy of nonpolar solvation of a point solute, β, were set to 0.00542 kcal mol−1 Å−2 and 0 kcal mol−1, respectively. SASA was calculated by a linear combination of the pairwise overlap (LCPO) model.

Result and discussion

Database screening

The crystal structure of SARS-CoV-2 Mpro (PDB ID: 6LU7) was extensively analyzed to understand their binding affinity, mode of binding and interacting residues (Figure 1). Structure analysis revealed the existence of 8 hydrogen bonds offered by Gly143, His163, His164, Glu166, Gln189 and Thr190 of Mpro with the N3-Inhibitor-Like-Peptide (N3-ILP) . The structure of N3-ILP was selected to screen the PubChem database. 409 compounds were filtered and collected based on 80% structural similarity to N3-ILP and Lipinski's rule of five (Lipinski, 2004) from the PubChem database (Kim et al. 2016). The screened 409 compounds are listed in Table S1.

Structure-based pharmacophore modeling

Structure-based pharmacophores derived from the three-dimensional structure of a target protein provide detailed and accurate information on ligand binding (Langer, 2010). The commonly used descriptors in the pharmacophore modeling are H-bond acceptors, H-bond donors, positive and negative ionizable groups, lipophilic regions and aromatic rings. The best 3 D structure-based e-pharmacophores (Loving et al., 2009; Salam et al., 2009) were generated using the receptor-ligand pharmacophore generation protocol implemented in PHASE (Dixon et al., 2006), based on the crystal ligand inside the active site and residues involved in ligand binding. The generated e-pharmacophore of the considered enzyme showed five main 3 D-features including, H-bond acceptor, H-bond donor, and aromatic rings. In each pharmacophore model, the red arrows represent hydrogen bond acceptor, the blue arrow represents hydrogen bond donor and the orange spheres represent an aromatic ring. The 3 D pharmacophore features and 2 D-chemical structure of N3-ILP are illustrated in Figure 2 showing five 3-D pharmacophore features as three hydrogen-bond donor, one hydrogen-bond acceptor and one aromatic ring sphere.

Figure 2

Structure-based pharmacophore modeling based on the analysis of N3-ILP interactions with the Mpro. 3 D pharmacophore features of N3-ILP in complex with Mpro and its chemical structure surrounded by excluded volumes. Red arrow: Hydrogen bond acceptor, blue arrows: Hydrogen bond donor, orange: aromatic ring.

Screening of compounds library

The obtained structure-based pharmacophore hypotheses were used to screen the PubChem database to explore new compounds. These compounds were screened considering the PHASE screen score and matched ligand site factors. A total of 409 compounds passed this filter and subsequently, 171 compounds were retrieved based on the created pharmacophore hypothesis. Molecules that have satisfied all the features of the pharmacophore model were considered as potential hits. The virtual screened compounds are presented in Table S2.

Docking-based virtual screening analysis

A total of 171 compounds obtained from virtual screening were docked using the Glide module (Friesner et al., 2006) of the Schrödinger package (Schrödinger Release 2020-1: Maestro) into the active site of Mpro. A stepwise filtering protocol was used, in the first stage, compounds were docked using HTVS where a total of 157 hits were obtained, Table S3. These compounds were further docked with Glide SP where a total of 81 hits were obtained (Table S4). Afterward, the hits from the previous step were subjected to Glide XP docking and only one pose per ligand was retained. Finally, a total of 20 hits were obtained as shown in Table S5.

Among the 20 hit compounds, the known inhibitor N3-ILP (EC50: 16.77 mM) (Figure S1) and top five compounds with PubChem IDs 54456426, 54152887, 54035018, 91366909, 57076946 with Glide GSCORE (XP SCORE) was estimated as −4.65 kcal/mol, −9.05 kcal/mol, −9.01 kcal/mol, −8.95 kcal/mol, −8.80 kcal/mol and −8.70 kcal/mol indicating a strong binding affinity. The docking based RMSD values for N3-ILP and other top-five compounds were 2.77 Å, 2.64 Å (54456426), 2.68 Å (54152887), 2.02 Å (54035018), 2.72 Å (91366909) and 2.81 Å (57076946) respectively. The criteria for selecting these top five compounds were based on the binding scores (above −8.50 kcal/mol), lowest RMSD values and most promising binding interactions generated after molecular docking calculations. The theoretical binding affinities of identified hits are better than well-known inhibitor N3-ILP, indicating potential future clinical use of these preclinical leads. The number of hydrogen-bonded interactions is a key factor for tighter and specific binding affinity. The interacting binding residues of Mpro enzyme forming non-covalent interactions with compounds 54035018, 54152887, 54456426, 57076946, and 91366909 are illustrated in Figure 3A and B.

Figure 3.

A: Docked poses of Mpro enzyme with the proposed inhibitors. Binding mode of and interactions of Mpro enzyme with the compound with PubChem ID 54035018 and 54152887. B: Docked poses of Mpro enzyme with the proposed inhibitors. Binding mode of and interactions of Mpro enzyme with the compound with PubChem ID 54456426, 57076946 and 91366909.

Compound 54035018 interacts with the binding site residues, Ser144, Glu166, Thr190 and Gln192 through four hydrogen bonds. Compound 54152887-Mpro formed similar four hydrogen bonds to Ser144, Cys145, Glu166, and Thr190. It also forms a pi-pi stacking with His41. The complexes 54456426-Mpro, 57076946-Mpro, 91366909-Mpro are stabilized by H-bonded interactions offered by His41, Leu141, Asn142, Ser144, Glu166 and His164.

Admet analysis

Pharmacokinetic and toxicity properties were predicted using the QikProp module of Schrodinger (QikProp) for N3-ILP, 54035018, 54152887, 54456426, 57076946 and 91366909 compounds. The results of pharmacokinetic and toxicity properties analysis are presented in Table 1. The selected properties are representatives of influence metabolism, cell permeation, bioavailability, and toxicity.

Table1.

In silico ADMET predictions of the selected compounds.

Compounds	aCNS	bQPlogKhsa	cSASA	dQPlogPo/w	eQPlogS	fQPlogBB	g% Human OralAbsorption	hRule OfFive
N3-Inhibitor	−2	−0.51	1150.17	3.14	−6.68	−3.28	51.17	2
54035018	−2	−0.53	700.82	3.45	−5.10	−2.23	82.74	0
54152887	−2	0.24	758.42	3.67	−5.26	−2.03	87.01	0
54456426	−2	0.37	632.56	0.78	−2.62	−2.60	53.92	0
57076946	−2	0.53	697.02	2.78	−3.42	−1.24	79.13	0
91366909	−2	0.13	668.87	4.14	−5.10	−1.50	96.62	0

Predicted central nervous system activity from –2 (inactive) to +2 (active).

Prediction of binding to human serum albumin (acceptable range: -1.5-1.5).

Total Solvent Accessible Surface Area: SASA (acceptable range: 300–1000).

Predicted octanol/water partition coefficient (acceptable range: −2–6.5).

Predicted aqueous solubility, S in mol/dm − 3 (acceptable range: −6.5–0.5).

Predicted brain/blood partition coefficient (acceptable range: –3.0 – 1.2).

Predicted percentage human oral absorption (<25% is poor and >80% is high).

Number of violations of Lipinski’s rule of five, Compounds that satisfy these rules are considered druglike (maximum 4).

The predicted central nervous system activity (CNS) of all the molecules depicted as inactive. The recommended binding to human serum albumin (QPlogKhsa), all compounds showed in the acceptable range, although the N3-ILP showed slight out of range. The estimated total solvent accessible surface area (SASA) of all PubChem compounds met the acceptable range: 300–1000 whereas the known N3-ILP showed out of range. The predicted octanol/water partition coefficient (QPlogPo/w) showed in the acceptable range from −2 to 6.5. The predicted aqueous solubility (QPlogS) and predicted brain/blood partition coefficient (QPlogBB) for all these compounds showed in the acceptable ranges whereas the known inhibitor was out of range. The percentage of human oral absorption for all the compounds met in the recommended range. Some violations of Lipinski’s rule of five for all the PubChem compounds significantly satisfied this rule with zero value than known inhibitor with the rule of five equal to two.

Post-dynamics trajectories analysis

Alternations inside the enzyme structure are directly associated with their biological activities. Any changes or disruption on enzymes' structural integrity might have a significant influence on its activity (Khan et al., 2019). The binding of small molecule inhibitors affects the mechanism of action of enzymes that are involved in disease pathways, thus there is a necessity to calculate the structural dynamics and conformational changes linked with the inhibitory activity of these inhibitors (Skjaerven et al., 2011). The calculation of a time variable concerning an RMSD of Cα atoms from produced trajectories was performed to discover the constancy and efficacy of the simulated Mpro in complex with N3-ILP and along with the top five hits. The 2 D structure of the top five PubChem compounds and N3-ILP-inhibitor are presented in (Figure S2).

The perturbations in the RMSD values as indicated in the plot (Figure 4A) during simulation time revealed a probable conformational deviation in the structure of the enzyme upon inhibitor binding. As Figure 4A revealed, all the systems were stabilized and achieved convergence after almost 30 ns of the simulation run. 54035018-Mpro exhibited the lowest average RMSD of 2.03 Å, while 54152887-Mpro and 54456426-Mpro demonstrated the average RMSD of 2.45 Å and 2.33 Å respectively. The N3-ILP-Mpro indicates a perturbation of 2.50 Å as shown in the plot. This analysis suggests that any other investigations made on the generated trajectories of all models were reliable. The RMSD plots suggest that 54035018-Mpro, 54152887-Mpro, and 54456426-Mpro complexes exhibit the lowest divergence of Cα backbone atoms indicating that binding of these three compounds imposed higher stability on Mpro enzyme relative to N3-ILP-Mpro complex. This could be implicated that the overall inhibitory activity of Mpro by 54035018, 54152887, and 54456426 compounds was illustrated by the stabilization of its structural conformation. The stability of compounds 54035018, 54152887, and 54456426 was further retrieved during the simulation time as their corresponding activities might inform the respective interactions that they cause with active site residues.

Figure 4.

Structural dynamics of Mpro enzyme-ligand complexes. A: C‐α backbone RMSD in Å of all the selected compounds bound to Mpro enzyme; B: Rg values after compound binding.

To further validate our results, we have calculated Rg value, a parameter directly associated with the overall conformational changes in the structure of the enzyme upon ligand binding. It also reveals the structure stability, compactness and folding behavior (Fakhar et al., 2017). We measured the compactness of all the selected compounds and the reference complex by calculating their Rg values. The average values of Rg for 54035018-Mpro, 54152887-Mpro, and 54456426-Mpro complexes were noted to be 38.87 Å, 40.39 Å, and 43.68 Å, respectively. Figure 4B plots revealed a slight alteration in the compactness of the three compounds. The compound showed the lowest Rg as compared to the other two complexes and also with the control N3-ILP-Mpro complex (41.26 Å), thus suggesting increased compactness and better binding with the Mpro enzyme. All these patterns of conformational evaluation are indicating greater stability, flexibility, and compactness of compound 54035018 with the Mpro enzyme.

After analyzing the conformational binding, we have also explored SASA to determine the function of hydrophobic and hydrophilic residues and forces exposed to the solvent during the simulation time (Khan et al., 2018). The fast and precise calculation of SASA is very beneficial in the energetic evaluation of biological molecules. The interactions between the hydrophobic native contacts within the enzyme structure is a significant intermolecular interaction that influences enzyme inhibition. Hydrophobic interaction produced among the non-polar residues confirms the stability of the enzyme structure in solution by protecting the non-polar residues inside the hydrophobic core distant from an aqueous solution (Chen & Panagiotopoulos, 2019; Gupta et al., 2020). As shown in Figure 5A, an average SASA for all the selected compounds has been calculated during the 100 ns MD simulation run. The average value of SASA for the compound 54035018-Mpro complex was 13745 Å2 which was exposed to the solvent system. A total mean of SASA of 14283 Å2 and 14439 Å2 was noted by 54152887-Mpro complex and 54456426-Mpro complexes, respectively. The variations in SASA values for all the complexes during the simulation time resemble the folding and unfolding of an enzyme. The overall SASA values in the control complex were 14009 Å2, higher than the 54035018-Mpro complex. The SASA estimation observed in compound 54035018 bound complex further confirmed that this compound has increased exposure to solvent and subsequently preferred the increased inhibitory activity of compound 54035018 over other complexes.

Figure 5.

Showing SASA and fluctuations in the backbone atoms of the Mpro enzyme shows a 100 ns MD simulation period. A: SASA values calculated during the 100 ns of MD trajectories, B. Values of RMSF in Å plotted against residue number for all the selected compounds bound to Mpro enzyme, and C. Fluctuation in the structure of Mpro enzyme after compound binding.

The flexible or rigid residues of a protein contribute specifically to the biological function of the target enzyme in different biological pathways. Thus, the binding of inhibitors to the enzyme may be analyzed through the change in flexibility in terms of RMSF values (Martínez, 2015). To discover the rigidity and flexibility of residues in Mpro upon binding of the selected compounds, RMSF values for Cα atoms were calculated from trajectories generated over 100 ns of MD simulations production run. As shown in Figure 5B and C, the 54035018-Mpro complex exhibited the least fluctuations in the residues with 11.12 Å. An average RMSF of 15.06 Å and 15.29 Å were noticed in complex 54152887-Mpro and 54456426-Mpro, respectively. The complex, N3-ILP-Mpro revealed an average of 11.25 Å that is slightly greater than 54035018-Mpro complex, suggesting a better binding in comparison to the N3-ILP-Mpro complex. This significant decrease could be associated with structural inactivation that is ensured as an outcome of the influential binding of this compound in the catalytic site of the Mpro enzyme. The limited fluctuation of residues could have preferred Mpro enzyme inhibition through compound 54035018.

Hydrogen bond analysis

For the overall conformation and stability of enzyme structure, we have assessed intramolecular and intermolecular hydrogen bond analysis (Figure 6). This analysis gives deeper insights into the binding mechanism of enzyme-ligand with specific consideration. An average number of intramolecular hydrogen bonds in 54035018-Mpro and 54152887-Mpro complexes were observed to be 135 and 141, respectively as shown in Figure 6A. As we have mentioned in the residue interaction map analysis these complexes contributed greatly to the binding throughout the simulation period. The number of intermolecular hydrogen bonds formed in the active site of the Mpro enzyme noted to be 2 and 3 hydrogen bonds in 54035018-Mpro and 54152887-Mpro complexes and between 3-4 hydrogen bonds in 54456426-Mpro and 57076946-Mpro complexes with higher fluctuations Figure 6B.

Figure 6.

Hydrogen bond analysis. A. Intramolecular and B. Intermolecular hydrogen bonds in Mpro enzyme with the selected compounds calculated after 100 ns MD simulation.

To understand the intermolecular interactions between ligand with the active site residues, we have calculated the ligand interaction maps (Figure 7A and B). Compound 54035018 forms two direct hydrogen bonds with the polar residues Thr26 and Asn142 with OH and NH group suggestive of a stable interaction. However, hydrogen bond formation of Asparagine was noted to be absent in other complexes. Another polar residue Gln192 formed two hydrogen bonds with the hydroxyl functional group of compound 54152887. Gln166 also formed a hydrogen bond with the hydroxyl group of 54152887 compound and residue His41 formed a π-π stacking within the complex. Compounds 54456426 and 57076946 also contributed significantly to the hydrogen bond network by forming major hydrogen bonds between the polar (Gln19, Thr26, Thr190), negatively charged (Gln166) and hydrophobic (Cys145) amino acid residues indicative of a stabilized binding with the enzyme. His41 forms a π-π stacking in compound 54456426 and 91366909. Compound 91366909 formed only one hydrogen bond with the negatively charged Gln166 residue, thus suggestive of weak binding. The hydrogen bond network plays an important role in the drug-binding mechanism, thus Thr amino acid residue contributed notably in the binding of the Mpro enzyme as shown in Figure 7.

Figure 7.

A: Showing H-bonded interaction map of Mpro enzyme with the selected compounds calculated after 100 ns MD simulation. B: Showing H-bonded interaction map of Mpro enzyme with the selected compounds calculated after 100 ns MD simulation.

Secondary structure analysis

Based on the conformational behavior of bound complexes noted in RMSD, RMSF, and Rg, we have calculated the secondary structure of 54035018-Mpro complex and apo form of Mpro enzyme as it is essential to observe the changes in enzyme structure without the inhibitor (Figure 8). As shown in Table 2, the average number of structural components are higher in the inhibitor bound Mpro complex as compared to the apo form of the Mpro enzyme. A slight increase was observed in the α-helix of the Mpro apoenzyme. The β-strands and 310-helix are marginally high in 54035018-Mpro complex indicative of a stable binding. However, no major change was observed in the secondary structure of the Mpro enzyme upon binding of the compound which shows strong stability, flexibility, and compactness of 54035018-Mpro complex. The secondary structural analysis characterized here might provide useful conformational information of Mpro which leads to the development of SARS-CoV-2 Mpro inhibitors. We presume that the design and development of selective inhibitors of Mpro enzyme using rational approaches may pave new routes in antiviral therapeutics.

Figure 8.

Secondary structure analysis; A: Compounds bound to Mpro enzyme; B: Apo form of Mpro enzyme, calculated after 100 ns MD simulation.

Table 2.

Percentage of residues contributing to the secondary structure in the Mpro enzyme.

	Motif Percentage of Protein Secondary Structure
System	a-Helix	b-Strands	310-Helix	Turn	Bend	Other
1-Mpro Complex	22	27	3	20	9	25
Mpro-Apo form	23	26	2	15	7	21

Mechanistic insights into binding affinity

The thermodynamic energy contribution of an inhibitor to the overall binding free energy of the complex corresponds to the structural stability of the inhibitor in the catalytic site of the enzyme. The intermolecular interactions in the catalytic/allosteric site residues participate considerably in the stability, binding affinity, and selectivity of the inhibitor. Thus, it was necessary to examine the binding affinity of all the selected compounds towards Mpro enzyme using MM/GBSA technique to establish the impact of the studied compounds. The results are presented in Table 3.

Table 3:

MM/GBSA-based binding energy profile of Mpro enzyme in complex with its inhibitors.

Complex	ΔEvdw	ΔEelec	ΔGgas	ΔGpolar	ΔGnopolar	ΔGsolvation	ΔGbind
54035018-Mpro	−44.79	−31.73	−43.83	28.63	−5.53	25.01	−37.40
54152887-Mpro	−39.09	−25.22	−56.21	30.76	−4.67	26.08	−37.18
54456426-Mpro	−32.12	−22.13	−54.25	33.57	−4.11	29.46	−24.79
57076946-Mpro	−35.44	−17.12	−67.17	34.25	−4.47	29.77	−18.82
91366909-Mpro	−27.24	−16.59	−70.02	38.37	−3.62	32.84	−24.79
N3-ILP-Mpro	−47.73	−22.65	−70.38	45.69	−6.19	39.49	−30.89

As estimated free binding energy (ΔGbind) of Mpro complex was calculated as the highest energy with an average value of −37.40 kcal/mol relative to 54152887-Mpro and 54456426-Mpro complexes with total mean values of −37.18 and −24.79 kcal/mol. The overall binding energy of control complex N3-ILP-Mpro appeared to be lower (-30.89 kcal/mol) than the 54035018-Mpro complex suggesting a better and improved binding of compound 54035018 to its target enzyme. We further estimated other components of the free binding energy associated with enzyme-inhibitor binding (Table 3).

We observed that the intermolecular van der Waals and electrostatic energies in complex 54035018-Mpro were favorable with average values of −44.79 and −31.73 kcal/mol whereas the gas-phase energy was higher (-56.21 kcal/mol) in complex 54152887-Mpro. There is a slight difference of 0.22 kcal/mol among the ΔGbind energies of 54035018-Mpro and 54152887-Mpro complexes. The ΔGsolvation contributed to the unfavorable binding of compound 54035018 as it is the least energy with a value of 25.01 kcal/mol among all complexes. Although the systematic movement of compound 54035018 and 54152887 from the solvent phase to the active site of Mpro stimulated van der Waals and electrostatic interactions with the catalytic residues, these interactions are not enough to contribute fully in the binding of Mpro enzyme as ΔGbind contributes to enhanced binding of these compounds. This analysis suggests that compound 54035018 bounds to Mpro enzyme with greater affinity in comparison with other complexes.

Conclusion

The necessity to control the emerging COVID-19 pandemic made us develop, a strategy to discover lead compounds that might be used in clinical trials to provide rapid success. Despite major efforts in the design and development of specific drugs or vaccines, not much proven to be efficacious against this SARS-CoV-2 infection. This curiosity leaves us to explore the drug designing approaches that could serve valuable to combat this disease. In this report, we have performed structure-based drug designing, virtual drug screening, and all-atom MD simulation approaches to discover highly selective compounds that possess significant binding affinity and presumably inhibition of the SARS-CoV-2 Mpro enzyme. We have identified the five best compounds, however, only three have been discussed in detail as they have shown favorable binding energy against Mpro enzyme. 54035018-Mpro complex exhibited most stable, flexible, and compact in all the complexes with the highest notable secondary structure elements and binding energy (ΔGbind −37.40 kcal/mol). Based on our overall observations, compound 54035018 could be recommended as a potential lead for the therapeutic management of COVID-19 patients after required clinical validations.

Click here for additional data file.