GRL0617 has been reported to inhibit the protease activity of PLpro from SARS-CoV-2 and crystallographic studies have determined the structural details supporting its activity[11]. Considering this promising background, other researchers have studied the activity of other compounds with chemical or structural similarity to GRL0617, obtaining candidates with higher affinity for PLpro[17]. However, groups of less than 500 thousand compounds have been explored and the selection has been made mainly based on physicochemical descriptors[23]. In the present work we performed a search for compounds structurally similar (Tanimoto Index 85%) to GRL0617 in the PubChem database containing ~ 109 million molecules.
Molecular Docking
To evaluate PLpro/ligand binding, molecular docking was used and the most stable conformation and ∆Gbind of the protein in complex with 17,889 compounds obtained from the first screening was determined. In addition, the same calculations were performed for the PLpro/GRL0617 (-10.0 kcal mol-1) and PLpro/12C (-8.8 kcal mol-1) complex, two compounds with proved in vitro activity and whose crystal structures have been published. Despite the large conformational space accessible to the protein/ligand complexes, the most stable conformations obtained by docking for PLpro/GRL0617 and PLpro/12C show a high similarity to their crystallographic structures. The structure of the best PLpro/GRL0617 docking is similar to the published (7CMD) crystal (rmsd of 2.4 Å) (Figure S1B). In the case of compound 12C, unlike what was observed with GRL0617, the PLpro/12C coupling calculation generates a different structure (Figure S2C) than the published (7E35), but at the same site and with a considerable ∆Gbind.
The energy ∆Gbind of the PLpro/GRL0617 complex (-10.0 kcal mol-1) was used as a threshold and 500 compounds with ∆Gbind between − 11 kcal/mol and − 14.1 kcal/mol were selected (Figure S2). In addition, the crystal structures of PLpro/GRL0617 and PLpro/C12 were used as controls allowing us a comparative analysis for each stage of the selection, mainly those associated with free energy calculations and MD simulation.
Physicochemical Descriptors and ADME Properties
The 500 compounds with structural similarity to GRL0617 and with the highest ∆Gbind were sorted according to descending LE to select the seven with the highest value of this physicochemical descriptor. According to this selection the order of the selected compounds from lowest to highest LE were D24, D28, D04, D59, D60, D99 and D06. By observing the chemical structures of these compounds (Figure S3), it is identified that molecules D24, D28, D59, D60, D99 and D06 present a naphthyl group as well as GRL0617 and 12C, while compound D04 presents a (9,10-dioxoanthracene-2-yl) amino]-1-oxopropan-2-yl group. This could indicate possible chemical bases for the affinity of these compounds for the PLpro binding site.
To evaluate the feasibility as a drug of each compound, six physicochemical properties were calculated: lipophilicity, size, polarity, solubility, flexibility, and saturation. As shown in Fig. 2, compounds D24, D28, D04 and D59 present properties within the range acceptable for a drug, as well as the control 12C. Similar to GRL0617, which exceeds one ADEM-Tox proprieties (unsaturation), 3 selected compounds exceed some parameters: compound D06 exceeds the flexibility range and compounds D99 and D60 exceeds flexibility and unsaturation parameter. As shown in Figure S4, all compounds, including the controls GRL0617 and 12C, show a high probability of gastrointestinal absorption, while the controls also show a high probability of brain penetration. In this theoretical framework the selected compounds could be drug candidates, however it is necessary to evaluate their affinity for PLpro using a methodology with higher predictive power. In this aspect, MD simulation allows to evaluate the protein-ligand interaction with a temporal and spatial resolution that is not possible to obtain so far through experimental studies.
Molecular Dynamics simulation and protein/ligand interactions
From the structures obtained by molecular docking, 200 ns MD simulations were performed in a TIP3P explicit water-box model. These simulations show a stable (Figure S5) PLpro/ligand interaction during the whole simulation time (Fig. 3A, 3B, 3D, 3E, 3F, 3G, 3H y 3I), except for the PLpro/D24 model that dissociates and binds to another region of the protein (Fig. 3C). The protein-ligand interaction is visualized in Fig. 4 (A to I), highlighting the residues that maintain a proximity of less than 3 Å for more than 80% of the simulation (red and orange). To identify residues that could participate in ligand stabilization, amino acids that remained within 3 Å (close residues) of the ligand for more than 80% of the simulation were identified (Fig. 4D to 4I). For compounds GRL0617 (Fig. 4A), 12C (Fig. 4B), D28 (Fig. 4D), D04 (Fig. 4F), D60 (Fig. 4G), D99 (Fig. 4H) and D06 (Fig. 4I) the nearby residues are Leu165, Asp167, Tyr267 and Gln272. Compounds D04, D59, D60, GRL0617 and 12C present the same nearby residues, which include Pro251, Tyr271 and Tyr276, in addition to those already mentioned. For the case of D24, the interaction with the protein is weaker, maintaining interactions in less than 40% of the simulation.
The ligand generates a local perturbation in the protein that can be observed by means of the root mean square fluctuation (RMSF). Figure 4I shows that compounds 12C (green line) and D28 (yellow line) generate higher fluctuation of the polypeptide chain compared to the protein without ligand (black dashed line). There are four regions of major perturbation in the protein, one near the N-terminus (residues 15 to 60), two near the binding site (135 to 175, and 180 to 240), which includes the Zinc finger domain and the compound binding region (255 to 300). To identify amino acids that can interact with the compounds, simulations were studied with the CPPTRAJ program Hbond, which uses geometric criteria to determine hydrogen bonds (Fig. 5). From the analysis it is observed that Asp 166 could act as a hydrogen acceptor for compounds D28 (Fig. 5D), D04 (Fig. 5F), D60 (Fig. 5G), D99 (Fig. 5H) and D06 (Fig. 5I), in addition to the control GRL0617 (Fig. 5A). Other residues of possible significance are: Tyr270 for GRL0617 (Fig. 5A); Gly165 and Gln271 for compound 12C (Fig. 5B); Thr77 and His75 for D24 (Fig. 5C); Tyr266, Tyr270, Gln271 and Tyr275 for D28 (Fig. 5D); Tyr275 for D04 (Fig. 5E); Tyr270 and Tyr275 for D59 (Fig. 5F); Thr303 and Asp304 for D60 (Fig. 5G) and Asp304 for D06 (Fig. 5I).
Non-covalent interactions
To obtain a representative conformation of the ligand/protein complex, a cluster analysis of 200 ns of simulation was performed with the DBSCAN algorithm. From the analysis, structures representing 84% of the simulation for the PLpro/12C complex and more than 91% of the simulation for the other complexes were obtained (Figure S6). These representative structures show the most relevant interactions for the stabilization of the PLpro/ligand complex. To characterize these interactions, NCI calculations were performed (Fig. 6), showing by means of an electron density profile, repulsive interactions in red, weak interactions in green and hydrogen bonds (HB) in blue. As can be seen, the NCI calculations show that weak (Van der Waalls) interactions between ligand and protein prevail, although some well localized stronger interactions can be identified. As seen in Fig. 6 (A, D-I), ligands can interact with the protein via HB. The nitrogen attached to the benzyl group of GRL6017 can form an HB with a hydrogen from the amide group of Gln272, and the carbonyl oxygen of this ligand can act as an acceptor for the hydroxyl proton of the side chain of Tyr267 (Fig. 6A). Two carbonyl oxygens of compound D28 can act as a hydrogen acceptor of the peptide bond of residue Glu164 (Fig. 6D). Asp167 can act as a proton acceptor forming a HB with a NH close to the cycloheptane of compound D04 (Fig. 6E). An HB can be established between the hydrogen of the peptide bond of residue Gln272 and a carbonyl oxygen of compound D59 (Fig. 6F). Residues Asp167, Tyr276 and Thr304 can act as proton acceptors to form HB with amino groups of compounds D60, in addition the peptide bond of residue Gln272 can act as a proton donor to form an HB with a carbonyl oxygen of this compound. The side chains of residues Asp167 and Tyr276 can form HB of compound D99 acting as proton acceptors. An amino group of compounds D06 can act as a proton donor to form an HB with the side chain of residue Asp167 and a carboxyl group as a proton acceptor to form an HB with the peptide bond of Gln272.
Free energy of binding by MMGBSA
It has been reported in vitro that compounds GRL0617 and 12C effectively inhibit PLpro protease (IC50 ~ 7 µM) and ubiquitinase (IC50 ~ 2 µM) activities. Based on this background we performed MMGBSA calculations to estimate ∆Gbind of the PLpro/ligand complex for each selected compound and of the PLpro/GRL0617 and PLpro/12C complexes (Table 1). The ∆Gbind allows estimating the affinity of a compound for a target protein, at more negative values the compounds are more affine for the protein. As seen in Fig. 7, the control compounds GRL0617 and 12C have similar affinity for PLpro, which correlates with the inhibitory activity measured in vitro. Compound D24 has less affinity for PLpro than the control compounds (GRL0617 and 12C) and compounds D28, D04 and D59 have a similar affinity as the controls. Compound D60 has a ∆Gbind ~6 kcal/mol lower than controls, in contrast D99 and D06 have a ∆Gbind ~10 kcal/mol lower than controls. Calculations of ∆Gbind performed using MMGBSA have smaller errors (Table S1) than the differences between compounds, supporting the observations. These differences are less evident in the ∆Gbin calculation performed using Molecular Docking (Fig. 7). For the complexes analyzed, Van der Waals interactions (ΔEvdw) generate the largest contribution to the total binding free energy (ΔGbind), whereas the electrostatic interaction (ΔEele) is counterbalanced by the polar desolvation energy (ΔGGB) (Table 1).
Table 1
Predicted binding free energies (kcal/mol) and individual energy terms calculated from molecular dynamics simulation following the MM-GBSA protocol for PLpro complexes.
| Calculated free energy of decomposition (kcal/mol) |
| ∆Gbind | ∆EvdW | ∆Eelect | ∆Ggas | ∆Gsolv |
GRL0617 | -32.6 | -37.9 | -21.6 | -59.52 | 26.9 |
12C | -32.8 | -38.9 | -22.1 | -61.0 | 28.2 |
D24 | -13.6 | -22.2 | -10.3 | -32.5 | 18.9 |
D28 | -30.2 | -38.2 | -31.9 | -70.1 | 39.8 |
D04 | -34.2 | -42.4 | -37,1 | -79.5 | 45.3 |
D59 | -30.5 | -42.2 | -25.5 | -67.7 | 37.2 |
D06 | -36.8 | -40.3 | -36.9 | -77.2 | 40.5 |
D60 | -42.0 | -44.7 | -42.2 | -86.9 | 44.9 |
D99 | -42.3 | -44.8 | -42.6 | -87.4 | 45.1 |
In this work we performed a massive search for compounds structurally like GRL0617, one of the few compounds capable of inhibiting the PLpro enzyme of SARS-CoV-2. By molecular docking we calculated ∆Gbind of ~ 18,000 compounds to PLpro which allowed us to select 500 compounds with ∆Gbind lower than GRL0617. After re-sorting these compounds based on LE we were able to select 7 compounds more affine by PLpro than GRL0617 with which molecular dynamics simulations were performed to evaluate PLpro/ligand stability and quantum mechanical calculations to identify their main interactions. These results show a high stability of the complexes, which can establish hydrogen bonds with the protein and maintain a high closeness. The free energy calculus has mean MMGBSA show that D60, D99 and D06 with mayor affinity for the protein.