The present study introduced a novel concept titled ‘peptide mold’ for hit discovery for a potential anti-SARS-CoV-2 target, PLpro enzyme. The study used biomimicry strategy based on tetrapeptides. The protease active site is tuned to identify specific amino acid sequences of the substrate protein for initial, highly-specific binding, followed by the precise catalytic events leading to selective peptide bond cleavage. We started off with a thought process that if we dock smaller peptides, e.g., tetrapeptides, in the PLpro active site, we would be able to identify potential binders based on their shape- and electrostatic complementarity with the active site residues. The tetrapeptides would contain the side chains of varied nature (polar/nonpolar, acidic/neutral/basic, hydrophilic/hydrophobic, aliphatic/aromatic, H-bond Don/H-bond Acc) which could also take up varied shapes due to the conformational freedom (dihedral angles θ and ψ). The thermodynamically-allowed conformational freedom and the nature of side chains of the constituent amino acids would make them ideal candidates for exploring the complimentary enthalpic and entropic contributions to the overall binding energy via ligand-receptor interaction forces, such as electrostatic, H-bonding, van der Waals, hydrophobic, cation-π, π-π stacking and few others. We attempted to mimic nature’s minimalistic design approach, i.e., we tried to imbibe biomimicry in our molecular design strategy.
Initial thought process was based on three options for smaller peptides, namely, tripeptides, tetrapeptides and pentapeptides, keeping in mind the overall size (molecular weight) of the resulting peptides. The molecular weights of the smallest tripeptide (Gly-Gly-Gly), tetrapeptide (Gly-Gly-Gly-Gly) and the pentapeptide (Gly-Gly-Gly-Gly-Gly) are 189.17, 246.22 and 303.27 Da, respectively. The corresponding values for the largest tripeptide (Trp-Trp-Trp), tetrapeptide (Trp-Trp-Trp-Trp) and pentapeptide (Trp-Trp-Trp-Trp-Trp) are 576.65, 762.87 and 949.08 Da, respectively. For a drug-like chemical space, the tetrapeptides were obviously the best option. Another reason to choose a tetrapeptide library was due to the preference of PLpro for a specific tetrapeptide sequence, Ac-LRGG-ACC, for recognition and cleavage; the same sequence is also recognized by other deubiquitinating enzymes. We, therefore, initiated the generation of tetrapeptide library consisting of all possible combinations of naturally-occurring L-amino acids (absolute stereochemistry S) (Chart 1).
The processed tetrapeptide library of 32 million conformations derived from 160,000 tetrapeptides was docked into the active site of PLpro. The native ligand was covalently bound to the active site nucleophilic Cys111 (cleavage site). For the docking tool to identify Het groups, thereby identifying the site for docking, the covalent bond with Cys111 was broken and the corresponding Hs added. The generated receptor unit was used for docking of the multiconformer tetrapeptide library. The task was broken into smaller subjobs, consisting of 8000 tetrapeptides (160,000 conformations); each docking run generated top-500 hits (including multiple conformations of same molecule). These 10,000 (500X20) docked poses were further docked using same workflow to select the top-500 poses. Table 1 lists the top-10 tetrapeptides with their corresponding FRED Chemguass4 scores.
As seen from Table 1, the top-10 tetrapeptides exhibited unique structural features such as – 1. Presence of Gly as the C-terminus residue, 2. Presence of Gly as the second residue in 9/10 cases, 3. presence of Ala/Ser as the N-terminus residue in 8/10 cases, and 4. presence a hydrophobic residue such as Phe (4/10), Leu (2/10), Thr (2/10) and Tyr (2/10 cases). Interestingly, the recognition sequence specific to PLpro cleavage site features were retained by the top-scoring tetrapeptides. The average relative size (mol wt) of the top-scoring tetrapeptides was 364.18, although bit skewed by peptide No. 9 (Table 1). The peptide mold workflow, thus, could successfully identify a unique set of tetrapeptides in the top-scoring poses.
The docked pose of top-scoring tetrapeptide 1 (Table 1) in the active site of PLpro is shown for reference (Fig. 2). The active site shape was adapted by 1, which interacted with PLpro via four H-bonds with Gly163 (two H-bonds with backbone –NH and –C = O groups), Gly271 and Tyr268 (one H-bond each with backbone –NH and –C = O groups, respectively). Additionally, the N-terminus –NH3 of tetrapeptide formed a strong salt-bridge interaction with anionic side chain of Asp164. The 3rd residue of 1 was exposed to solvent. The crystal structure ligand and 1 overlapped in the PLpro active site in similar fashion (pose not shown). The C-terminal Gly of 1 was located in proximity to the catalytic Cys111, confirming the appropriate occupancy of the catalytic site of PLpro.
The top-500 docked tetrapeptides were further docked in the PLpro active site using Glide Extra Precision (XP) mode. The top-scoring pose of each of the tetrapeptides was sorted as per Glide XP score. The top-100 docked poses were subjected to diversity analysis based on MACCS fingerprints. The top-10 most diverse tetrapeptides were then taken up for pharmacophore elucidation using Phase v3.5. A five-point pharmacophore comprising of AAADN features, was generated using default settings. Figure 3 depicts these pharmacophoric features and their geometric arrangement, including directions of the A and D features. The top-five hits from Phase pharmacophore screening were mapped onto the five-point pharmacophore and are shown in Figs. 4a-e. Other related pharmacophore hypotheses were not discussed here due to closeness (in terms of number of features and their nature) with the top-scoring hypothesis. The selected pharmacophore hypothesis (AAADN) was used for virtual screening of the multiconformer phytochemicals library using Phase v3.5 (default settings). In summary, the active site of PLpro was successfully mapped in the form of 5-point pharmacophore, ready to be deployed for the hit identification in the ensuing virtual screening exercise.
The molecular structures of top-five hits, 11–15 from Phase virtual screening of multiconformer phytochemicals library are shown in Table 2, along with representative calculated/predicted molecular, physicochemical and pharmacokinetic properties. These hits were docked into the catalytic site of PLpro using MOE 2022.02. A total of five poses of each of the hits were generated and further examined for their interactions with residues lining the PLpro catalytic site. Representative poses (hits 11 and 12, Table 2), along with the 2D interaction diagrams are shown in Figs. 5–6 (PDB ID 6WUU) and 7–8 (PDB ID 6WX4). The native ligands were docked in their respective crystal structures for optimizing the docking protocol settings. The docked poses and relative arrangement the 2D interaction diagrams of both the crystal structure ligands can be found in the Supporting Information section (Figs. 1S and 2S). The corresponding docking scores are shown in Table 3.
As seen from Fig. 5, the two top-scoring poses of 11 in the active site of PLpro (PDB ID 6WUU) were flipped; in pose 1, the –CH2OH donor of 11 formed a H-bond with Asp164 and secondary 4’–OH of the sugar moiety interacted with Tyr268 via H-bonding, while in pose 2, the 2’-OH of 11 formed a H-bond with Asp164. Given the polar nature of the proximal and distal portions of 11, the original and flipped poses did well in compensating the binding energy by forming multiple interactions with the active site residues. The salt-bridge interaction of the –COOH of 11 with nearby Arg166 was absent in both the top-scoring poses. This could be attributed to the fairly rigid framework of 11. The situation was bit different for the other crystal structure of PLpro (PDB ID 6WX4). In Fig. 7, both the top-scoring poses of 11 were found to be similar with subtle difference in the orientation of the terminal groups. Interestingly, the –COOH group of 11 formed a salt-bridge interaction with Lys157 (both poses). The polyol –OH groups interacted with PLpro active site residues with fewer H-bond interactions with Asp164, Thr301, and Tyr268. Overall, these multiple interactions led to slightly better docking score in the second crystal structure (PDB ID 6WX4). The docking scores of 11 were less than the corresponding scores of the native ligands of both the crystal structures (Table 3).
On similar lines as 11, both the top-ranking poses of 12 showed flipped orientation in the active site of PLpro (PDB ID 6WUU) (Fig. 6). The sugar moiety was observed to form couple of H-bonds with Met208 and Asp302 (pose 1), while the pose 2 formed a lone H-bond with Arg166. The propionyl moiety of pose 1 (guache conformation around C2-C3 bond) formed a salt-bridge interaction with Arg166 (Fig. 6b); a similar interaction was absent in pose 2 (anti conformation around C2-C3 bond) (Fig. 6c) due to flipping. In both poses, the propionyl moiety was exposed to solvent. Additionally, the benzofuran ring formed arene-H interaction with Asp164 only in pose 1 (Fig. 6b). The docked poses of 12 in the PLpro structure (PDB ID 6WX4) (Fig. 8) exhibited similar trend with respect to the flipped poses. While the propionyl side chain formed H-bonds with key residue, Cys111, the sugar moiety interacted with Thr301 via H-bond (Fig. 8b). In pose 2, the sugar moiety formed extensive H-bond network with Cys111, Cys270, Gly 271 and Tyr268 and the propionyl group was involved in arene-H interaction with Tyr264 (Fig. 8c). Overall, both poses 1 and 2 of 12 were distinct due to the moderate flexibility to the molecular structure. The docked poses of hits 13 to 15 in both crystal structures are given in Supporting Information section (Figs. 3S to 8S).
Table 2 lists all the five hits, their structures and relevant information (CAS number, molecular formula and weight) along with the calculated/predicted molecular, physicochemical, pharmacokinetic and toxicity properties. The hits were all derived from originally phytochemical library. Due to the presence of significant proportion of heteroatoms in the molecular composition, the hits were appreciably polar, violated the Lipinski’s Rule of 5 filter and showed poor ADME predictions. This was no surprise given the nature of the substrate-binding site of the PLpro enzyme, evolved to identify peptides. The idea is to further transform the hits into more drug-like leads and optimize them along the way. The interesting part is that the hits were derived from a novel strategy, i.e., peptide mold, which can potentially be applied to any protein target, leave aside the proteases such as PLpro. The first hit, 11 (Geniposidic acid, Table 2) belongs to iridoid glycoside class, found in significant number of plants. It possesses multiple pharmacological activities, with main activity being anticancer.61 A recent article reports its antiviral activity against white spot syndrome virus via reducing STAT (signal transducer and activator of transcription) gene expression, thereby reducing viral replication.62 Compound 11 has been studied as tool molecule for white spot syndrome virus replication inhibition. The second hit, 12, is a ZINC database member (ZINC39010597)63a, which has not been tested in any assay. The activity predictions using SwissTargetPrediction64, indicated that it was likely to be a transporter inhibitor, without any peculiar details on the transporter.
Compound 13 (Table 2), a dihydroxybenzamide derivative of Ala-Ser dipeptide, is present in ZINC database (ZINC2356677764)63b and has not been evaluated in any in vitro biological assay. The target prediction for 13 using SwissTargetPrediction64 yielded caspase-1 as the top-raking target. Further literature searches in this direction led to an interesting fact that few viruses are dependent on the host caspases for proteolytic processing of the viral components.65 A caspase inhibitor is likely to demonstrate antiviral effects. On the other hand, caspases are Zn2+-dependent proteolytic enzymes.66 The N-terminal modified Ala-Ser, i.e., 13, is likely to bind in the active site of caspases. Inhibition of host caspases by a potential ligand may disrupt virus particle processing. In short, 13, a virtual screening hit, can be a potential PLpro inhibitor by virtue of its pharmacophoric features, and it may also hit other proteolytic enzymes, such as caspases, which in turn, may indirectly alter viral spread and infection. The cyclic pentapeptide, 14 (WS-7338A), is unique due to its structural features. It is reported to exhibit potent endothelin receptor antagonist activity.67 A recent report discussed the potential utility of endothelin receptor antagonist as antiviral entities.68 The hit 15 has been reported as dipeptidyl aminopeptidase IV (DPP IV, EC 3.4.14.5) inhibitor long ago.69 Recently, it has been argued whether DPP IV inhibitors could be used for the development of novel viral entry inhibitors for emerging viral infections.70 To summarize, the hits from virtual screening of phytochemicals library have a potential to exhibit antiviral action, most likely, via PLpro inhibition. The possibility of these agents to hit other proteolytic enzymes is likely to exist due to their peptidic (13–15) or non-peptidic (11 and 12) nature. All the five high-quality hits have originally been screened during pharmacophore search, followed by molecular docking.
Table 3
Docking score of the five hits from pharmacophore screening along with the crystal structure ligands
Pharmacophore Screening Hit/Crystal structure liganda
|
Docking score (PDB ID 6WUU)
|
Docking score (PDB ID 6WX4)
|
6WUU ligand
|
-10.64
|
ND
|
6WX4 ligand
|
ND
|
-9.78
|
11
|
-7.15
|
-7.16
|
12
|
-7.56
|
-7.47
|
13
|
-6.82
|
-7.51
|
14
|
-7.53
|
-8.14
|
15
|
-7.62
|
-7.43
|
aDocking score of the top-ranking pose |
To obtain better insights into the interactions between 13 and 15 with SARS-CoV-2 PLpro, we performed all-atom molecular dynamics simulations in an explicit water box. Furthermore, the stability of the protein ligand was estimated using an end-state free-energy method like MM-GB/SA on the MD sampled conformations. We used the backbone RMSD to understand the progression of the MD simulations and backbone RMSF for the overall dynamics of the protein-ligand complexes. The SARS-CoV-2 PLpro protein was found to be very flexible as shown by the backbone RMSD (Fig. 9A) and the backbone RMSF (Fig. 9C). Upon closer inspection, we observed that large fluctuations were observed in the C-terminal region from amino acid numbers 170–320. During the simulation, the protein significantly deviates from the starting conformation until 100 ns. The normalized frequency distribution shows that a large population of the protein conformations fluctuates around RMSD between 2–3 Å which indicates significant deviation from the starting geometry. We hypothesize that it could be attributed to the ligand binding wherein the protein samples conformational states are distinct from the starting geometry. To understand the strength of the binding of the ligands, we used the MMGBSA method to compute the binding energies, due to computational cost we excluded the entropy calculations, and therefore, the binding energies we obtained shall be treated as effective binding energies (ΔGeff). Our MD simulations suggest that 15 binds with far superior binding affinity than 13 (Fig. 9B). Furthermore, we observed several unbinding events for 13, as a result, the mean ΔGeff for 13 is -22.61 kcal/mol less than 15. Compound 15 is well suited to occupy the SARS-CoV-2 PLpro binding site as the proline at the centre of the 15 creates a bend in the backbone that enhances the binding stability. This structural feature is absent in 13 which may have caused poor affinity towards SARS-CoV-2 PLpro.
Figure 10 shows the docked pose of 15 following MD simulations along with the 2D interaction diagram and the interaction surface for clear visibility of the substrate-binding site of PLpro. As seen in Fig. 5S (Supporting Information section) and Fig. 10, the interactions of 15 with binding site residues are conserved. The amino terminal group formed electrostatic interactions with Asp164 and Asp302, along with H-bonding interactions with Thr301. The other electrostatic interaction of anionic carboxylate of 15 with Arg166 (Fig. 5S) was absent in the post-MD simulation pose (Fig. 10). A new H-bond was formed between C-terminal –C = O of 15 and Tyr273 (Fig. 10) which was not observed in the pre-MD pose (Fig. 5S). Additionally, one of the methylene groups of pyrrolidine ring of 15 interacted with Tyr264 via arene-H interaction (Fig. 10). The superposed active site residues of the substrate-binding site residues along with the pre- and post-MD simulations are shown in Fig. 11 for comparison. As seen from Fig. 11, the MD simulations led to moderate changes in the backbone and side-chain atoms of the substrate-binding site residues, thereby optimizing the interactions with 15, which were stable over the simulation intervals, giving rise of gain in binding energy.
The direct docking of the phytochemicals library in the substrate-binding site SARS-CoV-2 PLpro (PDB ID 6WX4) and the subsequent binding free energy estimations were performed to compare the utility of the peptide mold strategy over the conventional structure-based methodology. The corresponding results for the top-five hits are shown in Table 4. Hits 16 and 20 were amphoteric, hits 17 and 19 were anionic and hit 18 was neutral. Their molecular weights were > 500 Da, except 16. Since the hits were arranged in the ascending order of binding free energy values, it was no surprise that the bigger molecules yielded higher binding free energy estimates. These hits were docked in similar manner as hits from pharmacophore screening in the substrate binding site of PLpro (PDB IDs 6WUU and 6WX4) and compared. The results are listed in Table 5. The top-ranking docked poses 1 and 2 along with their 2D interaction diagrams of hit 18 in the substrate-binding site of PLpro (PDB ID 6WX4) are shown for comparison in Fig. 10. Table 5 also lists the H-bonds formed between the top-ranking pose of the hits and PLpro substrate-binding site residue. The corresponding docked poses and the 2D interaction diagrams of hits 16, 17, 19 and 20 in PLpro site (PDB ID 6WX4) are given in Supporting Information section (Figs. 9S to 12S).
Hit 18 formed three H-bonds (pose 1) with Gly163, Gly271 and Cys270 (Fig. 10), while in pose 2, additional H-bond with Tyr268 was observed. Since 18 lacked anionic group (unlike other hits), no salt bridge interaction with PLpro cationic residues was seen. The tetrahydrofuro[2,3-d][1,3]dioxol-6-ol moiety was exposed to solvent at the proximal end of the substrate-binding site of the enzyme. Both the poses exhibited similar binding modes with subtle differences in the orientation of the benzyl group. No information of 18 could be retrieved from literature except the CAS No.: 1191113-39-9. It has not been previously used for any purpose at all. Hit 16 (CAS No. 455325-25-4) was reported in one research paper two decades ago.71 Hit 17 and 20 were not reported in the chemical literature at all. Though the CAS No. was available for hit 19 (1241902-07-7), no biological evaluation was ever carried out on it.
Table 5. Docking score of the five hits from direct docking of phytochemicals library along with
the crystal structure ligands
Direct docking hit/Crystal structure liganda
|
Docking score (PDB ID 6WUU)
|
PLpro residues involved in H-bonding with top-ranking pose of the ligand
|
Docking score (PDB ID 6WX4)
|
PLpro residues involved in H-bonding with top-ranking pose of the ligand
|
6WUU ligand
|
-10.64
|
Cys111, Asp164, Gly271, Tyr268
|
NDb
|
NAc
|
6WX4 ligand
|
NDb
|
NAc
|
-9.78
|
Gly163, Asp164, Cys270,
|
16
|
-5.93
|
Arg166
|
-6.34
|
Cys270
|
17
|
-6.73
|
None
|
-6.48
|
Asn109
|
18
|
-7.17
|
Cys111, Gly163
|
-7.54
|
Gly163, Gly271, Cys270
|
19
|
-8.44
|
Arg166, Cys270, Gly271
|
-6.30
|
None
|
20
|
-6.88
|
Arg166, Gln169
|
-5.85
|
None
|
aDocking score of the top-ranking pose; bNot determined; cNot applicable
Overall, the top-scoring hits scouted using the direct docking from the same phytochemical library were observed to be bigger and less drug-like, compared to the hits obtained using the complete process flow based on peptide mold concept. When the overlap of top-500 hits from both the peptide mold and direct docking strategies was investigated, we found less than 50 hits common both workflows (data not shown).