Re-Exploring the Ability of Common Docking Programs to Correctly Reproduce the Binding Modes of Non-Covalent Inhibitors of SARS-CoV-2 Protease Mpro

In the latest few decades, molecular docking has imposed itself as one of the most used approaches for computational drug discovery. Several docking benchmarks have been published, comparing the performance of different algorithms in respect to a molecular target of interest, usually evaluating their ability in reproducing the experimental data, which, in most cases, comes from X-ray structures. In this study, we elucidated the variation of the performance of three docking algorithms, namely GOLD, Glide, and PLANTS, in replicating the coordinates of the crystallographic ligands of SARS-CoV-2 main protease (Mpro). Through the comparison of the data coming from docking experiments and the values derived from the calculation of the solvent exposure of the crystallographic ligands, we highlighted the importance of this last variable for docking performance. Indeed, we underlined how an increase in the percentage of the ligand surface exposed to the solvent in a crystallographic complex makes it harder for the docking algorithms to reproduce its conformation. We further validated our hypothesis through molecular dynamics simulations, showing that the less stable protein–ligand complexes (in terms of root-mean-square deviation and root-mean-square fluctuation) tend to be derived from the cases in which the solvent exposure of the ligand in the starting system is higher.


Introduction
In the 1980s, with the first study provided by Kuntz et. al [1], the computational technique of molecular docking had its birth. The efficiency, speed, and robustness of this method make its presence a constant in every structure-based drug-discovery pipeline [2]. To give a brief explanation, molecular docking consists of a multistep computational process that aims to find the best conformation of a molecule to bind to another to form a stable complex [3]. In the field of medicinal chemistry, as is deductible, its main application is finding the best molecules to bind in a firm way to the desired target (a protein, a nucleic acid, etc.). The algorithm starts with the exploration of the conformations space of the ligands (exploiting the so-called "search algorithm"). The conformations (called "poses") are then classified by a "scoring function", which attributes a numeric value to the goodness of the interaction according to energetical and/or geometrical function.
After a series of iterations, the final conformations are presented from the program to the user and ranked by the internal scoring function [4].
In the last 30 years, many docking programs have been developed. Among them, GOLD [5] (a genetic docking algorithm developed by the Cambridge Crystallographic Data Center-CCDC), Glide [6] (a systematic docking program released under license by Schrödinger), AutoDock [7] (a non-commercial genetic algorithm from The Scripps Research Institute), AutoDock VINA [8] (created by the same organization and released for non-commercial use), and PLANTS [9] (an algorithm based on an "Ant Colony Optimization" method) have gained popularity.
The performance of molecular docking programs can be measured by evaluating their ability to reproduce the experimental structural data, such as the crystallographic coordinates of a ligand into its binding site [10]. This ability has been evaluated in several benchmarks [11,12] to rank the performance of different algorithms regarding a specific target, usually using as the key parameter the root-mean-square deviation (RMSD) between the coordinates of the different poses given by the program and the crystallographic ones.
The ability to reproduce a crystallographic conformation strongly relies on different factors. First, the geometrical characteristics of the binding site, such as extension and shape, play a very important role; it is known that the performance of the algorithms has been improved to dock molecules in "cavities" or "pockets", rather than surfaces of proteins [13]. Second, the nature and the dimensions of the ligand are also crucial. Indeed, very small ligands may explore different places in a binding site, and the interactions that they can establish are usually few in number, reducing the "synergism" which could induce a molecule to keep a peculiar shape in a pocket [14]. On the other hand, even if drug-like molecules generally have higher conformational freedom, their dimensions force them to be oriented into a site in a more conserved way, so they have less roto-translational freedom.
In this study, we examined the ability of three docking programs characterized by diverse conformational sampling algorithms to efficiently reproduce the crystallographic pose of different ligands bound in different sites of a protein. To accomplish this task, a target in which several crystal structures were solved with the ligands located in different sites of the macromolecule itself was needed. To this scope, we considered a very recent and relevant target in the current pharmaceutical scenario, namely the SARS-CoV-2 main protease (M pro ).
In the last couple of years, with the pandemic spreading of the SARS-CoV-2 virus, the world of medical sciences had found itself fighting a new and dangerous adversary [15,16]. This biological entity, which is part of the coronavirus family, has been demonstrated to cause a pulmonary infection which eventually leads to serious complications, as witnessed by the high number of deaths that have already been linked to it (more than 5 million, at the present day [17]). The replication cycle of this virus strongly relies on the activity of its main protease (known as M pro or 3CL pro , a crystallographic structure example is reported in Figure 1) [18]. Indeed, this protein is responsible for the cleavage of the propeptide transcribed by the viral genome. In this way, the formation of all the functional proteins for the building of new virions takes place, and so the viral infection can proceed. Even if many molecules have been shown to bind to M pro [19] and inhibit its activity, and even if a molecule is currently in phase III clinical trial for this purpose (PF-07321332, from Pfizer [20,21]), no drug has already been approved by the European Medicinal Agency for the treatment of SARS-CoV-2 (also called "COVID-19"). Computational methods have already proven to be beneficial in the research for new potential inhibitors for M pro , as the literature witnesses [22,23]. In this work, we decided to implement a molecular-dockingbased approach relying on the programs GOLD, Glide, and PLANTS. These algorithms are considered "orthogonal" because they are characterized by diverse placing and scoring algorithms to obtain the best solution to the "protein-ligand posing problem". Each of these programs was used to dock each of the different non-covalent ligands of the various crystal structures of M pro , and this allowed us to evaluate the factors which influence the variability in reproducing the crystallographic poses. A self-docking protocol similar to the one herein reported had already been developed by our laboratory, with the name "DockBench". This program was implemented with success in several workflows, as the literature assesses [24,25]. In this study, a slightly modified version of that tool was used, which exploits only three docking programs at the present moment but can expand the analysis of the results obtained. Looking at the docking benchmarking protocols on M pro , we see that a remarkable study has already been conducted and published by Zev et al. [26]. In that specific work, six different docking programs were evaluated in their performance in reproducing the M pro non-covalent ligands' crystallographic poses, and three of those algorithms have also been compared in their ability to correctly place M pro covalent ligands into their proper binding site. In our work, we decided to expand the considerations brought by that study, evaluating specifically how docking performance changes in respect of the crystallographic data that have to be reproduced.
Indeed, we considered in our calculations parameters such as the solvent exposure of the ligand and the influence of the crystallographic water molecules in docking calculations, focusing our evaluations just on non-covalent M pro ligands. We executed the experiment in two different scenarios, one which excluded the crystallographic waters from the calculation (which we will name "Scenario 1"), and one which induced the docking programs to consider them (called "Scenario 2"). After that, we compared the docking results with the percentage of solvent exposure of the crystallographic pose of the ligand, successfully confirming that a higher solvent exposure tendentially reflects a worsening in the ability to reproduce the crystallographic pose by the algorithms (that, as already mentioned, are better trained for "cavities" rather than "surfaces"). To further investigate this aspect, we expanded our computational analysis by performing a molecular dynamics (MD) experiment, in which each crystallographic ligand was left free to move for 5ns (three replicas per system). This approach (known as "MD post-docking") has already become part of our computational protocol [27,28] and is based on the fact that the conformations of the ligands, which are less prone to be displaced from their initial position during the simulation, are related to higher stability and binding strength with the target. In the case presented, this principle was applied directly to the crystallographic conformations of the ligands rather than to docking poses. This was performed because the goal was not to select the most promising molecules in binding to a specific region of the protein; instead, the objective was to elucidate which are the features of the crystallographic ligands that tend to guarantee a tighter binding with the receptor. Our evaluation demonstrated that the molecules bound to the orthosteric pocket of M pro keep their position much stronger than the molecules crystallized on other sites, further validating our solvent exposure-based theory. A representation of the M pro ligands crystallized in the various sites of the protein is given in Figure 2.

Scenario 1-Docking Calculations without Considering the Crystallographic Water Molecules
The results of our docking protocol for this section (which are numerically reported in the Supplementary Materials File "Selfdocking_scenario1.csv") are graphically represented with colormaps. All the colormaps present in this study are based on a colorimetric scale delineating the RMSD values, starting from 0 Å, which corresponds to a molecular docking pose exactly super posable to the crystallographic one (maximum docking performance, represented by the dark blue color), and reaching values of 5 Å or higher (minimum docking performance, all represented by the dark red color), corresponding to a very low level of overlap between the coordinates of the pose produced and the ones of the crystallographic conformation. The colormaps in Figure 3 show the self-docking results obtained on the different M pro crystal structures in the case in which water molecules are not considered in the calculation. As is depicted, the RMSD values were far lower for all the complexes in which the crystallographic ligand is located in the orthosteric pocket. To give a better resolution of this, we separated each map into two different colormaps, one grouping all the 78 proteins in which the ligand is located into the catalytic pocket, and one including all other cases (41 complexes).
We analyzed the data coming from the calculations, and we determined that, looking at all the complexes with all the different couples docking program-scoring functions, we see that the average values of all the RMSDs obtained was 5.76Å ("RMSD_average"). Looking at the average of the RMSDs coming from the poses which were scored as the best ones from the algorithms' scoring functions ("RMSD_scor_func"), we see that the value was 5.10Å. If the lowest RMSD values only are taken into account for each docking run ("RMSD_sorted"), the mean of the values was 3.70Å.
The average values were also calculated separately for all the complexes in which the crystallographic ligand is located in the catalytic pocket, and for all other cases. The colormaps for these different conditions are reported in Figures 4 and 5.
First, the analysis focused on the complexes having the crystallographic ligand located within the orthosteric pocket. For this set of systems, we calculated the average RMSD value of all the poses ("RMSD_average"), which was revealed to be 4.54Å. Then we computed the average of the RMSD values coming from the poses which were ranked with the best score from the scoring functions ("RMSD_scor_func"), and its value was 3.43Å. Finally, the average RMSD value of the poses with the lowest RMSD in each run was calculated ("RMSD_sorted"), and its measure was 2.45Å.
Second, the same steps were performed for the rest of the protein-ligand complexes, which are the ones in which the crystallographic ligand is located outside the orthosteric binding site. Moreover, in this case, the first passage involved the calculation of the average RMSD value of all the poses generated for these systems ("RMSD_average"), and its measure was 8.08Å. Then, the mean of the RMSD values coming from the poses which received the highest rank from the scoring functions was calculated ("RMSD_scor_func") and was revealed to be 8.29Å. In the end, the average value of the lowest RMSDs of each run was computed ("RMSD_sorted"), and its measure was shown to be 6.08Å.  The results obtained for Scenario 1 are summarized in Table 1.

Scenario 2-Docking Calculations Considering the Crystallographic Water Molecules
The outcomes of our molecular docking experiment for this section (which are reported in the Supplementary Materials File "Selfdocking_scenario2.csv") are graphically represented with colormaps, which were created with the same criteria listed in the previous paragraph. The results reported in the colormaps in Figures 6-8 reveal the self-docking performance obtained on the different M pro crystal structures in the case in which the crystallographic water molecules within 5 Å from the ligand were retained during the calculation. Moreover, in this case, it is easy to notice that the values result in being far better for the complexes in which the small molecule of interest is in the orthosteric binding site.   Similar to the first scenario, we divided each colormap into two sets, one with the 78 proteins having the ligand located into the catalytic pocket, and the other including all the remaining cases (41 proteins). Considering all the protein-ligand complexes with all the different pairs docking program-scoring function, the mean values of all the RMSDs obtained ("RMSD_average") was 5.64Å, but focusing only on the mean of the RMSDs derived from the poses which were given the highest rank from the algorithms ("RMSD_scor_func"), the value resulted to be 4.83Å. Looking only at the best RMSDs for each docking run ("RMSD_sorted"), we see that the average of the values was 3.68 Å.
As already performed for Scenario 1, also in Scenario 2, the analysis was divided between the complexes in having the crystallographic ligand crystallized into the catalytic pocket, and for all other situations.
We reported the colormaps which resulted from this evaluation, and those are represented in Figures 7 and 8.
We started from the complexes in which the ligand is located inside the catalytic pocket in the crystal. For those systems, the mean of the RMSD values coming from all the poses("RMSD_average") resulted in being 4.22Å. Then, the average of the RMSDs derived from the scoring function highest-ranked poses in all the docking runs ("RMSD_scor_func") was computed, and its value was 3.11Å. In the end, also the average value between the lowest of the RMSDs in each docking run was calculated ("RMSD_sorted") and was revealed to be 2.26Å.
Second, we repeated the analysis for all the complexes in which the crystallographic ligand is located outside the orthosteric pocket. For these systems, the average of the RMSD coming from all the poses collected in the docking runs ("RMSD_average") was calculated to be 8.32Å. Next, we computed the mean of the RMSD values derived from the poses which received the highest score (from the scoring functions) in each run ("RMSD_scor_func"), and this value was 8.11Å. Last, also the average value between the lowest of the RMSDs in each docking run was calculated ("RMSD_sorted"), giving 6.36 Å.
The results obtained for Scenario 1 are summarized in Table 2.

RMSD_Average (Å) RMSD_Scor_Func (Å) RMSD_Sorted (Å)
All Just analyzing the numbers coming from the average values allows us to see how the performance of the docking programs dramatically increases when the ligand is docked inside the catalytic pocket rather than on the surface of the protein, in line with the fact that the molecules have a limitation in the conformation that they can explore into a binding site. Together with this, the ligands can exploit their accessible surface area to interact with the protein more efficiently, following the principle of "complementarity" [29,30].

Solvent Exposure Analysis
The results of the docking calculations were then analyzed in light of the data coming from the solvent exposure analysis. For each docking program-scoring function pair, the best RMSDs given by the docking calculation were evaluated against the solvent exposure of the ligand in its crystallographic pose. The results were reported in different plots, one for each couple docking program-scoring function, also in this case dividing the graphs in respect to the "scenario" from which the docking result was coming. To give an example, we reported in this article the plots for the pair GOLD-goldscore for each of these cases (Figures 9 and 10).  The plots arising from all other docking program-scoring function pairs, both in Scenario 1 and Scenario 2, are reported in Supplementary Materials Figure S1. From these graphs, we can easily see how the best RMSDs values tended to be derived from proteinligand complexes in which the solvent exposure of the ligand is low, and, most of the time, this means that the small molecule is crystallized in the binding pocket (indicated with the red dots in the plots). There are some cases in which the mean RMSD values were suboptimal also for this kind of ligands, and this can be due to several reasons. In some situations, of which the complexes 5REH, 5RE9, 5RGK (represented in Figure 11), and 7AVD are an example, the solvent exposure was tendentially higher in respect to the other orthosteric ligands, while, in other cases, the increase in RMSD can be attributable to the small dimensions of the ligand itself, making it harder for the docking algorithms to reproduce the reference pose in a pocket of such considerable volume (the complexes 5R82 and 5RG0 are an example for this) [31]. Figure 11. Representation of the crystallographic complex conformation of 5RGK, one of the proteinligand complexes in which the crystallographic ligand is located inside the orthosteric binding site, but the docking calculation results in high RMSD values. This is mainly due to the high level of solvent exposure that characterizes this ligand, which locates just a small portion of its structure inside the pocket, leaving the rest in an outer zone. The ligand is represented with stick representation (C-atom are colored in magenta), and the catalytic dyad (Cys145 and His41) is highlighted, as well as the His163 and the binding site residue interacting with the ligand. To give a better representation, the surface of the protein in a 5 Å radius from the ligand is represented and colored in blue.
On the other hand, there are also some cases in which the best RMSD given by the protocol was pretty low, even if the crystallographic ligand was not placed inside the orthosteric pocket. This is the case, for example, of 7LFP (the crystallographic pose is reported in Figure 12); the ligand was placed at the interface between the monomers, and so its solvent exposure and RMSDs values were low, even if was marked to be "outside the catalytic pocket". A similar situation is observed on 5RF0, where the ligand, even if not located into the orthosteric pocket, is not situated in the peripheral part of the protease. Figure 12. Representation of the crystallographic pose of 7LFP, which is one of the protein-ligand complexes in which, even if the crystallographic ligand is located outside the orthosteric binding site, the RMSD values between the original coordinates and the ones given from the docking runs are considerably low. The reason for this can be found in the very low solvent exposure of this ligand, which is located in the interface between the monomers, and so is shielded by them. The ligand is represented in orange, and the catalytic dyad (Cys145 and His41) of both monomers is highlighted. To give a better representation, the surface of the protein in a 5 Å radius from the ligand is represented and colored in blue.

Molecular Dynamics Simulations
For each of the 119 crystallographic complexes, three different molecular dynamics simulations (MD) of 5 ns each were collected to examine the behavior of the ligands in a dynamic environment. The trajectories were wrapped, aligned to the first frame and the root-mean-square fluctuation (RMSF) of the ligand, as well as the RMSD between its crystallographic and final coordinates ("RMSD_final"), and were calculated for every single experiment. For each protein, the values coming from the average of the RMSFs and "RMSD_final" derived from the replicas were considered. Considering all the simulations collected, the average of all the ligand RMSF values was calculated to be 5.28 Å, while the average of the RMSD values between the coordinates of the crystallographic conformation of the ligand and the ones coming from the last frame of the trajectory ("RMSD_final") was of 8.89 Å.
As already performed for the docking results analysis, we first focused on the complexes in which the crystallographic ligand is originally located inside the orthosteric pocket. For these systems, the average of all the RMSFs coming from the simulations was 2.19 Å. The mean value of the RMSDs of the ligands in the last frame of each trajectory ("RMSD_final") was instead calculated to be 4.43 Å.
Second, we concentrated on the systems in which the crystallographic position of the ligand (and so its initial location) is outside the catalytic pocket. For these systems, the average value of all the ligand RMSFs during the trajectories was calculated to be 11.15Å. Then the RMSD value between the final coordinates of the ligands and their crystallographic ones ("RMSD_final") were considered. The average of these values, for all the trajectories collected for these complexes, was 17.66Å. The output files of the molecular dynamics simulation geometric analysis are available in Supplementary Materials "MD_data.csv".
As already performed for the docking experiments, also for MD results, the average values of RMSF and "RMSD_final" were plotted against the percentage of solvent exposure of the crystallographic conformation of the ligand, and the plots that were obtained are reported in Figures 13 and 14.
As expected, the complexes in which the ligand is crystallized in the orthosteric site (marked with the red dots in the scatter plot) tended to fluctuate much less than the ligands which are complexed in the external parts of the protease (represented with the blue dots in the graphs). As depicted, MD analysis confirms that the best values in terms of RMSF and "RMSD_final", which are correlated to a more energetically stable situation for the protein-ligand complex, come from the systems in which the crystallographic ligand is localized inside the catalytic pocket and are characterized by a low percentage of solvent exposure. These outcomes further support the already-mentioned hypothesis about the correlation between the improvement of the docking performances in the case in which the binding site is a pocket rather than a surface.
The overall results obtained with molecular dynamics simulations are summarized in Table 3. A graphical representation of the molecular dynamics simulations is reported in Supplementary Materials "Video_S1.mp4". In this video, the ligands crystallized into the catalytic pocket are colored in magenta, while the other ligands are colored in cyan.

Software Overview
The molecular modeling operations were executed with the Molecular Operating Environment (MOE) suite (version 2019.01) [32]. The molecular docking calculations were carried out with CCDC GOLD (version 2020), Schrodinger Glide (from the Schrödinger suite 2021.3), and PLANTS. The solvent exposure calculation was performed with a series of SVL commands (exploiting "moebatch" of the MOE suite) implemented into a python script. The systems for molecular dynamics simulations were prepared by using tleap [33] and VMD [34]. The simulations were then carried out with ACEMD3 [35](version 3.3.0), a licensed program based upon OpenMM [36] (version 7.4.0). The modeling and docking calculations were performed on a 12 CPU (Intel Xeon E5-1620 3.50 GHz) Linux Workstation, while the MD simulations were carried out on a GPUs-cluster based composed of 20 NVIDIA GPUs.

Structure Preparation for Docking Calculations
The different crystal structures of M pro were collected from the Protein Data Bank [37]. Among these, the proteins which did not present any ligand, or which were complexed with a covalent ligand, were excluded. This way, only the non-covalent protein-ligand complexes were retained, and the complete list of all the 119 complexes used is available in Supplementary Material Table S1. The structures were grouped into a database and were prepared with MOE "QuickPrep" tool. With this tool, each complex was properly prepared to recreate the small missing loops in the structure, assigning the proper conformation to the residues with alternate orientation (based on occupancy) and adding the hydrogens to the system (this last passage was performed with the MOE "Protonate 3D" tool). The hydrogens added this way were then minimized by using the AMBER10:EHT force field implemented in MOE [38].
After these preliminary but crucial steps, each complex was manually examined and treated to eliminate every molecule, except for the crystallographic waters and the main ligand. Each complex was then independently saved.

Docking Calculations
For each of the complexes prepared, the crystallographic ligand was separated from the protein and self-docked into its binding site. For each docking program, all the scoring functions available were used for separate runs, and in each run, 5 poses were collected for the ligand. GOLD supports 4 different scoring functions: goldscore, chemscore, asp, and plp; Glide supports two main functions for docking, which are Glide-SP and Glide-XP, while PLANTS implements plp and chemplp.
For each docking-program-scoring function couple, the docking calculation was carried out in two different scenarios: one in which the crystallographic water molecules were not considered (which we refer to as Scenario 1) and one in which also the water molecules 5 Å or nearer from the ligand atoms were taken into account into the computation (which we refer to as Scenario 2).
When all the docking calculations were executed, the ligand root-mean-square deviation (RMSD) between the coordinates of each one of the poses and the crystallographic conformations were computed. The data of major interest were the RMSD in respect to the pose which is marked with the highest score by the program (RMSD_scor_func), the lowest RMSD of the docking run (RMSD_sorted), and the average of the RMSDs of all the poses generated (RMSD_average). The output files of the self-docking experiments executed are available in the Supplementary Materials ("Selfdocking_scenario1.csv" and "Selfdocking_scenario2.csv").

Solvent Exposure Calculation
For each M pro complex, the solvent exposure of the main crystallographic ligand was calculated with an SVL script based on MOE "moebatch". The output of such calculation was the percentage of the ligand surface which is exposed to the solvent in the protein-ligand crystallographic complex. All the percentages obtained are presented in Supplementary Materials Table S2.

Molecular Dynamics Simulations Setup and Execution
All the protein-ligand M pro systems were independently prepared for molecular dynamics simulations. The program tleap was used for the creation of the simulation box, which was set to be cubic and characterized by a 15 Å padding. The solvation model used was the explicit TIP3P, and the neutralization of the system was performed by adding Na + and Cl − ions until the salt concentration inside the box reached the value of 0.154 M.
The systems then underwent a two-passage equilibration. In the first one, both protein and ligand atoms were subjected to a harmonic position restrain of 5 kcal/mol. The length of this step was set to 0.1 ns, and the ensemble used was the canonical one (NVT). During the second equilibration step, the ensemble was set to NPT (isothermal-isobaric), the length was 0.5 ns, and the harmonic restrains (always 5 kcal/mol) were applied both on the ligand and on the alpha-carbons of the protein backbone.
After these preliminary steps, 3 replicas of 5 ns each were collected for each system; the ensemble was again the NVT one, and no restraints were applied. At the end of the simulations, the average root-mean-square fluctuation (RMSF) of the ligand during the trajectory, as well as the RMSD betweencrystallographic coordinates of the ligand and the ones coming from the last frame of the trajectory, were collected.

Conclusions
In this study, we evaluated the performance of three orthogonal docking algorithms in reproducing the crystallographic pose of several ligands located in different parts of the same target, which, in our case, was the SARS-CoV-2 main protease (M pro ). Our analysis revealed how, even if the docking programs used operate in different ways to give the final conformations to the user, all of them perform much better in the case in which the ligands are located in a binding pocket rather than crystallized outside of it. Specifically, we reported that their performance tends to decrease with the increment of the exposure of the crystallographic pose to the solvent. This was confirmed both from the experiments executed without considering the crystallographic water molecules in the docking calculations and from the ones taking into account the waters 5 Å or nearer to the ligand. Molecular dynamics simulations further give credit to our study, demonstrating how the less-fluctuating ligands (and so the most stable) through the trajectories were the ones crystallized inside the orthosteric binding site of M pro .