In Silico Elucidation of Potent Inhibitors and Rational Drug Design against SARS-CoV-2 Papain-like Protease.

Global public health has been a critical problem by the sudden increase of the COVID-19 outbreak. The papain-like protease (PLpro) of SARS-CoV-2 is a key promising target for antiviral drug development since it plays a pivotal role in viral replication and innate immunity. Here, we employed the all-atom molecular dynamics (MD) simulations and binding free energy calculations based on MM-PB(GB)SA and SIE methods to elucidate and compare the binding behaviors of five inhibitors derived from peptidomimetic inhibitors (VIR250 and VIR251) and naphthalene-based inhibitors (GRL-0617, compound 3, and compound Y96) against SARS-CoV-2 PLpro. The obtained results revealed that all inhibitors interacting within the PLpro active site are mostly driven by vdW interactions, and the hydrogen bond formation in residues G163 and G271 with peptidomimetics and the Q269 residue with naphthalene-based inhibitors was essential for stabilizing the protein-ligand complexes. Among the five studied inhibitors, VIR250 exhibited the most binding efficiency with SARS-CoV-2 PLpro, and thus, it was chosen for the rational drug design. Based on the computationally designed ligand-protein complexes, the replacement of aromatic rings including heteroatoms (e.g., thiazolopyridine) at the P2 and P4 sites could help to improve the inhibitor-binding efficiency. Furthermore, the hydrophobic interactions with residues at P1-P3 sites can be increased by enlarging the nonpolar moieties (e.g., ethene) at the N-terminal of VIR250. We expect that the structural data obtained will contribute to the development of new PLpro inhibitors with more inhibitory potency for COVID-19.

Introduction

The coronavirus disease 2019 or COVID-19 epidemic, which is caused by the severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) in January 2020, has already become a global problem.[1,2] Fever, cough, and shortness of breath are common signs of SARS-CoV-2 infection found in patients. Severe cases may develop pneumonia, SARS, organ failure, and probability of death.[3,4] SARS-CoV-2 has a long incubation time and a rapid transmission speed, more than other coronaviruses as SARS-CoV and Middle East respiratory disease (MERS-CoV). Furthermore, those who are asymptomatic can transmit the virus to others.[5,6] Nowadays, the number of SARS-CoV-2-infected people has been increasing. COVID-19 had spread to more than 222 nations as of February 03, 2021, with more than 104.9 million infections and more than 2.27 million deaths.[7] Currently, the number of COVID-19 patients throughout the world is still high. As a result, there is a significant need for effective therapies and medications.

SARS-CoV-2 belonging to Betacoronavirus shows a genome identical to a bat coronavirus with 96% similarity and shares 79.6% sequence identity to SARS-CoV.[7] This virus contains a large positive-stranded ssRNA genome encapsulated within a membrane envelope.[8] The CoV-2 spike glycoprotein regulates the virus entry into host cells. Subsequently, the CoV-2 genomic sense ssRNA replication mechanism is used as mRNA to produce 16 nonstructural proteins (Nsps) from two polyproteins, Pp1a and Pp1ab.[9] Proteolytic processing is achieved by two viral cysteine proteases, which are encoded by the SARS-CoV-2 genome including chymotrypsin-like main protease (3CLpro, Nsp5 domain) and papain-like protease (PLpro, Nsp3 domain). 3CLpro cleaves 11 sites of viral polyproteins with sequence consensus X-(L/F/M)-Q↓(G/A/S)-X.[10] Many SARS-CoV-2 main protease interactions with drug/inhibitors mechanisms have been reported.[11−13] The second protease, less examined compared to 3CLpro, is PLpro, which cleaves three sites with recognition sequence LXGG↓ (X = N or K for SARS-CoV-2).[14−16] PLpro is another attractive target because it plays an essential role in not only the cleavage of viral polyproteins for the proliferation cycle but also the disruption of host responses.[14,17] Note that only 54 residues of SARS-CoV-2 PLpro are different from SARS-CoV PLpro, while they share similar residues at the active site.[18] The catalytic triad residues C111, H272, and D286 are located at the interface between the thumb and palm subdomains, β14−β15 loop (P1 site), whereas the palm subdomains (P2, P3, and P4 sites) are substrate-binding clefts of PLpro (Figure A).[19] Furthermore, the BL2 loop residues Y268 and Q269 in the palm subdomain are involved in the formation of hydrogen bonds with inhibitors.[14,20]

Figure 1

(A) 3D structure of SARS-CoV-2 PLpro bound with the focused inhibitors: VIR250, VIR251, GRL-0617, compound 3, and compound Y96, where the catalytic triads are shown in the inset. (B) Atomic labels of each inhibitor explained for further discussion.

The SARS-CoV-2 PLpro inhibitors have been studied as follows. The peptidomimetic inhibitors Ac-Abu(Bth)-Dap-Gly-Gly-VME (VIR250, IC50 of 1.60 μM) and Ac-hTyr-Dap-Gly-Gly-VME (VIR251, IC50 of 1.65 μM) can inhibit PLpro activity via covalent bond formation between their vinylmethyl ester (VME) and the catalytic residue C111 at the P1 site.[21] By considering the catalytic efficiency (Kcat/Km), VIR250 (600 M–1 s–1) is more selective against PLpro than VIR251 (1000 M–1 s–1). Meanwhile, the naphthalene-based reversible inhibitors had the peptidase activity of PLpro (IC50 of 0.6–2.3 μM for GRL-0617 and 6.4 μM for compound 3).[14,16] GRL-0617 also showed potency for inhibiting SARS-CoV-2 replication in Vero E6 cell infection (EC50 of 1.4 μM). The crystal structures of SARS-CoV-2 PLpro and naphthalene-based hydroxyamino-benzamide derivative (compound Y96, PDB code: 7KOL) are released; however, the inhibition activity of this compound has not been reported yet. In this work, we aimed to elucidate an insight into the binding pattern and compare the binding affinity of these five inhibitors VIR250, VIR251, GRL-0617, compound 3, and compound Y96 (Figure B) against SARS-CoV-2 PLpro using molecular dynamics (MD) simulations and free energy calculations. We expect that the obtained structural information could be used as introductory guidance for further development for novel anti-PLpro drug candidates for SARS-CoV-2 treatment.

Computational Details

System Preparations

The crystal structures of SARS-CoV-2 PLpro in the apo form (PDB code: 6W9C) and PLpro in complex with five different inhibitors VIR250 (6WUU), VIR251 (6WX4), GRL-0617 (7JRN), compound 3 (7JIW), and compound Y96 (7KOL) were retrieved from the Protein Data Bank.[22] The missing residues (225–227) of the apo form at the finger lobe were built using the SWISSMODEL server.[23] The protonation states of all ionizable amino acids and inhibitors were predicted at pH 7.4 using the PDB 2PQR server[24] and ChemAxon software,[25] respectively. The sulfhydryl group of cysteine residues 189 and 224, which were coordinated with Zn2+, was deprotonated (CYM type, charge −1), while the other cysteines were neutral. The histidine residues 73, 89, and 175 were set as the neutral form with a protonated N-delta position (HID type), while histidine residues 17, 47, 50, 255, 272, and 275 were set as the neutral form with a protonated position at N-epsilon (HIE type). Note that the catalytic residue H272 was set as HIE type according to the reaction mechanism of papain-like proteases.[26] For inhibitors, the total charge of VIR250 and VIR250 was +1 at the amino group, while that of GRL-0617, compound 3, and compound Y96 was 0 (Figure ). All of the ligands were optimized with the HF/6-31g(d) level of theory using the Gaussian 09 program,[27,28] and then, the electrostatic potential (ESP) charges and restrained ESP (RESP) charges were generated.[29] The force fields of FF14SB[30] and generalized AMBER force field version 2 (GAFF2)[31] were applied for proteins and inhibitors, respectively. The tleap module was used to add all missing hydrogen atoms. Then, any steric hindrance or improper geometry involved to hydrogen atoms was removed by minimization using 1000 iterations of steepest descent (SD) followed by 3000 iterations of conjugated gradient (CG). The system was soaked in a 10 Å octahedral box of explicit waters using the TIP3P model and was subsequently neutralized by adding chloride ions. The water molecules were minimized using 1500 steps of SD followed by 3000 steps of CG, while the remaining parts of the system were restraint using a force constant of 500 kcal/mol·Å2. Subsequently, using the same protocol, the whole complex was fully minimized without any restraint before running the MD simulation.

Molecular Dynamics Simulations

SARS-CoV-2 PLpro in apo and complex forms was simulated by all-atom MD simulations with three different initial velocities using AMBER16[29] under periodic boundary conditions. Nonbonded interactions were considered using the short-range cutoff of 12 Å, whereas long-range electrostatic interactions were treated using Ewald’s method.[32] Berendsen algorithm was used to regulate the pressure.[33] All covalent bonds containing hydrogen atoms were constrained using the SHAKE algorithm.[34] The simulated models were heated to 310 K with 100 ps for relaxation time. The time step was set as 2 fs.[35−38] A Langevin thermostat with a collision frequency of 2.0 ps was used to control the temperature. Finally, a 100 ns unrestrained NPT simulation at 310 K and 1 atom was carried out. The temperature versus simulation time plot is shown in Figure S1. The MD trajectories were collected every 500 steps for later analysis. Protein backbone atoms (CA, C, O, and N) without water molecules and ions were used for RMSD calculations. The dynamic properties, intermolecular hydrogen bond occupation, the number of contacts, the distance between the inhibitor and active site residues, and the solvent-accessible surface area (SASA) were evaluated by the CPPTRAJ module.[39] Noted that the distance and the angle between the hydrogen bond donor (HD) and hydrogen acceptor (HA) of ≤3.5 Å and ≥120°, respectively, were used as a criterion for hydrogen bond calculations.

Binding Energy Calculations

To calculate the free energies of protein–ligand complexes, the molecular mechanics/Poisson–Boltzmann surface area (MM/PBSA) and molecular mechanics/generalized Born surface area (MM/GBSA) methods[40] were performed using the 100 snapshots extracted from the last 20 ns MD simulations of each system to ensure that all of the simulated systems actually reached the equilibrium state. On top of that, as shown in the snapshot per time from the last 20 ns of simulations (Figure ), all studied inhibitors were found to be stable at the active site of SARS-CoV-2 PLpro. A detailed description used for the MM-PB(GB)SA calculation is as follows: the internal dielectric constant was set to 1, while the external dielectric constant was 80. The grid size, the surface tension, and the solvent probe radius were set to 0.5 Å, 0.0072 kcal/mol·Å2, and 0.14 Å, respectively. It should be noted that MM-PB(GB)SA calculations are dielectric constant-sensitive.[41] The binding free energy (ΔGbind) of MM-PB(GB)SA between a ligand (L) and a receptor (R) to form a complex RL is calculated as follows[40,42]The changes in the gas-phase MM energy, the solvation free energy, and the conformational entropy are represented by ΔEMM, ΔGsol, and −TΔS, respectively. The ΔEMM includes ΔEinternal, which is calculated from the bond, angle, and dihedral energies of the atom, and ΔEelectrostatic and ΔEvdW are electrostatic and van der Waals energies, respectively. The solvation free energy (ΔGsol) is the sum of the polar and nonpolar energy terms. The polar contribution is determined using either the GB or PB model (ΔGPB/GB), while the nonpolar energy is estimated using the solvent-accessible surface area (SASA, ΔGSA). Normal-mode analysis on a set of conformational snapshots collected from MD simulations is used to compute the conformational entropy change −TΔS.

Figure 4

Superimposition of inhibitor(s) between the 20 snapshots derived from the last 20 ns of simulations compared to the X-ray structure represented in the green stick model within SARS-CoV-2 PLpro. The representative frame of the most frequent binding of each complex from the run1 system is given in Files S1–S5.

In addition, the solvated interaction energy (SIE) approach was utilized to predict the ΔGbind using the same sets of MD snapshots from MM-PB(GB)SA calculations.[43] SIE is a physics-based end-point scoring function that shares elements from the MM-PB(GB)SA methods.[43,44] The key distinction between both methods is that the former one was calibrated utilizing a diverse data set of 99 protein–ligand complexes in solution, while the MM-PB(GB)SA methods combine molecular mechanics calculations and continuum solvation models.[42] The free energy of binding between the ligand and protein is calculated by summation of Evdw (van der Waals interaction) and ECoul (Coulomb interaction) in the bound state, ΔGRF (reaction free energy between the bound and free states), and ΔSA (molecular surface area upon binding) as followswhere the coefficients ρ (AMBER van der Waals radii linear scaling coefficient), α (global proportionality coefficient relating to the loss of configurational entropy upon binding), γ (molecular surface area coefficient), Din (solute internal dielectric constant), and constant C were optimized with the set of 99 protein–ligand complexes. ρ = 1.6, α = 0.105, γ = 0.013 kcal mol–1 Å–2, Din = 2.25, and C = −2.89 kcal/mol are the optimized values.

Results and Discussion

System Stability of Simulation Models

The stability of inhibitor(s)/PLpro complexes was determined using root-mean-square displacement (RMSD, Figure ) plotted along a 100 ns simulation period. The RMSD outcome from the three independent simulations of each system showed consistent patterns. The RMSD values of the VIR250 system continuously increased in the first 20 ns and then fluctuated between ∼2.0 and 3.0 Å until the simulation ended. Whereas the RMSD values of other systems, VIR251, compound 3, and compound Y96, increased at 40 ns and fluctuated between ∼2.0 and 3.0 Å. In addition, the RMSD values of GRL-0617 increased during the first 40 ns and reached an equilibrium afterward. The distance between the center of mass (Cm) of inhibitor(s) and the Cm of active site residues of SARS-CoV-2 PLpro,[45] including the catalytic residues C111, H272, and D286 and the substrate-binding residues L162, M208, T301, P247, P248, Y264, Y268, and Q269, was further calculated to determine the protein–ligand stability. It was found that all studied inhibitors were relatively stable at the active site of SARS-CoV-2 PLpro from the beginning to the end of the simulation, supported by the MD snapshots of inhibitor(s)/PLpro complexes along the simulations time compared to the X-ray structures (Figure S2) and the movies generated from run1 (Videos S1–S5) showing that the binding mode of each type of ligand bound to SARS-CoV-2 PLpro was quite similar to the X-ray structures.

Figure 2

RMSD plot of protein backbone atoms (CA, C, O, and N), the distance between the Cm of inhibitor(s) and the Cm of active site residues, #H-bonds, and #atom contacts of inhibitor(s)/SARS-CoV-2 PLpro complexes. The minimized structure of each system was used as a reference structure for RMSD calculations.

Moreover, the number of hydrogen bonds (#H-bonds) and the number of atom contacts within 3.5 Å (#atom contacts) of the inhibitors within the SARS-CoV-2 PLpro binding site were also calculated. In the binding site, the values of #H-bonds and #atom contacts of all systems showed stability at the beginning until the end of the simulation; however, the last 20 ns MD trajectories of each system were used for further analysis. It was found that the peptidomimetic inhibitors (VIR series) gave the #H-bonds as well as #atom contacts more than naphthalene-based inhibitors (GRL-0617, compound 3, and compound Y96), suggesting that peptidomimetic inhibitors could interact with SARS-CoV-2 PLpro better than the other inhibitors (discussed in more detail in the next section).

Key Binding Residue Interactions

The key binding amino acids of inhibitor binding within SARS-CoV-2 PLpro were investigated. The ΔGbindresidue calculations based on the MM/GBSA method were performed using 100 snapshots derived from the last 20 ns of simulations. The energy contributions of <−1.0 kcal/mol of inhibitor orientation at the binding pocket of SARS-CoV-2 PLpro are displayed in Figure . The negative ΔGbindresidue value represents inhibitor stabilization, while positive values represent inhibitor destabilization. The average ΔGbindresidue values of each inhibitor binding to SARS-CoV-2 PLpro were obtained from the three independent simulations.

Figure 3

Per-residue decomposition free energy (ΔGbindresidue) of the inhibitor(s) within the SARS-CoV-2 PLpro active site domain. The data were derived from the last 20 ns of the average from three different simulations. The lowest and highest energies range from dark magenta to red, respectively. The orientations of the inhibitors were quite similar to the binding mode observed in the X-ray structures.

There are two major patterns of ΔGbindresidue of inhibitor(s) based on the types of inhibitors, peptidomimetics, and naphthalene-based inhibitors as follows: (i) the C-terminal of peptidomimetic inhibitors, VIR250 and VIR251, interacted with hydrophobic residues at the P2 site, L162 (green, −2.38 ± 0.02 and −2.61 ± 0.10 kcal/mol for VIR250 and VIR251, respectively) and G163 (dark blue, −3.37 ± 0.07 kcal/mol for VIR250 and 3.51 ± 0.16 kcal/mol for VIR251), and (ii) the naphthalene-based inhibitors (GRL-0617, compound 3, and compound Y96) showed low interaction bindings (>−1.0 kcal/mol) in these regions. Lengths of the peptidomimetic inhibitors are larger than those of naphthalene-based inhibitors, resulting in accommodating within the P1 site of SARS-CoV-2 PLpro. It was found that vinylmethyl ester (VME) at the C-terminal of both peptidomimetics located well at the catalytic site of the P1 site (Figure ).

In addition, both peptidomimetic inhibitors, VIR250 and VIR251, were stabilized by similar hydrophobic interactions with six hotspot residues including (i) P247 and P248 at the P4 site, (ii) Y264 and Y268 at the P3 site, and (iii) G271 and Y273 at P1 site. Interestingly, at the P4 site of SARS-CoV-2 PLpro, we found that the benzothiazole ring of VIR250 (yellow, −1.38 ± 0.04 kcal/mol) forms hydrophobic interactions with P248 slightly better than the phenol ring of VIR251 (yellow, −1.29 ± 0.02 kcal/mol) as well as interacts with Y264 (dark orange, −1.67 ± 0.02 kcal/mol for VIR250 and yellow, −1.10 ± 0.21 kcal/mol for VIR251) and BL2 loop Y268 residue (dark green, −2.04 ± 0.13 kcal/mol for VIR250 and orange, −1.31 ± 0.12 kcal/mol for VIR251) more than VIR251. In addition, vinylmethyl ester (VME) of VIR250 (yellow, −1.00 ± 0.05 kcal/mol) hydrophobically interacted with catalytic residue H272 better than VIR251 (−0.87 ± 0.05 kcal/mol).

The H272 residue is involved in SARS-CoV-2 PLpro activity to process the newly polypeptide chain of the virus. This residue paired with D286 functions as a general acid–base, which promotes deprotonation of C111 for catalytic mechanisms stabilizing with oxyanion hole Y106 residue.[19,26] Both peptidomimetic inhibitors interacted with the H272 catalytic residue. Therefore, they can inhibit the activity of SARS-CoV-2 PLpro following reduced propagation of the virus. However, from binding comparisons between VIR250 and VIR251, it could be noted that VIR250 stabilized within the pocket site more than VIR251. These findings agree well with a previous study reporting that the VIR251 (Kcat/Km = 1000 M–1 s–1) compound showed less selective inhibition against SARS-CoV-2 PLpro than the VIR250 compound (Kcat/Km = 600 M–1 s–1).[21]

For naphthalene-based inhibitors, compounds GRL-0617, 3, and Y96 formed hydrophobic interaction with four-spot residues (P247, P248, Y264, and Y268), which are the substrate-binding residues. Their naphthalene core is stabilized within the pocket with P247 and P248 at the P4 site, while their carbonyl group points to Y264, Y268, and Q269 at the P3 site. The flexible β-hairpin BL2 loop Y268 and Q269 residues are important for controlling viral protein substrate binding.[46] Three naphthalene-based inhibitors, GRL-0617 (deep pink, −3.98 ± 0.04 kcal/mol), compound 3 (deep pink, −4.20 ± 0.06 kcal/mol), and compound Y96 (deep pink, −4.06 ± 0.06 kcal/mol) formed strong hydrophobic interactions with Y268 to block substrate binding, which is consistent with the previous reports on GRL-0617. Other inhibitors/drugs that can inhibit PLpro activity by occupying this site include rac3j (N-benzyl-1-[1-(naphthalen-1-yl)ethyl]piperidine-4-carboxamide), amprenavir, pexidartinib, and indinavir.[47,48] Three naphthalene-based inhibitors (compounds GRL-0617, 3, and Y96) are small molecules resulting lack of P1 site binding (Figure ). These binding patterns obtained by MD simulations were similar to the X-ray structure. As a result, these inhibitors did not interact with catalytic residues C111 and H272. However, they can inhibit the SARS-CoV-2 PLpro activity, supported by previous reports on GRL-0617 and compound 3,[19,49] and interact well with the Y268 residue in the BL2 loop of SARS-CoV-2 PLpro better than peptidomimetic inhibitors, VIR250 and VIR251. Altogether, all inhibitors bind to the substrate peptide motif LXGG at P2–P3 sites.[14] The benzothiazole ring of VIR250, the phenol ring of VIR251, and the naphthalene core of GRL-0617, compound 3, and compound Y96 were stable at the P4 site (Figure ). Besides, the benzothiazole ring in run1 of VIR250 showed slight fluctuation at the P4 site. The positions of the benzothiazole ring of VIR250 and the phenolic ring of VIR251 obtained by MD simulations were closer to the pocket P4 site when compared to the X-ray crystal structure. It should be noted that the P4 site is preferable for hydrophobic interaction caused by the hydrophobic nature of the P4 pocket that is largely formed by residues P247 and P248.

The average energy contributions in terms of electrostatic (ΔEele + ΔGpolar) and vdW (ΔEvdW + ΔGnonpolar) interactions of the important residue for inhibitor(s) binding from three independent simulations are illustrated in Figure . The main contributor for all inhibitors’ stabilization was vdW interactions, which gave negative values (red lines) of up to ∼−4.50 kcal/mol. These vdW interactions showed values less than the electrostatic interactions (∼− 2.50 to ∼1.60 kcal/mol represented by black lines). Our finding was consistent well with other inhibitors of SARS-CoV-2 PLpro, for example, phthalazinone derivative (ADM13083841),[50] Ebselen, and 4′-O-methylbavachalcone.[51,52] Among the residues of peptidomimetic inhibitors/SARS-CoV-2 PLpro complexes, the substrate residue binding Q269 in the BL2 loop of both systems showed the lowest of vdW interactions. For naphthalene-based inhibitors, Y268 of SARS-CoV-2 PLpro is the key important residue via vdW interactions (Figure ).[15] GRL-0617 could specifically bind to SAR-CoV-2 PLpro (IC50 of 2.3 μM) but not to MERS-PLpro. Mutation of SARS-CoV-2 PLpro Y268 to T268 and G268 showed that GRL-0617 was unable to inhibit PLpro. Therefore, it was concluded that in MERS-PLpro, GRL-0617 could not inhibit MERS-PLpro, possibly because the amino acid at position 268 was a threonine residue instead of tyrosine.[49,53]

Figure 5

Electrostatic and van der Waals (vdW) energy contributions of inhibitor(s) binding within SARS-CoV-2 PLpro. The data was derived from the last 20 ns of the average three independent simulations.

Although the vdW interactions were the main contributors for all inhibitors, some residues were also stabilized by electrostatic interactions. For instance, the backbone amide group of the G163 residue could form strong H-bonds (>90%) with backbone amide group of both VIR250 and VIR251 (Figure ), resulting in stronger electrostatic interaction (∼− 2.50 kcal/mol represented by black lines) than vdW interaction (∼−1.00 kcal/mol represented by red lines), as shown in Figure . In addition, the electrostatic interaction of Q269 against GRL-0617 (−1.13 kcal/mol) was stronger than that against compound 3 (−0.23 kcal/mol) and compound Y96 (−0.14 kcal/mol) since hydrogen bond formation of Q269 with GRL-0617 was stronger than those with the two naphthalene-based inhibitors (Figure ). Notably, among the studied naphthalene-based inhibitors, GRL-0617 showed the highest binding ability with SARS-CoV-2 PLpro (discussed in more detail later).

Figure 6

Hydrogen bond occupation (%) of donor···acceptor atoms with the important residues of inhibitor(s) within the substrate-binding region of SARS-CoV-2 PLpro. The data were derived from the last 20 ns of three different simulations. Data are shown as means ± the standard deviation (SD).

Hydrogen Bonding of SARS-CoV-2 PLpro and Inhibitor(s)

Hydrogen bonds (H-bonds) are crucial in biological systems for a wide range of factors, such as helping for stabilizing the ligand within proteins and protein folding. The average intermolecular H-bonds formed between the inhibitors and residues of SARS-CoV-2 PLpro from the last 20 ns of simulations are plotted in Figure , and only the strong H-bonds (>90%) of five inhibitors within SARS-CoV-2 PLpro are displayed. For peptidomimetic inhibitors, the backbone amide group of both VIR250 and VIR251 forms H-bonds with substrate-binding residues, including N4H···G163@O at the P2 site (99.65 ± 0.23% for VIR250 and 99.72 ± 0.18% for VIR251) and N5-H···G271@O at the P1 site (99.25 ± 0.64% and 95.03 ± 3.02% for VIR250 and VIR251, respectively). These MD results are similar to previously reported SARS-CoV-2 PLpro in complex with amprenavir, an HIV-1 protease inhibitor, and semagacestat, a drug in phase III clinical trials for Alzheimer’s disease from the glide docked structure, and both VIR250 and VIR251 from the X-ray crystal structure.[21,47] Moreover, the backbone oxygen atom (O4) in the carbonyl group of VIR250 could form H-bond with G163@N–H (95.65 ± 1.59%), which is the substrate-binding site located near the P1 catalytic site.

The decomposition free energy (Figure ) and energy contribution (Figure ) of both peptidomimetics showed high stabilization within the pocket with the G163 residue at the P2 site in a similar range. However, VIR250 could form strong H-bonds with three atoms including G271@O, G163@O, and G163@N–H, which was slightly more than the H-bond formation of VIR251 with two atoms, G271@O and G163@O. Therefore, VIR250 is suitable for SARS-CoV-2 PLpro inhibitor than VIR251.

The naphthalene-based inhibitors (99.70 ± 0.26% for GRL-0617, 99.38 ± 0.13 for compound 3, and 99.82 ± 0.03% for compound Y96) were found to form strong H-bonds with the substrate-binding residue of Q269@N–H···O1 at the BL2 loop in the P3 site. This Q296 residue also formed H-bonds with gallocatechin gallate[54] and disulfiram,[52,55] which are inhibitors of SARS-CoV-2 PLpro. The P3 site is a wide pocket adjacent to the oxygen atom of the carbonyl groups of Y268 and Q269 residues, where hydrogen bonds can be formed at this site, as well as hydrophobic interactions with the side chain of the Y268 residue.[21] The peptidomimetic inhibitors (VIR250 and VIR251) showed formation of a higher number of hydrogen bonds with SARS-CoV-2 PLpro residues than naphthalene-based inhibitors (GRL-0617, compound 3, and compound Y96), in good agreement with #H-bonds along with the simulations (Figure ). However, GRL-0617, compound 3, and compound Y96 could stably accommodate within P4 and P3 pocket sites (Figures and 4). Altogether, among all studied inhibitors, VIR250 is the best one for SARS-CoV-2 PLpro inhibition.

Binding Affinity Prediction

The binding affinities of five inhibitors against SARS-CoV-2 PLpro were predicted by MM-PB(GB)SA and the solvated interaction energy (SIE) method derived from 100 snapshots that were extracted from the last 20 ns of simulations (Table ). It was found that the electrostatic interactions play a major role in complexation between peptidomimetic inhibitors (ΔEele’s of −122.91 ± 1.38 for VIR250 and −118.15 ± 1.38 for VIR251) and SARS-CoV-2 PLpro rather than vdW interactions (ΔEvdW’s of −56.84 ± 0.39 and −54.27 ± 0.36 for VIR250 and VIR251, respectively). This is because peptidomimetic inhibitors, VIR250 and VIR251, have an amino group side chain that could be protonated at physiological conditions (pH 7.4), resulting in providing the predominant electrostatic interactions. On the other hand, vdW interactions were the main force for naphthalene-based inhibitors’ binding to SARS-CoV-2 PLpro (ΔEvdW’s of −37.70 ± 0.26 for GRL-0617, −39.28 ± 0.26 for compound 3, and −38.10 ± 0.27 for compound Y96).

Table 1

Energy Components (ΔGbind, kcal/mol) of SARS-CoV-2 PLpro/inhibitor(s) Complexes in Terms of MM-PB(GB)SA and SIE Methods Compared to the Experimental Result Previously Reporteda

energy component (kcal/mol)	VIR250	VIR251	GRL-0617	compound 3	compound Y96
Gas Term
ΔEele	–122.91 ± 1.38	–118.15 ± 1.38	–19.59 ± 0.51	–21.57 ± 0.89	–16.90 ± 0.48
ΔEvdW	–56.84 ± 0.39	–54.27 ± 0.36	–37.70 ± 0.26	–39.28 ± 0.26	–38.10 ± 0.27
ΔEMM	–179.75 ± 1.46	–172.42 ± 1.51	–57.29 ± 0.55	–60.85 ± 0.91	–54.99 ± 0.53
TΔS	–32.53 ± 0.89	–32.17 ± 0.97	–20.72 ± 1.00	–22.01 ± 1.06	–19.47 ± 0.74
Solvation Term
PBSA
ΔGsol(PBSA)nonpolar	–7.62 ± 0.02	–7.29 ± 0.02	–5.32 ± 0.01	–5.68 ± 0.02	–5.44 ± 0.02
ΔGsol(PBSA)ele	141.18 ± 1.16	137.63 ± 1.30	32.87 ± 0.44	38.35 ± 0.90	32.70 ± 0.55
ΔGsol(PBSA)	133.55 ± 1.15	130.33 ± 1.29	27.55 ± 0.44	32.67 ± 0.89	27.26 ± 0.54
GBSA
ΔGsol(GBSA)nonpolar	–7.16 ± 0.02	–7.11 ± 0.02	–4.38 ± 0.01	–4.63 ± 0.02	–4.49 ± 0.02
ΔGsol(GBSA)ele	139.36 ± 1.15	135.51 ± 1.27	29.77 ± 0.40	34.62 ± 0.83	28.24 ± 0.44
ΔGsol(GBSA)	132.20 ± 1.15	128.40 ± 1.26	25.39 ± 0.39	29.99 ± 0.82	23.75 ± 0.44
Gas Term + Solvation Term
ΔEvdW + ΔGsol(PBSA)nonpolar	–64.46 ± 0.39	–61.56 ± 0.36	–43.02 ± 0.26	–44.96 ± 0.26	–43.54 ± 0.27
ΔEvdW + ΔGsol(GBSA)nonpolar	–64.00 ± 0.39	–61.38 ± 0.36	–42.08 ± 0.26	–43.91 ± 0.26	–42.59 ± 0.27
ΔEele + ΔGsol(PBSA)ele	18.27 ± 1.80	19.48 ± 1.90	13.28 ± 0.67	16.78 ± 1.27	15.80 ± 0.73
ΔEele + ΔGsol(GBSA)ele	16.45 ± 1.80	17.36 ± 1.88	10.18 ± 0.65	13.05 ± 1.22	11.34 ± 0.65
Binding Free Energy
ΔGtotal(PBSA)	–46.21 ± 0.54	–42.10 ± 0.57	–29.75 ± 0.30	–28.19 ± 0.32	–27.73 ± 0.31
ΔGtotal(GBSA)	–47.56 ± 0.50	–44.02 ± 0.5	–31.90 ± 0.27	–30.87 ± 0.24	–31.24 ± 0.25
ΔGbind(MM/PBSA)	–13.68 ± 0.67	–9.92 ± 0.71	–9.02 ± 0.53	–6.18 ± 0.61	–8.26 ± 0.44
ΔGbind(MM/GBSA)	–15.02 ± 0.64	–11.85 ± 0.45	–11.18 ± 0.99	–8.86 ± 0.11	–11.77 ± 0.42
SIE	–10.04 ± 0.05	–9.55 ± 0.09	–7.08 ± 0.03	–7.06 ± 0.03	–7.00 ± 0.03
IC50 (μM)	1.6021	1.6521	2.349	6.449
ΔGexpa	–8.22	–8.20	–8.00	–7.37
Kcat/Km (M–1 s–1)	60021	100021

Binding free energies from the experiment (ΔGExp) were converted by the Cheng–Prusoff equation of ΔG =RT ln(IC50).[58] Data from triplicate simulations are shown as mean ± standard error of mean (SEM).

The summation of the gas term and solvation free energy (ΔEvdW + ΔGsolnonpolar and ΔEele + ΔGsolele) of peptidomimetics VIR250 and VIR251 from PB and GB models was calculated. It is worth noting that the vdW interaction was favorable to the overall binding free energies of two peptidomimetic inhibitor PLpro complexes. This evidence is in agreement with other protease inhibitors such as faldaprevir/HCV NS3/4A protease, which is a drug for hepatitis C virus (HCV),[56] and the fonsecin/SARS-CoV-2 PLpro complex, which is a fungal secondary metabolite.[57]

Among five compounds, we found that the MM-PB(GB)SA of peptidomimetic inhibitor VIR250 shows the highest binding affinity (ΔGbind’s of −13.68 ± 0.67 and −15.02 ± 0.64 kcal/mol for PBSA and GBSA, respectively) followed by VIR251 affinity (ΔGbind’s of −9.92 ± 0.71 and −11.85 ± 0.45 kcal/mol for PBSA and GBSA, respectively), which is rather than two naphthalene-based inhibitors, GRL-0617 (ΔGbind’s of −9.02 ± 0.53 kcal/mol for PBSA and −11.18 ± 0.99 kcal/mol for GBSA) and compound Y96 (ΔGbind’s of −8.26 ± 0.44 kcal/mol for PBSA and −11.77 ± 0.42 kcal/mol for GBSA). In contrast, the remaining naphthalene-based inhibitor of compound 3 showed the highest binding affinity within the SARS-CoV-2 PLpro binding pocket (ΔGbind’s of −6.18 ± 0.61 and −8.86 ± 0.11 kcal/mol for PBSA and GBSA, respectively). In addition, SIE-based calculations of all inhibitors gave a somewhat similar trend of MM-PB(GB)SA binding affinity calculations. These predicted ΔGbind values are somewhat consistent with the reported experimental data ranked in order VIR250 > VIR251 > GRL-0617 > compound Y96 > compound 3[21,49] (Table ).

Solvent Accessibility of the SARS-CoV-2 PLpro/inhibitor(s) Complex

Solvent accessibility within 4 Å in the binding pocket site of SARS-CoV-2 PLpro without the inhibitor (apo form) and SARS-CoV-2 PLpro/inhibitor(s) complexes (holo form) along 100 ns of the simulation was calculated by the solvent-accessible surface area (SASA). The SASA result is illustrated in Figure , and the average SASAs from the last 20 ns of MD simulations are listed in the purple text. The holo forms consisting of VIR250, VIR251, GRL-0617, compound 3, and compound Y96 (877–1019Å2) showed lower SASA values than the apo form (1212 Å2) since the holo form of the ligands could well occupy within the binding pocket. This phenomenon is corresponding to the previous reports on phthalazinone derivatives[50] and gallocatechin gallate[54] bound to SARS-CoV-2 PLpro, which show a decrease in water accessibility when inhibitors bind to the enzyme active site. The water accessibility of all inhibitors bound SARS-CoV-2 PLpro ranged in the order VIR251 (877.84 ± 42.13 Å2) < VIR250 (886.78 ± 44.84 Å2) < GRL-0617 (947.61 ± 18.93 Å2) < compound Y96 (967.30 ± 21.19 Å2) < compound 3 (1019.22 ± 48.18 Å2). This implies that the peptidomimetic inhibitors could bind to the binding pocket of PLpro better than naphthalene-based inhibitors since naphthalene-based inhibitors lack interaction with the P1 site, while both VIRs250 and 251 peptidomimetic inhibitors could interact with all four sites (P1–P4) of SARS-CoV-2 PLpro. Therefore, the peptidomimetic inhibitors were selected to modify for enhancing the binding efficiency with SARS-CoV-2 PLpro.

Figure 7

SASAs of the residues surrounding atoms within the 4 Å sphere of inhibitor(s). The purple text represented the average values of SASAs derived from the last 20 ns of three independent simulations.

Rationale for SARS-CoV-2 PLpro Inhibitor Design

The peptidomimetic VIR250 inhibitor was selected as a template for developing a novel SARS-CoV-2 PLpro inhibitor since it exhibits the highest vdW contributions with residues Y248, Y264, and 268 (Figure ). In addition, VIR250 shows the strongest hydrogen bond formation with three atoms of G163 (P2 site) and G271 (P3 site) residues (Figure ) and gave the highest binding efficiency to SARS-CoV-2 PLpro (Table ). To improve the efficiency of binding of VIR250, we suggested that the functional group in the side chains should be modified as follows (Figure A): (i) changing the benzothiazole ring in the P4 site to increase the hydrophobic interaction with residues P247 and P248; in addition, heteroatoms that could form hydrogen bonds with the P2 site (R166 and E167 residues) should be fused in this ring (e.g., flavone, azaindole, quinoxaline, phthalazine, pyrrolopyridine, and thiazolopyridine); and (ii) enhancing the nonpolar moieties at the N-terminal of VIR250 in the P3 site, which interacted with the BL2 loop, helix 5 (Y264 and Y268),[46] and proximately hydrophobic residues at P1 (Y273 residue) and P2 (L162 residue) sites. In addition, small functional groups that could form hydrogen bonds with the Y268 and Q269 should be added (e.g., methyl, ethyl, and ethene). However, the NH2 side chain located at the solvent-exposed regionshould be conserved. Furthermore, vinylmethyl ester (VME) at the C-terminal at the P1 site should be preserved as hydrophobic interaction with catalytic H272 residue (Figure ).

Figure 8

Rational design of the SARS-CoV-2 PLpro inhibitors based on the peptidomimetic VIR250 inhibitor. (A) 3D and 2D structures of VIR250 with ligand modifications, (B) 2D structure of modified VIR250 and their binding free energy prediction comparison with VIR250 against SARS-CoV-2 PLpro derived from MM-PB(GB)SA calculations, (C) hydrogen bond occupation, and (D) per-residue decomposition free energy (ΔGbindresidue) of modified VIR250 and the A2/SARS-CoV-2 PLpro complex. Calculations are obtained from one snapshot of the complex after system minimization and solvation in the TIP3P model.

Three novel VIR250 derivatives were proposed by modification of the VIR250 crystal structure (Figure B). The VIR250 derivative structures were optimized at the HF/6-31g(d) level. Subsequently, these complexes were solvated in the TIP3P model water and then minimized the system to stabilize the complex structure. The MM-PB(GB)SA method was used to evaluate their binding affinity for PLpro from one snapshot comparison with VIR250 from structure minimization. It was found that the binding affinities from MM-PBSA calculations of three compounds including A1 (−19.31 kcal/mol), A2 (−18.91 kcal/mol), and A3 (−18.01 kcal/mol) showed stronger binding within PLpro than VIR250 (−12.65 kcal/mol). In addition, the binding from MM-GBSA, all new derivatives, A1 (−22.99 kcal/mol), A2 (−27.16 kcal/mol) and A3 (−24.86 kcal/mol) also gave stronger binding than VIR250 (−13.14 kcal/mol). Among novel VIR250 derivatives, compound A2 was selected to study the binding behavior since it showed the lowest binding free energy from the MM-GBSA method. Interestingly, hydrogen bond formation between the thiazolopyridine ring (N atom) and R166@N-HE at the P2 site was observed (Figure C). In addition, the hydrophobic interactions based on the MM/GBSA method of compound A2 with residues Y264 (ΔGbindresidue of −3.72 kcal/mol, deep pink) and Y268 (ΔGbindresidue of −5.33 kcal/mol, violet red) at the P3 site (Figure D) were greatly increased as suggested.

Conclusions

In this study, all-atom MD simulations were used to explain the main binding of the two peptidomimetic inhibitors (VIR250 and VIR251) and three naphthalene-based inhibitors (GRL-0617, compound 3, and compound Y96) against SARS-CoV-2 PLpro. Peptidomimetic inhibitors, VIR250 and VIR251, form hydrophobic interaction with four important residues sites including (i) P1 site, G271, and Y273; (ii) P2 site, L162, and G163; (iii) P3 site, Y264, and Y268; and (iv) P4 site, P247, and P248. In addition, VIR250 hydrophobically interacted with the H272 catalytic residue more than VIR251, while the three naphthalene-based inhibitors showed low interactions with the P1 site. However, three inhibitors form strong hydrophobic interaction with Y268 at the BL2 loop and the naphthalene core stabilizes with P247 and P248. Both peptidomimetics stabilize the binding by forming hydrogen bonds with three sites (P1–P3), whereas strong hydrogen bonds at the P2 site of naphthalene-based inhibitors were found. Besides, the vdW interactions are the main drivers for five inhibitors interacting within PLpro. The binding affinities for PLpro of all inhibitors agreed with the experimental results in previous reports, which were ranked in order VIR250 > VIR251 > GRL-0617 > compound Y96 > compound 3. Thiazolopyridine enhanced binding with R166, P247, and P248 residues of the P2 and P4 sites for the rational drug design of VIR250. In addition, the hydrophobic interaction at the P3 site increased by interacting with the BL2 loop and nearby residues at P1 (Y273) and P2 (L162) sites. The obtained information from this work is useful for understanding the inhibition of SARS-CoV-2 PLpro at the atomic level, which could be the basis for further antiviral drug development.