Benchmarking the ability of novel compounds to inhibit SARS-CoV-2 main protease using steered molecular dynamics simulations.

BACKGROUND: The SARS-CoV-2 main protease (Mpro) is an attractive target in the COVID-19 drug development process. It catalyzes the polyprotein's translation from viral RNA and specifies a particular cleavage site. Due to the absence of identical cleavage specificity in human cell proteases, targeting Mpro with chemical compounds can obstruct the replication of the virus.
METHODS: To explore the potential binding mechanisms of 1,2,3-triazole scaffolds in comparison to co-crystallized inhibitors 11a and 11b towards Mpro, we herein utilized molecular dynamics and enhanced sampling simulation studies. RESULTS AND
CONCLUSION: All the 1,2,3-triazole scaffolds interacted with catalytic residues (Cys145 and His41) and binding pocket residues of Mpro involving Met165, Glu166, Ser144, Gln189, His163, and Met49. Furthermore, the adequate binding free energy and potential mean force of the topmost compound 3h was comparable to the experimental inhibitors 11a and 11b of Mpro. Overall, the current analysis could be beneficial in developing the SARS-CoV-2 Mpro potential inhibitors.

Introduction

Coronaviruses (CoVs) are considered as an alarming threat to humans. The 2019-nCoV is extremely homologous with SARS-CoV [1]. CoVs are viruses that carry a single positive-stranded enveloped RNA and induces a broad array of neurological, respiratory, and gastrointestinal disorders in the host [2,3]. However, SARS-CoV-2 is more transmissible than the other two coronaviruses. They are the parts of the Coronavirinae subfamily in the Nidovirales order of Coronaviridae group. From the findings of recent researches, this virus has a compact genomic arrangement, including various non-structural and structural proteins, such as proteases, spike glycoprotein, and nucleocapsid [4]. The CoVs Main protease (Mpro) is one of the vastly investigated targets to develop probable inhibitors [5]. The Mpro is also recognized as 3-Chymotrypsin-like protease (3CLpro) [6,7]. The Mpro is liable to process the polyproteins produced by the viral RNA transcription. The Mpro separates the polyproteins by its proteolytic action to create polypeptides (non-structural) [8]. These polypeptides are needed to generate four necessary structural proteins (Membrane, Spike-RBD, nucleocapsid, and Envelope) and other subordinate proteins [9,10]. The inhibition of Mpro would fundamentally block the expansion of the virus in the human body [6,8].

The Mpro is present in the form of a dimer (A and B protomers) in active biological conditions. Every protomer includes three domains. The domains I and II (8–101 and 102–184 residues) carry a catalytic dyad made by His41 and Cys145 residues. These parts contains 6 anti-parallel β-barrel stranded structures. Domain III (201–303) comprises 5 α-helices building a globular cluster (anti-parallel) [8]. The domain III controls the arrangement of the physiologically dynamic two monomeric conformation by furnishing salt bridge interchanges linking Arg4 and Glu290 of each protomer. The aperture of domains I and II resembles a substrate-binding site. The binding site comprises four (S1, S1’, S2, and S4) subsites formed by nearby residues [11]. The nonexistence of homologs of Mpro in humans fits it to be an appropriate target without cytotoxicity [6,12]. Several possibilities for COVID-19 therapy were considered from the starting of the pandemic. These include monoclonal antibodies, oligonucleotide-based treatments, peptides, interferon therapies, vaccines, drug re-purposing, and small molecules. To date, only Remdesivir is authorized for emergency use. The lack of specific FDA-approved medicines against Mpro of SARS-CoV-2 still exists. Several in-silico investigations have proposed prospective inhibitors from diverse source against distinct proteins of SARS-CoV-2 [5,[13], [14], [15], [16], [17]]. In earlier studies, the 1,2,3-triazole scaffolds were identified to possess antiviral, anti-inflammatory, antibacterial, and anti-cancerous properties [[18], [19], [20], [21]]. In addition, the 1,2,3-triazole scaffolds were also shown to have an antiviral effect on the Mpro in experimental studies [[22], [23], [24]]. Based on the above-information, we aimed to identify the Mpro inhibitors with high-affinity binding and stability from our selected molecules in comparison to co-crystallized inhibitors by employing molecular dynamics (MD) and enhanced sampling simulations.

Material and methods

Selection of target protein and compounds

The 3D co-crystal structure of the Mpro (PDB ID: 6M0K) with 1.50 Å resolution was obtained in PDB format with co-crystallized molecules 11a and 11b (IC50: 0.053 ± 0.005 and 0.040 ± 0.002 μM) [25]. The earlier synthesized seven novel molecules (3f, 3g, 3h, 3i, 3j, 5d, and 5e) out of eighteen 1,2,3-triazole scaffolds were chosen for the study (Fig. S1) [26]. The structure of protein was prepared by the “prepare protein” tool of the Discovery studio [27]. Gaussian16 DFT (energy minimization) protocols were used to optimize the ligand geometry [28].

Molecular docking study

The molecular docking was done using the (CDOCKER method) BIOVIA Discovery Studio (DS) 2018 [29]. We assigned the binding pocket of co-crystallized molecule 11b for the docking of 1,2,3-triazole scaffolds. After selecting 11b, the SBD-site sphere showed the attributes X: 11.622; Y:11.881; Z:68.704 with a 12.00 Å radius for sphere object covering the binding site. The 6M0K was assigned as input receptor in parameter values setup and a set of chosen molecules defined as input ligands. The orientations to refine, random conformations, parallel processing, and simulated annealing parameters were specified as default. After completing calculations, the top hit complexes were selected on the basis of CDOCKER interaction energy and visualized the protein-ligand interactions. All the complexes were then prepared for MD simulations to examine the stability of the protein-ligand complexes.

Molecular dynamics simulations study

We selected docked 1,2,3-triazole scaffolds-Mpro complexes along with 11a and 11b (standard Mpro inhibitors) for 100 ns MD simulation applying GROMOS96 43a1 force field to validate the docking outcomes and estimate their stability [[30], [31], [32]]. The topology files for ligands were acquired from the PRODRG2 server [33]. We selected the simple-point charge water model to solvate the systems [34]. We created a 954.61 nm3 volume cubic box for all the designated systems. After running the gmx genion command, 4 Na + ions were added, and the energy minimization was accomplished by employing 50,000 steps of the steepest descent. The simulations were executed following the periodic boundary conditions with NVT accompanied by an NPT ensemble. Berendsen's coupling algorithm was carried throughout the restraint MD runs to keep the 300K temperature and 1 bar pressure steady, respectively [35,36] The electrostatic associations were specified by the PME algorithm, with a 1.2 nm C cutoff [37]. The LINCS algorithm was designated to constraint the bond lengths [38].

Binding free energy study

The molecular mechanics Poisson-Boltzmann surface area (MM-PBSA) approach was utilized to evaluate all the compound's binding free energies with Mpro to determine complexes stability [39]. The g_mmpbsa package was employed to assess the binding free energies of 1,2,3-triazole scaffolds and experimental inhibitor complexes. Evaluation of binding free energy implicates enumeration of non-polar solvation, van der Waals, solvent-accessible surface area, electrostatic, and polar energies.

Entropy calculation by Quasiharmonic method

Schlitter proposes the Quasiharmonic (QH) technique to estimate the relative and absolute entropies. These entropies were calculated on the basis of the assessment of the covariance matrix of Cartesian coordinates by MD simulations. Schlitter's approach directs only a determinant's computation, which was the premise for its endorsement in computational biology. In Schlitter's framework, the absolute entropy calculation describes the upper limit of the quantum mechanics entropy [40].

SMD and umbrella sampling simulations

The time-dependent external pulling force is applied in steered MD simulation to separate the compound from the protein's binding site [41]. The steered MD simulations were performed using the GROMACS software v4.6.7 [42]. The energy of all the identified protein-ligand structures was minimized before starting the pulling process. We applied the external pulling force on the z-axis along the simulated pathway to discharge the compound from the protein's binding site. The simulation for 500 ns was executed at the 0.01 nm/ps pull rate and allocated a 500 kJ/mol/nm2 spring constant. We constructed the force profiles of each system, determining the steered MD trajectories. Moreover, we assembled the distance summary to produce umbrella sampling slots [43]. Slot spacing of 0.1 nm is used to initiate the center of mass (COM) distance of 3 nm. After that, we generated sampling windows at each spacing distance of 0.2 nm. A 10 ns simulation run was allotted for every sampling window, finishing in 1170 ns of simulation period for each chosen protein-ligand complexes. Then the trajectories were constrained to the WHAM computation procedure to construct potential mean force (PMF) graphs. The PMF graph is utilized to compute the binding free energies.

Results and discussion

Analysis of binding poses

Molecular docking is a method of structure-based drug design to determine the necessary residues involved in protein-ligand interactions, as well as the binding affinity between protein and ligand [13,44]. The 11a, 11b (co-crystallized molecule), and 1,2,3-triazole scaffolds received after geometry and energy minimization (ligand preparation) were docked at the binding pocket of Mpro. We estimated the CDOCKER interaction energy of all the docked compounds, as represented in Table 1 . It has to be stated that a lower (more negative) interaction energy score signifies more favorable binding of protein-ligand [45]. The binding interaction of top protein-ligand poses were represented in Fig. 1 and S2. The examination of binding poses disclosed that all the docked compounds moderately interacted with the Mpro binding site, thereby showing their feasible inhibitory tendencies against Mpro of SARS-CoV-2. After perceiving the interactions of docked compounds with Mpro, it was noticed that the compounds formed interactions with both catalytic site residues Cys145 and His41. The 11a-Mpro complex was stabilized via six hydrogen bonds, including Glu166, His163, Ser144, Asn142, Gln189, Asp187 residues, and four hydrophobic interactions with Met165, Pro168, Met49, and His41 residues. Leu167, Gly143, Cys145, Phe140, His172, Leu141, His164, Val186, Arg188, and Gln192 residues formed van der Waal interactions. Compound 11b interacted with Mpro via hydrogen bonds with six residues (Glu166, His163, Ser144, Asp187, Asn142, Gln189) and formed hydrophobic interactions with His41, Met49, Met165, and Pro168 residues. Moreover, residues Leu167, Gly143, Cys145, Phe140, His172, Leu141, His164, Val186, Arg188, and Gln192 formed van der Waal interactions.

Table 1

CDOCKER Interaction Energy data of co-crystallized and selected compounds with SARS-CoV-2 Mpro.

S. No.	Complexes	CDOCKER Interaction Energy (kcal/mol)
1.	Mpro-11a	−47.750
2.	Mpro-11b	−50.934
3.	Mpro-3f	−38.016
4.	Mpro-3g	−36.571
5.	Mpro-3h	−37.091
6.	Mpro-3j	−38.550
7.	Mpro-5d	−35.620
8.	Mpro-5e	−34.781
9.	Mpro-3i	−41.041

Fig. 1

3-D binding poses showing SARS-CoV-2 Mpro complexes with co-crystallized and 1,2,3-triazole scaffolds.

The 3f-Mpro complex achieved stability by making hydrogen bonds with Asn142, Glu166, and Phe140 residues. The residues Cys145, Met49, Met165, and Leu141 contributed to the hydrophobic interactions. The residues His172, Ser144, His163, His164, Asp187, Val186, Gln192, Arg188, Gln189, Ser46, and His41 were involved in van der Waal interactions. The 3g-Mpro complex acquired stability via hydrogen bonds contributed by residues Asn142, Phe140, and Glu166. The hydrophobic interactions were made with Met165, Met49, Cys145, Leu141, and Glu166 residues. Ser144, His163, His164, Asp187, Val186, Arg188, Gln189, His41, and Ser46 showed van der Waal interactions. In the 3h-Mpro complex, the hydrogen bonds created by the residues Leu141, Asn142, Phe140, and Glu166 stabilized the binding site. The hydrophobic interactions were formed by Cys145, Met49, Met165, and Glu166 residues. The residues Ser144, His163, His164, Val186, Asp187, His41, Arg188, and Ser46 showed van der Waal interactions. The 3j-Mpro complex was stabilized via three hydrogen bonds (Asp187, His164, His41 residues), and hydrophobic interactions with Met49, Met165, Gln189, and His41 residues. The residues Val186, Arg188, Cys145, Glu166, Leu167, and Asn142 formed van der Waal interactions. Compound 5d interacted with Mpro via hydrogen bonds with Leu141, Asn142, and hydrophobic interactions with Glu166, Met165, and Met49 residues. However, residues His172, His164, Phe140, Asp187, Val186, Arg188, Gln189, His41, and Ser46 exhibited van der Waal interactions. The 5e-Mpro complex attained stability via hydrophobic interactions bestowed by His163, Met49, and Met165 residues. The van der Waal interactions were shown by Ser144, Gly143, Phe140, Leu141, Asn142, His164, His41, Gln192, Val186, Asp187, Arg188, Gln189, Glu166, Cys145, and His172 residues. Compound 3i interacted with Mpro via five hydrogen bonds with residues Ser144, His163, Asn142, His41, Gly193, and Phe140. The residues participating in hydrophobic interactions were Glu166, Met165, Cys145, and Met49. Moreover, residues His172, Leu27, Thr26, Thr25, Ala191, Thr190, Arg188, Gln189, Leu141, Gln192, and Leu167 exhibited van der Waal interactions. Overall, all the docked compounds showed hydrophobic interactions and hydrogen bonds with catalytic residues His41 and Cys145. Most of the 1,2,3-triazole scaffolds have displayed equivalent binding interactions in the binding pocket of Mpro as compared to the co-crystallized compounds 11a and 11b (see Fig. 1).

Structural dynamics of protein-ligand complexes

As the protein-ligand interactions obtained from molecular docking are static entities, thus performing MD simulations are deemed to be a crucial part of any in-silico analysis [[46], [47], [48]]. It provides detailed data concerning to the protein-ligand interactions with a dynamic aspect [49,50]. MD simulation was implemented to obtain more profound insights into the impact of protein flexibility and structural alteration in the interaction profile of the complexes. The system is stabilized and equilibrated if it acquires lower RMSD with constant variations in the complete simulation [51]. RMSD was measured for C-α atoms of protein backbone for all eight complexes. The RMSD graphs of Mpro complexed with 11a, 11b, 3f, 3g, 3h, 3i, 3j, 5d, and 5e were displayed in Fig. 2 . We saw that the maximum RMSD values achieved by all the Mpro complexes were very low, ranging from 0.2 nm to 0.45 nm. In contrast, the average RMSD for apo-protein was ∼0.34 nm (Fig. S3). Such low RMSD values strongly indicated conformational stability that the protein had obtained in the presence of compounds. All the complexes represented stable trajectories with steady and lesser fluctuations, signifying that the protein backbone experienced minor structural disturbance.

Fig. 2

The Root Mean Square Deviation (RMSD) graph of the Mpro C α-atoms in complex with 11a (black), 11b (blue), 3f (green), 3g (red), 3h (dark green), 3j (cyan), 5d (magenta), 5e (violet), and 3i (yellow).

RMSF, the measure of flexibility of Cα atoms of a protein is a crucial parameter in determining the stability of a protein-ligand complex. The residues play an indispensable role in attaining a stable conformation for the strong binding of a protein-ligand complex [52]. RMSF for all the Mpro complexes were measured and displayed on a graph, shown in Fig. S4. The fluctuations of binding site residues Glu166, Cys145, Asn142, Phe140, Met165, Met49, His41, and Leu141 were minor in all the Mpro complexes. Small fluctuations of these binding site residues recommended high stability of Mpro complexes. Additionally, the data intimated that the residues involved in the binding were also liable for the stable RMSD of the Mpro complexes.

Furthermore, the radius of gyration (Rg) is applied to determine the compactness of the protein-ligand complexes [53]. The Rg can estimate the unfolding and folding of a protein upon ligand binding. The high Rg value bestows unfolded (less compactness) protein-ligand complex. The average Rg values of Mpro complexes with selected compounds were between 2.09 and 2.13 nm, while between 2.13 and 2.10 nm for co-crystallized compounds (Fig. 3 ). As a result, every complex presented an approximately similar expression of compactness and uniform Rg values compared to the co-crystallized compounds. This implied that the selected compounds helped to sustain the complex's stability and compactness, similar to co-crystallized compounds.

Fig. 3

Radius of gyration (Rg) plot depicting the changes observed in the conformational behavior of Mpro complexes 11a (black), 11b (blue), 3f (green), 3g (red), 3h (dark green), 3j (cyan), 5d (magenta), 5e (violet), and 3i (yellow).

Additionally, we also executed ensemble cluster analysis to elucidate the protein's considerable flexibility towards the compounds and the degree of stability as a complex [54]. We conducted an ensemble cluster investigation over the MD trajectories of all the Mpro complexes, as shown in Fig. 4 . This interpretation furnished a substantial number of clusters (ranging from 46 to 57) for Mpro complexes with all the selected compounds. The average RMSD values ranging from 0.169 to 0.191 nm and energy of matrix between 3.00 and 3.75 nm were admitted for the complex's conformational stability. The less number of members in a cluster indicated that the mean conformation of the complexes in three-dimensional space was compact, with substantially less variance between clusters, referring that the protein-ligand complexes were conformationally stable.

Fig. 4

Pictorial representation of conformational flexibility using ensemble cluster analysis of Mpro in complex with (a) 11a, (b) 11b, (c) 3f, (d) 3g, (e) 3h, (f) 3j, (g) 5d, (h) 5e, and (i) 3i.

Binding free energy and conformational stability

The energy of interaction indicates the strength of the protein-ligand complex. Therefore, to validate the interaction energy of molecular docking analysis, a thorough investigation was conducted to estimate the binding free energy using the MM-PBSA (Table 2 ) approach. The data indicated that the selected compounds showed comparable binding free energies with Mpro as shown by co-crystallized compounds (Fig. S5). These results inferred that all the selected compounds made a strong binding with the Mpro. The van der Waals, SASA, and electrostatic energy furnished favorable energy to the total binding free energy for all the protein-ligand complexes. In contrast, the polar solvation energy demonstrated an unfavorable contribution to total binding free energy. The MM-PBSA results recommended that the predicted energies of all the selected compounds were thermodynamically favorable. In addition, we have also explored the per residue contribution energy to analyze the residues responsible for the difference in binding free energies between 11a, 11b, and 3h (Table S1). Furthermore, the configurational entropy was calculated for the chosen protein-ligand complexes with molecules 11a, 11b, and 3h (topmost) by a QH approach. The entropic effects are vital in molecular association, protein folding, stability, and ligand binding [55]. The entropy due to the QH calculation unveiled that all the selected protein-ligand complexes had comparable entropic values (Table S2). These outcomes showed that the complex with 3h molecule acquired similar entropy and ordered conformation as shown by 11a and 11b complexes. Moreover, the free energy landscape (FEL) for the PC1 and PC2 (two principal components) were generated utilizing projected eigenvectors values of PCA to verify the binding impacts of the compound at the conformational distribution's. The FEL plots demonstrating different patterns and stable conformations for the selected compounds 3h (top most), 11a, and 11b were shown in Fig. 5 . The 11a-Mpro complex had single minima with a flattened end, and 11b had two minima basins with one canonical and one flat end. In 3h bound complex, the basins were segregated into three fragments with two deep canonical ends and one flat end showing the lowest energy minima. The energy minima with a canonical end suggested a stable conformation, while the flat end indicated the lack of minima-energy conformations. The comparison of FELs for all the selected complexes revealed that the 3h-Mpro complex exhibited a comparable and stable free energy surface to 11a and 11b Mpro complexes.

Table 2

Binding free energy values of selected Mpro complexes using MM-PBSA.

Mpro Complexes with	ΔE binding (kcal/mol)	ΔE Van der Waal (kcal/mol)	SASA (kcal/mol)	ΔE Electrostatic (kcal/mol)	ΔE polar solvation (kcal/mol)
11a	−53.72	−64.13	−4.77	−9.83	25.02
11b	−53.29	−65.52	−5.10	−9.92	27.24
3f	−42.09	−69.88	−5.06	−28.28	61.13
3g	−35.04	−59.40	−4.47	−37.39	66.22
3h	−53.69	−67.00	−4.73	−6.15	24.19
3j	−38.98	−54.58	−4.52	−13.29	33.41
5d	−24.48	−40.88	−3.15	−35.20	54.76
5e	−18.80	−49.73	−4.14	−40.95	76.03
3i	−40.57	−66.59	−5.15	−43.06	74.23

Fig. 5

Free energy landscape of Mpro as a function of first two principal components PC1 and PC2 (a) 11a, (b) 11b, and (c) 3h.

Umbrella sampling simulations

Umbrella sampling simulation is a prominent technique for computing the potential mean force (PMF) to examine protein-ligand binding/unbinding procedures [56,57]. The reaction coordinate (ξ) along the free-energy profile was examined using the “gmx wham.” The received PMF value of every studied system is illustrated as a reaction coordinates (ξ) function. We calculated the PMF required for the unbinding of the topmost molecule (3h) obtained from MM-PBSA analysis compared to co-crystallized inhibitors 11a and 11b of Mpro. The plotted binding free energy can be outlined as the variation between the highest and least PMF curves, as represented in Fig. 6 . The co-crystallized inhibitors 11a and 11b exhibited binding free energy values of −64.70 ± 2.15 kJ/mol and −79.97 ± 1.11 kJ/mol, respectively, while the 3h delivered convenient binding free energy (−78.91 ± 2.26 kJ/mol) value. The free energy initiates, fall to the lowest value, and eventually extends to a stable value when ξ acquires 3.0–3.5 nm (Fig. 6). This span resembles the condition where the protein-ligand interaction is lost, and the ligand detached from the binding cavity of Mpro along the unbinding pathway. These free energy profiles complement our MM-PBSA calculation, where binding free energy of 3h was comparable to the co-crystallized molecules 11a and 11b. In the end, the selected 3h compound could be considered as lead to develop the Mpro inhibitor.

Fig. 6

Potential of mean force (PMF) curves obtained after umbrella sampling for the co-crystallized and the best selected compound in complex with Mpro SARS-CoV-2.

Conclusion

The urgency to prevent COVID-19 has suggested discovering potential lead molecules that can be utilized to target the SARS-CoV-2 Mpro. This research has implemented molecular docking, MD simulation, MM-PBSA, and enhanced sampling methods toward Mpro to recognize the prospective compounds with higher affinity and stability compared to co-crystallized compounds. MD simulation studies validated the docking outcomes and suggested that the selected compounds have the ability to bind and block the SARS-CoV-2 Mpro similarly as co-crystallized compounds. In addition, the umbrella sampling analysis also reported that the force required for 3h to leave the binding site of the Mpro was similar and comparable to the experimental inhibitors 11a and 11b. Moreover, the 3h compound concluded from this study will be helpful for designing new SARS-CoV-2 inhibitors.

Consent for publication

All the authors have read and approved the manuscript in all respects for publication.

Declaration of competing interest

There are no competing interests declared by the authors.