Computational study on peptidomimetic inhibitors against SARS-CoV-2 main protease.

The emergence outbreak caused by a novel severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) has received significant attention on the global risks. Due to itscrucial role in viral replication, the main protease 3CLpro is an important target for drug discovery and development to combat COVID-19. In this work, the structural and dynamic behaviors as well as binding efficiency of the four peptidomimetic inhibitors (N3, 11a, 13b, and 14b) recently co-crystalized with SARS-CoV-2 3CLpro were studied and compared using all-atom molecular dynamics (MD) simulations and solvated interaction energy-based binding free energy calculations. The per-residue decomposition free energy results suggested that the key residues involved in inhibitors binding were H41, M49, L141-C145, H163-E166, P168, and Q189-T190 in the domains I and II. The van der Waals interaction yielded the main energy contribution stabilizing all the focused inhibitors. Besides, their hydrogen bond formations with F140, G143, C145, H164, E166, and Q189 residues in the substrate-binding pocket were also essential for strengthening the molecular complexation. The predicted binding affinity of the four peptidomimetic inhibitors agreed with the reported experimental data, and the 13b showed the most efficient binding to SARS-CoV-2 3CLpro. From rational drug design strategies based on 13b, the polar moieties (e.g., benzamide) and the bulky N-terminal protecting groups (e.g., thiazole) should be introduced to P1' and P4 sites in order to enhance H-bonds and hydrophobic interactions, respectively. We hope that the obtained structural and energetic information could be beneficial for developing novel SARS-CoV-2 3CLpro inhibitors with higher inhibitory potency to combat COVID-19.

Introduction

The novel coronavirus (SARS-CoV-2) was first identified as the source of the illness with the infection trace in a seafood market in Wuhan city, Hubei province, China [1,2]. Nevertheless, infectious diseases of related coronaviruses at pandemic level occurred in 2002 and 2012 for SARS-CoV [[3], [4], [5]] and Middle East Respiratory Syndrome (MERS-CoV) [6,7], respectively. Coronaviruses are a large family of viruses that cause symptoms of the common cold and a serious respiratory illness that can lead to mortality [8]. More importantly, this current contagious disease is spreading significantly and there is still no vaccine or antiviral drug available.

SARS-CoVs are single stranded positive RNA viruses belonging to the Coronaviridae family. They consist of a 5′-untranslated region (UTR), a replicate complex (ORF1ab) encoding non-structural proteins (NSPs), a spike protein (S), an envelope protein (E), a membrane protein (M), a nucleocapsid protein (N), 3′-UTR, and several unidentified non-structural open reading frames [[9], [10], [11]]. In NSPs of ORF1ab, SARS-CoVs contain two proteases that cleave different sites of the replicase polyproteins. The 3 chymotrypsin-like cysteine protease (3CLpro), also called main protease (Mpro) cleaves 11 sites, whereas the papain-like protease (PLpro) cleaves three sites in the polyproteins [12,13]. Consequently, SARS-CoV 3CLpro plays a key role in proteolytic processing of the replicate polyproteins and has been recognized as a potential target for anti-viral drugs.

The sequence identity between SARS-CoV 3CLpro and SARS-CoV-2 3CLpro are 96%, and the amino acid residues involved in the catalytic and substrate-binding sites are 100% conserved [13,14]. Catalytic dyad is composed of the residues H41 and C145. Besides, Q189 residue commonly found in the main protease of both viruses involves in hydrogen bond formations with inhibitors [[14], [15], [16]]. The 3CLpro requires dimerization for its proteolytic activity [17]. The catalytic/substrate binding cleft of 3CLpro is located between the antiparallel β-barrels of domains I and II (Fig. 1 ) at the P1 and P2 positions [18], while the domain III is notably required for enzymatic reaction [19].

Fig. 1

Superimposition of SARS-CoV-2 3CLpro co-crystal structures. (A) Protomer A (blue) in complex with N3 (magenta stick, 6 LU7.PDB), 11a (red orange stick, 6LZE.PDB), and 13b (green stick, 6Y2F.PDB) and protomer B (pink) in unbound form. Domains I, II, and III are demonstrated. (B) A close-up view of the binding orientation of N3, 11a, and 13b inhibitors in the substrate-binding cleft. (C) The covalent bond (C–S) between peptidomimetic inhibitors and C145.

Broad-spectrum inhibitors targeting SARS-CoV 3CLpro have been designed and evaluated [[20], [21], [22]]. A recent discovery has shown that N3 is a potent Michael acceptor inhibitor to SARS-CoV-2 3CLpro (EC50 of 16.77 μM in SARS-CoV-2-infected Vero cells and k obs/[I] of 11,300 ± 880 M−1 s−1 for SARS-CoV-2 3CLpro) [23]. Meanwhile, L. Zhang et al. have developed a potent α-ketoamide (13b) inhibitor with the IC50 value of 0.67 ± 0.18 μM for anti-SARS-CoV-2 3CLpro and it exhibits a higher plasma half-life (~3-fold) and a lower rapid clearance than its parent lead compound [24]. A loss of the Boc group in compound 14b weakens the inhibitory potency and becomes almost inactive. Lately, the crystal structure of SARS-CoV-2 3CLpro in complex with the peptidomimetic aldehyde inhibitor 11a has been determined, indicating the achievement of a covalent bond formation between the aldehyde group and the catalytic residue C145 [25]. This compound shows excellent inhibitory potency with the IC50 value of 0.053 ± 0.005 μM and also possesses good pharmacokinetic properties with low toxicity. Accordingly, these different types of peptidomimetic inhibitors are promising for anti-viral drug developments against SARS-CoV-2. In this work, the molecular complexation between these four recently reported inhibitors (N3, 11a, 13b, and 14b; Fig. 2 ) and SARS-CoV-2 3CLpro was elucidated and compared by means of molecular dynamics (MD) simulations and binding free energy calculations.

Fig. 2

2D structure of the four studied compounds (Michael acceptor inhibitor N3, aldehyde inhibitor 11a, and two α-ketoamide inhibitors 13b and 14b). The canonical binding fragments are identified demonstrating binding pocket moieties (P5, P4, P3, P2, P1 and P1’). Atomic labels of each inhibitor are clarified for further discussion.

Computational details

Structural preparation

The co-crystal structures of SARS-CoV-2 3CLpro with the three inhibitors N3 (6LU7), 11a (6LZE), and 13b (6Y2F) were downloaded from the Research Collaboratory for Structural Bioinformatics (RCSB) Protein Data Bank [26]. All inhibitors were covalently bonded with C145 residue of SARS-CoV-2 3CLpro [[23], [24], [25]]. This covalent bond was removed using Accelrys Discovery Studio 2.5 [27]. The warhead of N3, 11a, and 13b at P1’ position was then restored to the α,β double bond, aldehyde group, and α-keto group, respectively, while the vacancy side chain of C145 in all structures was protonated. The two missing amino acid residues E47 and D48 in the 13b/3CLpro complex were built using SWISS-MODEL server [28]. The 13b/3CLpro complex was used to prepare the 14b system by eliminating the butyloxycarbonyl moiety at P4 site of 13b using the Accelrys Discovery Studio 2.5.

According to the standard protocols [[29], [30], [31]], the electrostatic potential (ESP) charges of ligands were calculated with the HF/6-31G(d) level using the Gaussian 09 program [32]. The antechamber and parmchk modules in AMBER program version 16 [33] were respectively used to generate the restrained ESP charges and prepare the parameters of all inhibitors. The bonded and non-bonded parameters for the protein and ligand were treated with the AMBER ff14SB force field [34] and generalized AMBER force field version 2 (GAFF2) [35], respectively. The protonation states of all ionizable amino residues were assigned at pH 7.4 using PDB2PQR server version 2.0.0 [36]. Hydrogen atoms were added using the LEaP module. The systems were firstly minimized to remove any steric hindrance or inappropriate geometry. Prepared structures were placed in a 10 Å truncated box filled with TIP3P water molecules [37]. Sodium ions were randomly added to neutralize the system. The added hydrogen atoms and water molecules were minimized using 500 steps of steepest descent (SD) followed by 1500 steps of conjugated gradient (CG) methods before running the MD simulations, while the remaining parts of the system were held fixed. The protein and ligand were further minimized by SD (500 iterations) and CG (1500 iterations) methods with constrained solvent molecules. Finally, the whole complex was fully minimized using the same procedure.

Molecular dynamics simulations

All molecular dynamics (MD) simulations were carried out using AMBER16 [33]. The prepared inhibitors/3CLpro complexes were heated from 10 to 310 K for 200 ps with harmonic positional restraint of 30.0 kcal/mol Å2 to Cα atoms of protein. Each complex was then subjected to four steps of restrained MD simulations at the same temperature with harmonic restraints of 30, 20, 10, and 5 kcal/mol Å2 for 1300 ps in total and another 200 ps was performed without any restraint. Finally, the entire system was simulated under the NPT ensemble (310 K, 1 atm) until reaching 100 ns. Non-bonded interactions were set to 10-Å cutoff distance, whereas the particle mesh Ewald summation method [38] was employed to treat the electrostatic interactions. The SHAKE algorithm [39] was used to constrain all covalent bonds involving hydrogen atoms. A 2-fs simulation time step was used throughout the MD simulations. Temperature and pressure were controlled by the Langevin thermostat [40] with a collision frequency of 2 ps−1 and the Berendsen barostat [41] with a pressure-relaxation time of 1 ps, respectively. The MD trajectories were saved every 10 ps. The results from MD simulations were analyzed using the same criteria as previously described [14]. The CPPTRAJ [42] was used to calculate the structural dynamics properties and intermolecular hydrogen bonding of the inhibitors/3CLpro complexes. The MM/GBSA per-residue decomposition free energy (ΔG bind residue) calculations [43] were carried out on the 100 snapshots extracted from the last 20-ns MD simulations using the MMPBSA.py [44] implemented in AMBER16. The binding free energy of protein–ligand binding (ΔG bind) was estimated from the solvated interaction energy (SIE) approach [45] as the summation of the van der Waals (E ), electrostatic (E ), reaction field (G ), cavity (γ∆SA(ρ)), and constant (C) as the following equation:where D is the solute interior dielectric constant and coefficients used in the calculations are α = 0.105, γ = 0.013, and C = −2.89. E and E are the intermolecular van der Waals and Coulomb interaction energies in the bound state, respectively. ΔG RF is the electrostatic contribution of the solvation free energy to binding, and ΔG cavity (γ∆SA) is the nonpolar contribution of the solvation free energy to binding.

Results and discussion

System stability

The system stability of all four inhibitors/SARS-CoV-2 3CLpro complexes was determined using all-atom RMSD upon the simulation time. Fig. 3 shows that the RMSD values of all systems are increased in the first 40 ns and then maintained at the fluctuation of ~2.0–3.0 Å until the end of simulation. Thus, the number of intermolecular H-bonds (# H-bonds) and the number of atom contacts (# atom contacts) of inhibitor plotted along the simulation were also taken into account in a consideration of equilibrium phase. It can be seen that all systems reache equilibrium at 40 ns; however, in this work the MD trajectories from 80 to 100 ns were extracted for further analysis. Over the last 20 ns, the # H-bonds of N3, 11a, 13b, and 14b are 3.33 ± 1.04, 2.97 ± 1.11, 3.15 ± 0.93, and 3.13 ± 0.98, respectively. The # atom contacts of these inhibitors are observed in the range of ~18–21, and among them, the α-ketoamide inhibitor 13b shows the highest # atom contacts. This result suggests that 13b could interact more favorable to the SARS-CoV-2 3CLpro than the remaining inhibitors (discussed in more detail later).

Fig. 3

All-atom RMSD, # H-bonds, and # atom contacts of N3, 11a, 13b, and 14b in complex with the SARS-CoV-2 3CLpro plotted along the 100-ns MD simulations.

Predicted inhibitory efficiency

The binding affinity of all studied inhibitors/3CLpro systems was estimated using solvated interaction energy (SIE) method on 100 snapshots extracted from the last 20 ns of simulations. The obtained results are listed in Table 1 , whereas the experimental binding free energies (ΔG Exp) were converted from the reported IC50 for 11a [25] and 13b [24] and K i for N3 [21]. The molecular mechanics energy (ΔE MM) shows that van der Waals (vdW) interaction plays a vital role in molecular complexation with SARS-CoV-2 3CLpro (ΔE vdW of −60.37 ± 0.34, −57.63 ± 0.55, −62.64 ± 0.36, and − 56.18 ± 0.28 kcal/mol for N3, 11a, 13b, and 14b, respectively), whereas the electrostatic attractions are observed in the range of ~ − 14 to −18 kcal/mol. It is noticed that vdW interaction is ~3.5–4.2-fold stronger than the electrostatic attraction, in line with several reported SARS-CoV-2 3CLpro inhibitors, including ritonavir and lopinavir [14], Cmp-17 [46], and α-ketoamide 13b [47]. From SIE-based ΔG bind calculations, we found that, among all studied compounds, the α-ketoamide inhibitor 13b shows the highest binding affinity (ΔG bind of −10.35 ± 0.04 kcal/mol), which is to some extent greater than α-ketoamide inhibitor 14b (−9.64 ± 0.03 kcal/mol) and the two remaining inhibitors (−9.92 ± 0.04, and − 9.68 ± 0.07 kcal/mol for N3 and 11a, respectively). The predicted ΔG bind of these inhibitors are somewhat in correspondence with the reported experimental data. However, it is worth noting that experimental values of IC50 and K i were obtained from different laboratories. Although all studied inhibitors are not covalently attached to SARS-CoV-2 3CLpro but the results show that these inhibitors are strongly bound in the complex form, which are comparable well with X-ray crystallographic data.

Table 1

Energy components (kcal/mol) of focused inhibitors/3CLpro complexes calculated with solvated interaction energy (SIE) method.

Energy component	N3	11a	13b	14b
SIE-based binding free energy
EvdW	−60.37 ± 0.34	−57.63 ± 0.55	−62.64 ± 0.36	−56.18 ± 0.28
Eele	−14.33 ± 0.30	−13.95 ± 0.24	−17.85 ± 0.28	−15.04 ± 0.31
ΔGRF	18.42 ± 0.21	15.90 ± 0.18	19.14 ± 0.20	15.66 ± 0.18
ΔGcavity	−10.81 ± 0.06	−9.10 ± 0.09	−9.85 ± 0.06	−8.93 ± 0.03
ΔGbind	−9.92 ± 0.04	−9.68 ± 0.07	−10.35 ± 0.04	−9.64 ± 0.03
Experimental binding free energy
ΔGExpa	−6.83* [21]	−9.92 [25]	−8.42 [24]	N/A

Data are shown as mean ± standard error of the mean (SEM). aΔGExp was converted from the reported IC50 values using the Cheng−Prusoff equation of ΔGExp = RT•ln(IC50) [48]. *ΔGExp of N3 was calculated from Ki.

Key residues upon inhibitors binding

The key residues involved in the binding of ligands to the SARS-CoV-2 3CLpro were investigated using ΔG bind residue calculation based on MM/GBSA method. The obtained results are shown in Fig. 4 , whereas the binding orientation of each inhibitor is depicted in the middle panel, in which the hot-spot residues are colored according to their ΔG bind residue values. It should be noted that only residues that exhibit the energy stabilization of ≤ −1.0 kcal/mol at least in one system are labeled and discussed.

Fig. 4

ΔGbindresidue of inhibitors/3CLpro complexes (left). The contributing residues involved in ligand binding are colored according to their ΔGbindresidue values, where the highest to lowest free energies is shaded from white to blue, respectively (middle). Binding pattern of inhibitors/3CLpro complexes is drawn from the last MD snapshot (right). Black dashed lines represent H-bond formation. H-bonds observed in the crystal structures and remained after MD simulations are denoted as “*”.

In consistent with the ΔG bind results (Table 1), the potent α-ketoamide inhibitor 13b shows the highest hot-spot residues (10 residues, including H41, M49, L141, N142, S144, C145, H164, M165, E166, and P168), whereas there are nine (L27, M49, G143, C145, M165, E166, P168, Q189, and T190) and eight (M49, F140, L141, C145, H163, H164, M165, and E166) important residues found in the Michael acceptor inhibitor N3 and the aldehyde inhibitor 11a, respectively. Note that these three potent compounds have the same binding residues as follows: M49, C145, M165, and E166. Although the eight residues (H41, M49, L141, G143, C145, H164, M165, and Q189) are found to be contributed for binding of the inactive inhibitor 14b, the interactions with the residues E166 and P168 positioned closed to the P4 site are dramatically reduced due to eradication of the Boc group.

Since the canonical fragments at P4 and P5 positions were designed to favor the S4 hydrophobic site [21], N3 is able to occupy the binding pocket in this region better than the other inhibitors and shows high affinity with Q189 (ΔG bind residue of −2.15 kcal/mol) and T190 (ΔG bind residue of −1.79 kcal/mol) residues. Notably, the hydrophobic contribution from M165 is found to be essential in all the studied inhibitors with the ΔG bind residue of −2.9, −3.1, −3.0, and − 2.0 kcal/mol for N3, 11a, 13b, and 14b, respectively. This finding is in accordance with previous studies demonstrating that M49 and M165 stabilized inhibitors through hydrophobic interactions [14,49]. The G143 and C145 residues interacted with 13b (ΔG bind residueof −0.95 and − 2.27 kcal/mol, respectively) and 14b (ΔG bind residueof −1.16 and − 2.58 kcal/mol, respectively) better than N3 and 11a inhibitors, since this oxyanion hole stabilization is commonly found in α-ketoamide inhibitors [22,24]. The high contribution from the residue P168 was detected in the inhibitors containing the P4 site systems (of −1.9 and − 1.2 kcal/mol for N3 and 13b), whereas this contribution is dramatically decreased in 11a and 14b. This result could be the reason why removing the Boc group (14b) weakens its inhibitory efficiency [24].

Superimposition over the last 20 ns of MD simulations (Fig. 5 ) shows that the active inhibitors N3, 11a, and 13b are well-aligned in the substrate-binding pocket, and among them, 13b is the most stable one. It should be noted that the S1 site is dramatically shifted to S1* for N3 and moves the benzyl ester portion up together (S1’ site) due to lacking of H-bond formations at the catalytic center residues. Additionally, 14b also has large fluctuations at S1’ and S3 sites, in good agreement with the reported low inhibitory activity toward SARS-CoV-2 CLpro [24].

Fig. 5

Superimposition of N3, 11a, 13b, and 14b inhibitors at the substrate-binding pocket of SARS-CoV-2 3CLpro derived from the last 20-ns MD simulations (magenta stick representation) compared with X-ray crystal structures (cyan stick representation). The canonical sites are labeled in blue text.

The electrostatic (ΔE ele + ΔG polar) and vdW (ΔE vdW + ΔG nonpolar) energy contributions of the residues important for inhibitor binding are depicted in Fig. 6 . We found that the main contribution for stabilizing the inhibitors is vdW interactions (ΔE vdW + ΔG nonpolar) up to −4.7 kcal/mol as seen by negative values (green grids), whereas there are rather low electrostatic attractions observed (~ − 1.0 to 2.6 kcal/mol shaded by light green to red grids), correlating well with the ΔE MM calculations based on SIE method (Table 1). The hydrophobic contributions from residues H41, M49, N142, C145, M165, and E166 are found to stabilize all the studied inhibitors. For contribution from H41, one of the catalytic dyad, the canonical moiety at P2 of all inhibitors stacks with the imidazole ring of H41 at the S2 site (Fig. 4). Interestingly, this catalytic residue shows high vdW contribution (< −2.4 kcal/mol) to the α-ketoamides (13b and 14b) over the remaining types of inhibitors, suggesting that cyclopropyl ring is preferentially occupied in this site, which commonly presents in preclinical/clinical drugs such as efavirenz, abacavir, and roflumilast [50]. It should be noted that the highest vdW contributions from the M165 (−3.11 kcal/mol) and E166 (−4.74 kcal/mol) residues are detected in the potent α-ketoamide inhibitor 13b, as previously reported [47]. However, 13b lacks the P5 site, resulting in a lower vdW contribution with P168 and Q189 residues (−1.31 and − 0.68 kcal/mol) compared to that of the Michael acceptor inhibitor N3 (−2.05 and − 3.49 kcal/mol).

Fig. 6

Electrostatic (ΔEele + ΔGpolar, red) and vdW (ΔEvdW + ΔGnonpolar, green) energy contributions from individual residues in the binding pocket of SARS-CoV-2 3CLpro toward the inhibitors N3, 11a, 13b, and 14b.

Protein-drug hydrogen bonding

In biological system, H-bonds play important role in many aspects, e.g., helping a complexation between ligand and protein. The percentages of the intermolecular H-bonds formed between the inhibitors and SARS-CoV-2 3CLpro residues are plotted in Fig. 7 , while the representative structures taken from the last snapshot of MD simulations are depicted in Fig. 4 (right panel). As expected, there are not many strong hydrogen bonds (> 80%) stabilizing the inhibitor binding, i.e., the E166 backbone nitrogen only interacts with the O1 atom of all three potent inhibitors at P3 site, while the N3 atom of the two α-ketoamides is highly stabilized by the H164 backbone oxygen.

Fig. 7

Percentage of H-bond occupation of the SARS-CoV-2 3CLpro residues contributed to N3, 11a, 13b, and 14b binding.

We found that there are three H-bonds formations detected in the Michael acceptor inhibitor N3 at P1 and P3 sites, including OE1 (Q189)···N-N4 (P1) at 37.3%, NE2 (Q189)···N-N4 (P1) at 13.7%, and O1 (P3)···H-N (E166) at 90.7%, which are similarly found in the initial structure from X-ray crystallographic data [23]. In the case of aldehyde inhibitor 11a, it mainly formed H-bonds with O1 (P3)···H-N (E166), O4 (P1’)···H-N (C145), and O (E166)···H-N1 (P3) at 88.7%, 58.4%, and 44.1%, respectively. For the α-ketoamide inhibitors, 13b and 14b share H-bond formation pattern as follows: O (H164)···H-N3 (P1) at 93.8% and 95.8%, O4 (P1’)···H-N (C145) at 35.3% and 57.7%, and O5 (P1’)···H-N (G143) at 70.1% and 55.8%. The stabilization of oxyanion hole from residues G143 and C145 via two H-bonds is commonly found in α-ketoamide inhibitors [22], whereas the warhead of the aldehyde [51] and Michael [52] inhibitors formed only one H-bond with the catalytic center of the target protease. It should be noted that the inactive 14b inhibitor has low level of H-bond formation (8.2%) with E166 compared to the other inhibitors.

Solvent accessibility of binding pocket

The ability of water accessibility at the active site of SARS-CoV-2 3CLpro with and without inhibitor bound (protomers A and B in Fig. 8A) was estimated by means of solvent-accessible surface area (SASA) calculation. The obtained results are shown in Fig. 8B, while the SASA averaged from the last 20 ns of MD simulations is listed in Table 2 . Note that the SASA calculations were performed on the residues located within 4 Å of the inhibitor. The obtained results showed that the SASAs for protomer B are 1289 ± 84, 1193 ± 76, 1087 ± 81, and 1033 ± 84 Å2 for the N3, 11a, 13b, and 14b systems, respectively. Once the inhibitor bound to the active site of protomer A, the water accessibility of this protomer is diminished to 734 ± 60, 879 ± 71, 919 ± 63, and 880 ± 60 Å2 for N3, 11a, 13b, and 14b, respectively. A reduction of water accessibility upon the binding of inhibitors to the enzyme active site is also found in the reported anti-HIV drugs, lopinavir and ritonavir, bound to the SARS-CoV-2 3CLpro [14] as well as in other protein-ligand complexations [[53], [54], [55]]. Since N3 is the largest molecule among studied inhibitors, accommodating well at the S1’ and S1-S5 sites, a lower water accessibility to the substrate-binding cleft is observed. Although the α-ketoamide inhibitor 13b shows the greatest binding affinity, the SASA value is not lower than the inactive compound 14b. This suggested that further structural modification on 13b should be made, e.g., enlarging the functional group at P4 site to better occupy the hydrophobic S4-S5 sites of SARS-CoV-2 3CLpro.

Fig. 8

(A) SARS-CoV-2 3CLpro homodimer, in which (A) each inhibitor is bound at protomer A binding pocket (blue) and protomer B is without an inhibitor (red pink). The surface of residues within 4 Å from an inhibitor of each protomer is shown. (B) SASA plots along the simulations time of N3, 11a, 13b, and 14b systems.

Table 2

Average SASA (Å2) during the last 20 ns MD simulations of N3, 11a, 13b, and 14b inhibitors. Data are shown as mean ± standard deviation.

System	Average SASA (Å2)
	Protomer A	Protomer B
N3	734 ± 60	1289 ± 84
11a	879 ± 71	1193 ± 76
13b	919 ± 63	1087 ± 81
14b	880 ± 60	1033 ± 84

Rational drug design

According to the aforementioned results, the 13b α-ketoamide inhibitor was selected as a template to propose newly designed protease inhibitors (Fig. 9 ) owing to the following reasons: (i) it has the highest binding efficiency to SARS-CoV-2 CLpro (Table 1), (ii) it shows the highest stability in the binding pocket (Fig. 5), and (iii) it exhibits the highest vdW contribution with the residues M165 and E166 (Fig. 6). However, 13b has lower interactions with S4-S5 sites, and its molecular complexation is mainly driven by vdW forces. Additionally, the phenyl ring of 13b at S1’ site is relatively unstable compared to the other moieties (Fig. 5).

Fig. 9

Rational drug design of the SARS-CoV-2 3CLpro inhibitors. (A) 2D structure of 13b with possible modified fragments. (B) Predicted ΔGbind of modified 13b compounds from SIE method. (C) ΔGbindresidue of D1–1/3CLpro complex. (D) The contributing residues involved in D1–1 compound binding are colored according to their ΔGbindresidue values, where the highest to lowest free energies is shaded from grey to blue, respectively. (E) Binding pattern of D1–1/3CLpro complexes drawn from the last MD snapshot. Black dashed lines represent H-bond formation.

To enhance the binding efficiency, especially electrostatic contributions, to SARS-CoV-2 CLpro, the structural modifications of P1’ and P4-P5 sites of 13b are suggested as follows: (i) the polar moieties (e.g. (1) benzamide and (2) benzoate fragments) should be added to the phenyl ring at P1’ region to facilitate H-bond formations with the surrounding residues T25, H41, and C44 and (ii) the bulky N-terminal protecting groups (e.g., (1) thiazole, (2) isoxazole, and (3) MPC moieties) should be introduced to the P4 position to enhance hydrophobic interactions with P168 and T190 residues. However, the cyclopropyl ring at P2 position should be conserved as it strongly stacks with the catalytic H41 via hydrophobic interaction [50] (Fig. 6).

Six new designed 13b derivatives were proposed and taken to perform MD simulations for 100 ns. The binding affinity of them toward 3CLpro was estimated using SIE method on 100 snapshots extracted from the last 20 ns of simulations. The obtained results are shown in Fig. 9B. There are four out of six compounds which show higher binding affinity than 13b (−10.35 ± 0.04 kcal/mol), including D1–1, D1–2, D2–2, and D2–3 with the ΔG bind of −11.61 ± 0.08, −10.72 ± 0.06, −11.28 ± 0.06, and − 10.73 ± 0.06 kcal/mol, respectively. Among these new design compounds, the modification of P1’ to benzamide and P4 to thiazole (D1–1) shows the highest binding efficiency. The ΔG bind residue of D1–1 is improved for N142, S144, M165, E166, and P168 compared to 13b. The binding pattern of D1–1/3CLpro complex is shown in Fig. 9E. In comparison with 13b, H-bond formation of O (H164)···H-N3 (P1) of D1–1 is mostly preserved (93.8%). Interestingly, the H-bond formation of benzamide fragment with C44 (H–N (P1’)···O (C44) is magnificently enhanced (69% occupation) as proposed.

Conclusions

In the present study, the key binding and susceptibility of the reported peptidomimetic inhibitors (N3, 11a, 13, and 14b) against SARS-CoV-2 3CLpro have been investigated using all-atom MD simulations. The ΔG bind residue calculations revealed that residues (i) H41, M49, L141, N142, S144, C145, H164, M165, E166, and P168, (ii) L27, M49, G143, C145, M165, E166, P168, Q189, and T190, and (iii) M49, F140, L141, C145, H163, H164, M165, and E166 are important for the binding of 13b, N3, and 11a, respectively. In the case of an inactive inhibitor 14b, residues H41, M49, L141, G143, C145, H164, M165, and Q189 are contributed. The molecular complexation between SARS-CoV-2 3CLpro and all the studied inhibitors is driven mainly by vdW interactions. Besides, hydrogen bond formations at F140, G143, C145, H164, E166, and Q189 residues were also important for their binding efficiency. The predicted ΔG bind somewhat agreed with the ΔG Exp, and among them, the 13b exhibited the most efficient binding to SARS-CoV-2 3CLpro. The bulky N-terminal protecting groups (e.g., thiazole, isoxazole, and MPC moieties) should be introduced at the P4 position in order to improve binding interactions with P168 and T190 hydrophobic residues, whereas the polar moieties (e.g. benzamide and benzoate fragments) should be added to the phenyl ring at P1’ site to facilitate H-bond formations with the surrounding residues T25, H41, and C44. This information can be helpful for further drug development of novel peptidomimetic inhibitors targeting SARS-CoV-2 3CLpro.

Declaration of Competing Interest

The authors declare no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.