3.3. Docking and Binding Energy Calculations—PLpro

For docking to PLpro and calculations of binding energy, we used BIOVIA Discovery Studio v20.1.0.19295 [56]. The structures retrieved from PDB were prepared using the Prepare proteins protocol with pH set to 7.4 and the CHARMm forcefield. We used spherical gridboxes with a radius of 15 Å that were created around a ligand if the PLpro complexes featured one, or around the PDB ID: 3e9s ligand after superimposing the proteins by Cα atoms. We prepared the potential inhibitors using the Prepare ligands protocol and changed their ionization for pH 7.5 ± 1 without generating tautomers and isomers. Afterward, we employed the CDOCKER protocol, a tool performing grid-based docking to a rigid protein with ligand conformational sampling utilizing molecular dynamics. We retrieved only the best pose for each ligand by setting the option top hits to 1. After docking, we assessed the obtained poses by scoring functions in the Score ligand poses protocol and utilized the Calculate binding energies protocol with generalized Born implicit solvent model (MM–GBSA). Afterward, we created a MLR model by merging the Jain scoring function and the MM–GBSA binding energy. For this purpose, we utilized Statistica 13.1 software [57].

Initially, we validated the ability of the software to predict the correct inhibitor poses by employing the redocking and cross-docking techniques. We retrieved various PLpro structures (PDB IDs: 7jir, 7jit, 7jiv, 7jiw, 6wuu (all chains), 6wx4, 6w9c (chain A), 6xa9 (chain A), 7jn2, 7jrn (all chains), 7cjm, 7cmd, and 7cjd (all chains)) and redocked their original inhibitors as well as cross-docked inhibitors from other PLpro complexes. Additionally, we docked a few ligands found in complexes of SARS-CoV PLpro (PDB IDs: 3mj5, 4ovz, and 4ow0). The obtained protein–ligand complexes were superimposed by Cα atoms on the original PLpro complexes with the corresponding ligands in order to calculate ligand heavy-atom RMSD values. Based on the validation results we chose the PLpro structures with the lowest cross-docking RMSD values and used them to test the ability of the docking procedure to predict binding affinities of inhibitors. Two of the chosen structures (PDB IDs: 7jiv and 7jit) had a C111S mutation so we created additional structures where the serine residue was mutated back into cysteine. We created a set of 25 PLpro inhibitors with known IC₅₀ values (Supplementary Table S7). After docking them to the chosen PLpro structures, we assessed the obtained poses by scoring functions (CDOCKER Energy, CDOCKER Interaction Energy, LigScore1, LigScore2, PLP1, PLP2, Jain, PMF, PMF04) and MM–GBSA calculations. Subsequently, we checked Pearson correlation coefficients and the corresponding p-values between pIC₅₀s and the values obtained from scoring functions and binding energy calculations. If an inhibitor had its IC₅₀ measured for a racemate, we took an average of scoring functions and binding energies calculated for both enantiomers. Through this validation, we found the best PLpro structure, which was used for the subsequent drug screening. Additionally, we prepared a MLR model for this structure and the best-performing scoring function (Jain) and MM–GBSA. To validate the model, we examined p-values for the model itself as well as for all the terms of the equation. All p-values came below 0.01. Due to the simplicity of the model and to the occurrence of only two independent variables, we did not conduct cross-validation.

We conducted an additional step of validation in order to confirm the best-performing docking procedure’s (PDB ID: 7jn2, Jain, and MM–GBSA) ability to correctly predict potential inhibitors’ binding poses and affinities toward PLpro. First, we docked inhibitors cocrystallized with PLpro structures available in PDB that we did not utilize in the main cross-docking validation (PDB IDs: 7koj, 7kok, 7lbr, 7lbs, 7llf, 7llz, 7los, and 7e35). We superimposed the complexes and calculated the ligand’s heavy atom RMSD analogically as described above. Next, we created an additional set of PLpro inhibitors with known IC₅₀ values (Supplementary Table S8). We docked them to 7jn2 and conducted analogical scoring and binding energy calculations as described above. We determined Pearson correlation coefficients and the corresponding p-values between the experimental pIC₅₀ values and the values of Jain function, MM–GBSA binding energies, and pIC₅₀ values estimated using the MLR model for the combined two sets of test ligands.

In the screening part, we prepared the potential inhibitors and docked them using the same procedure as during validation, and then assessed the obtained poses by the Jain scoring function and MM–GBSA calculations.