Fragment dissolved molecular dynamics: a systematic and efficient method to locate binding sites

1 Fragment Dissolved Molecular Dynamics: A systematic and efficient method to locate binding sites Cristian Privat1, José M. Granadino-Roldán2, Jordi Bonet1, Maria Santos Tomas3, Juan J. Perez4, Jaime Rubio-Martinez*1 1 Departament de Ciència dels Materials i Química Física, Universitat de Barcelona (UB) and the Institut de Quimica Teorica i Computacional (IQTCUB), Martí i Franqués 1, 08028 Barcelona, Spain 2 Departamento de Química Física y Analítica, Facultad de Ciencias Experimentales, Universidad de Jaén, Campus “Las Lagunillas” s/n, 23071, Jaén, Spain 3Department of Architecture Technology, Universitat Politecnica de Catalunya, Av. Diagonal 649, 08028 Barcelona, Spain 4Deparment of Chemical Engineering, Universitat Politecnica de Catalunya, Av. Diagonal 647, 08028 Barcelona, Spain * Correspondence and requests for materials should be addressed to J.RM. (email: [email protected]) Page 1 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 2 ABSTRACT Diverse computational methods to support Fragment-based drug discovery (FBDD) are available in the literature. Despite their demonstrated efficacy to support FBDD campaigns, they exhibit some drawbacks such as protein denaturation or ligand aggregation that have not been yet clearly overcome in the framework of biomolecular simulations. In the present work, we discuss a systematic semi-automatic novel computational procedure, designed to surpass these difficulties. The method, named fragment dissolved Molecular Dynamics (fdMD) utilizes simulation boxes of solvated small fragments, adding a repulsive Lennard-Jones potential term to avoid aggregation, which can be easily used to solvate the targets of interest. This method has the advantage of solvating the target with a low number of ligands, thus preventing this way denaturation of the target, while simultaneously generating a database of ligandsolvated boxes that can be used in further studies. A number of scripts are made available to analyze the results and obtain the descriptors proposed as a means of trustfully discard spurious binding sites. To test our method, four Test cases of different complexity have been solvated with ligand boxes and four molecular dynamics runs of 200 ns length have been run for each system, which have been extended up to 1 µs when needed. The reported results point that the selected number of replicas are enough to identify the correct binding sites irrespective of the initial structure, even in the case of proteins having several close binding sites for the same ligand. We also propose a set of descriptors to analyze the results, among which, the average MMGBSA and the average KDEEP energies emerge as the most robust ones. Page 2 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 3 1. INTRODUCTION Fragment based drug discovery (FBDD) is a mature methodology currently used to accelerate the drug discovery process. With more than 40 drug candidates in the clinic and four marketed drugs including vemurafenib, venetoclax, erdafitinib and pexidartinib, FBDD was designed to circumvent shortcomings of High Throughput Screening (HTS) including limitations of the chemical space covered or its concomitant low hit rate1, 2. In addition, FBDD shows a lower attrition rate than HTS in clinical trials due the improved pharmacokinetic profile of the compounds tested. FBDD typically starts by screening a library of a few thousand of low molecular weight compounds –fragmentsthat exhibit specific features summarized in the “rule of three”3. Fragments exhibit low affinity -usually in the µM-mM rangeso screening requires the use of biophysical methods including NMR or X-ray crystallography initially and more recently, mass spectroscopy, thermal shift analysis or surface plasmon resonance, among others. Once fragments for a specific target are identified, several techniques can be used for lead generation using a combination of biophysical, biochemical and computational methods4, 5. Unfortunately, FBDD has also shortcomings including the occurrence of false positives or false negatives since detection techniques are pushed to their limits6 or sometimes, the lack of a detailed knowledge of the fragment-target interactions prevents the process to finish successfully7. To avoid these difficulties and at the same time being efficacious in FBDD campaigning, a combined use of the diverse technologies available is recommended. Computational methods, used either to prescreen libraries or to complement the diverse information provided by biophysical techniques can be helpful in FBDD campaigns8. In order to develop computational methods that can help the FBDD process it is necessary to understand the structural traits behind FBDD success. Analysis of ligandprotein complexes shows that binding pockets exhibit a small polar region where ligands are shown to establish two hydrogen bonds on average9. Interestingly, these hydrogen bonds are conserved among diverse ligands binding to the same pocket10. These regions are called ligand-binding hotspots and contribute to a large fraction of the corresponding binding free energy11. The first computer based methods for hotspot identification as GRID12 or SuperStar13 used atomic probes to identify areas of favorable Page 3 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 4 interaction on target protein surfaces. In contrast, more recent methods like FTMap14 use a set of small molecule probes, both hydrophobic or containing one or two polar functional groups to identify hotspots. This procedure was inspired in the experimental observation that when organic co-solvents are used to solve crystallographic structures they cluster around hotspots15. All these methods were designed to screen static crystallographic structures, although target flexibility can be incorporated by using the methods in a recursive fashion on an ensemble of protein conformations provided by a molecular dynamics (MD) study16, 17. Nowadays, the ever-increasing accessibility of computing power gives the opportunity to carry out a direct identification of hotspots. Inspired in the multiple solvent crystal structures (MSCS) method15, several procedures use the analysis of the solvent distribution produced in long MD trajectories of a protein in an organic cosolvent environment to identify hotspots18. Different strategies have been proposed in recent years. Thus, in MDmix19 organic co-solvent mapping is performed from the analysis of three 20 ns MD simulations of the macromolecule in water, ethanol/water and acetamide/water at 20%, respectively. The miscibility of the solvents avoids the need to introduce artificial forces to prevent phase separation or self-aggregation. In SILCS20 organic co-solvent mapping is performed with a solvent system including propane, benzene, methanol, formamide, acetaldehyde, methylammonium, and acetate. Here, a repulsive potential between fragments is used to prevent fragment association. Finally, in MixMD21 organic co-solvent mapping is performed using acetonitrile, isopropanol, and pyrimidine. Although these methodologies have demonstrated to be useful in FBDD, they have the drawback to identify a large number of extra hotspots that may be misleading in prospective applications. Moreover, these methods present limitations regarding the computation of the binding free energy since it depends on system setup and the concentration of probes used to perform the simulations18. To overcome these limitations, dynamic docking of fragments can be a plausible alternative. Diverse examples published in the literature demonstrate that long enough MD simulations allow performing ligand docking successfully exploring cryptic sites22-26. The purpose of the present study is to introduce a systematic and semi-automatic MD based procedure aimed at identifying hotspots and/or cryptic binding sites using Page 4 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 5 fragments without any previous knowledge of their affinity for the protein or having any description of their binding site. This method is designed to be used in early stages of a FBDD campaign, to prescreen a fragment library. The method, referred as fragment dissolved Molecular Dynamics (fdMD), consists of performing multiple 200 ns MD simulations of the target protein surrounded by many copies of the fragment of interest27. Specifically, simulation systems contain several copies of a fragment depending on their size and are constructed easily using pre-generated small boxes containing one copy of the fragment, soaked in an equilibrated box of TIP3P28 water molecules that are used as building blocks. These boxes will form part of a repository of fragment-solvated boxes that can be used for subsequent studies. Thus, solvation of the protein using this pre-equilibrated box renders simulation systems with diverse fragment molecules ranging from 6 to 26, avoiding the problem of protein denaturation. Finally, to understand the most suitable binding site as well as the best pose, discarding spurious binding sites, the method uses a set of descriptors obtained from the analysis of the simulations that permits to rank order the identified hotspots/cryptic sites. In the present report, the method is tested with four protein-fragment test cases of diverse complexity, using specific fragments for which there is experimental information available. The results suggest that the selected number of replicas are enough to identify the correct binding sites, irrespective of the initial structure used, even in the case of proteins having several close binding sites for the same ligand. 2. METHODS 2.1. fdMD procedure. The fdMD procedure is summarized in the flowchart shown in Figure 1 (a more detailed explanation is described in the SI). It consists of six steps that go from the preparation of the system to the analysis of the results derived from diverse computations. Technical details of the calculations are provided in sections 2.3 and 2.4. Steps 1 and 2 correspond to fragment and protein preparation, respectively. In step 1, the fragment of interest is soaked in a box of equilibrated TIP3P28 water molecules using the LEaP module of Amber18 (details are provided in section 2.4)29. After energy minimization, the system is subjected to a MD run (details are provided in section 2.4), and the final snapshot saved for step 3. Step 2 corresponds to protein preparation. In Page 5 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 6 step 3, the protein is solvated using several copies of the box containing the fragment produced in step 1. In step 4, the system is subjected to energy minimization and equilibration, followed by four 200 ns MD simulation generated using different initial velocities. Sampling is subsequently analyzed according to an automatized procedure involving steps 5-9 for each of the four MD simulations done. Thus, in step 5 water molecules and counterions are removed from the solvated trajectory and divided into N protein-ligand non-solvated trajectories each containing only the protein and one copy of the ligand, where N is the number of ligands in the simulation box. In addition, a pdb file of the last snapshot of each of the N trajectories is generated. This process is carried out using the script fdMD_OneTraj.cmd. In step 6, each of the N trajectories is analyzed using script fdMD_ReactiveTraj.py to assess if they have been reactive or not, in the sense that the ligand interacts with a prospective binding site during a pre-established simulation time. Subsequently, scripts fdMD_MMGBSA_Send.cmd and fdMD_MMGBSA_Analize.py are used to generate the inputs needed to run a MMGBSA calculation30 and to extract and graph data from those calculations, respectively. Finally, for each of the prospective sites occupied we compile information of the time the fragment remains bound (the residence time, RT) and the MMGBSA binding energy for the last 20 ns of trajectory. Moreover, the binding energy with KDEEP31 can also be computed using the last snapshot of each trajectory. KDEEP is a fast machine-learning approach for predicting binding affinities using 3D-convolutional neural networks. Analysis of the trajectories permits to identify prospective fragment binding sites. These sites are rank ordered using a consensus of a set of six descriptors including the number of reactive trajectories (i.e. how many of the four trajectories end up with one ligand bound to the prospective site); the longest RT; the average MMGBSA binding energy computed as the sum of the MMGBSA binding energy from the diverse trajectories divided by the total number of trajectories nt ( four in our simulations ) ; the average KDEEP binding energy computed as the sum of the KDEEP binding energy of the diverse trajectories divided by nt; the highest MMGBSA binding energy among the diverse reactive trajectories; the highest KDEEP binding energy among the diverse reactive trajectories. The six descriptors are used altogether and evaluated following a consensus approach, not as a scoring function. That is, the binding pose most highly “voted” is Page 6 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 7 selected as the most probable. If there is no interaction of the fragment with that binding site in a trajectory, the value of the descriptor for this trajectory is set to zero. After analyzing all prospective binding sites identified, the one that exhibits the best descriptors is taken as reference. Then, the descriptors for the remaining binding sites are compared to those of the reference. We propose as a consensus criterion to consider a binding site as false positive if more than three descriptors are worse than those of the reference binding site. In contrast, a binding site can be considered prospective if two or more descriptors are better than those of reference binding site. When more than one prospective binding sites are identified, the MD simulation time should be extended, the descriptors re-calculated, and the analysis repeated. It is important to note that there may be more than one experimental binding site. 2.2. Test cases. Four different protein-ligand test cases addressing diverse aspects of the methodology were selected to illustrate the robustness of the fdMD procedure (see Table 1). Test case I was designed to check whether the fdMD methodology is robust enough to identify the correct binding site for a fragment, irrespective of the initial structure of the target used for the study. For this purpose, we selected 4 fragments out of the 38 whose crystallographic structure bound to endothiapepsin was recently reported by N. Radeva et al.32. In the present study, we studied the binding of these four fragments to endothiapepsin using two different holo crystallographic structures as initial structure, and compared the results with the corresponding crystallographic structures of the complex ligand-protein32. Endothiapepsin exhibits a large, solvent exposed pocket that can be compartmentalized into nine sub-pockets that can be labeled according to the Schechter and Berger nomenclature33. Specifically, we chose fragments f031, f035, f207 and f240 whose chemical structures are shown in Table 1, spanned over the protein pocket. As mentioned above, two different initial protein structures were used for each ligand. One corresponds to the experimental structure of the ligand bound to that protein (labelled with the name of the ligand), while the other one corresponds to the protein structure in its holo form bound to the randomly selected f278 ligand (PDB code 5DR1)32 and labelled as NameLigand_5DR1. Page 7 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 8 Test case II was designed to evaluate the efficiency of the fdMD methodology in three different proteins that exhibit a much less solvent exposed pocket than that of endothiapepsin, using their crystallographic structures in the holo form (see Table 1). The systems were selected from a set of 12 protein-fragment complexes reported by M. Congreve et al.34. The set was also used in a theoretical article that assayed a combined Multiple Copy Simultaneous Search (MCSS) and MMGBSA approach35, avoiding the existence of cofactors in the structures. Test case III deals with the typical situation found in a FBDD campaign, where the available structure of a protein is the apo form. Moreover, test case III also investigates the performance of the fdMD methodology to identify the right binding site for larger ligands. Actually, one of the drawbacks of the fdMD methodology comes from the fact that the repulsion term used to prevent aggregation of ligands limits the size of the ligand that can be studied, because the cut-off distance for the van der Waals interactions is not large enough. To overcome this drawback, in the present work the ligand was split into two parts. Thus, test case III is actually composed of two ligands: 3chlorobenzo[b]thiophene-2-carboxylate (class1) and 3-(4-chloro-3,5dimethylphenoxy)propanoate (class2), which are known to bind myeloid cell leukemia 1 (MCL-1) with Ki of 131 µM and 60 µM, respectively. No experimental structure of any of the two fragment/MCL-1 complex is available. However, when linked together, the two fragments render a compound with Ki = 0.32 µM for which a crystallographic structure is available36 (PDB code 4HW3). Despite the availability of this structure, the structure of the target used in the present study is the apo form of the MCL-1 protein with PDB code 4OQ537. Finally, test case IV is used as a benchmark to compare the performance of fdMD on three systems previously studied with the ColDock method26. The starting structure used for these systems corresponds to the apo form of the target protein (see Table 1), even though their corresponding holo structures are available. One of these three fragments is dimethylsulphoxide (dmso), highly miscible in water and with a Ki = 20 mM to the FK506 binding protein (FKBP)38. The other two fragments, methyl sulphinylmethyl sulphoxide (dss) and -aminocapric acid (aca), show Ki = 250 µM with the FKBP 𝜀 Page 8 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 9 protein [40] and IC50 = 40-105 µM with the human plasminogen kringle 4 (HPK4)39, respectively. 2.3. Protein Preparation. The 3D structure of each protein was downloaded from the Protein Data Bank40 (PDB codes are shown in Table 1). Water molecules and organic molecules were removed from the structures, disulfide bonds were generated when required, and residues were protonated with the Protein Preparation Wizard41 of Maestro 2016-242 at the experimental pH, when available. Missing residues in the 2OHK PDB file were added by homology modelling, using the crystallographic structure of memapsin 243 (PDB code 1FKN) as template, by means of the Prime module44, 45 of Maestro 2016-2. Each protein structure was then loaded into the LEaP module of Amber1829, where counterions, if needed, were added, and each protein was solvated using TIP3P water molecules28 with a minimum distance from the edge of the box of 15 Å, and removing those water molecules closer than 2.0 Å from any protein atom. The ff14SB46 force field was used and the solvated protein was minimized with pmemd47, 48 for 10000 steepest descent (SD) steps using a cut-off of 10 Å for non-bonded interactions, being electrostatic interactions modelled with the particle-mesh Ewald method49. The apo structure of MCL-1 used for test case III was generated using a specific procedure. Specifically, starting with the crystallographic structure (PDB code 4OQ5), which was prepared in the same way as explained above, the minimized, ligand-free structure, was subjected to a heating and equilibration protocol using the CUDA version of pmemd47, 48. All bonds involving hydrogen atoms were constrained to their equilibrium value using the SHAKE algorithm50. A two steps preparation protocol with an integration time step of 2 fs was used which included increasing the temperature of the system to 300 K at a constant rate of 1.5 K ps-1 during 200 ps in the NVT ensemble and 500 ps at 300 K in the NPT ensemble to adjust density. After 200 ns of MD 300 K in the NVT ensemble, a clustering analysis with the cpptraj module embedded in Amber1829, allowing to extract the representative of the most populated cluster as an apo structure for MCL-1. 2.4. Ligand BOX Preparation. Each ligand was prepared with Maestro 2016-242 and optimized with semiempirical Austin Model 1 (AM1)51. Parameters and charges for the Page 9 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 16 seen in Tables 5 and S9, none of the four MD trajectories of 200 ns length was reactive with our proposed criteria. This is because the binding energy of dmso to the FKBP protein is very small and has a clear dynamic binding. The ligand interacts with a particular binding site of the receptor, remains in this site for a short period of time and goes out to another protein region, and so on continuously. Thus, a statistical approach which counts the number of times that a dmso molecule interacts with a specific point can be used to detect druggable pockets. This is the philosophy of the cosolvent approach and also the one used by ColDock. Our method does not use a statistical approach, but it requests that the ligand has a reasonable binding energy and thus keeps bonded in a specific binding site for a long period of time. It is important to note that all our four MD runs detect the binding of some dmso molecules to the experimental binding site, but always leave it (see Figure S3). The results for the other two systems are described in Table 6 (see also Table S9). As it can be seen, for both ligands the method detects only the experimental binding site with good agreement with the corresponding binding pose (Figure 5). 4. CONCLUSIONS In the present work, we introduce a computational semi-automated method designed as an aid in FBDD campaigns. Our method, named fragment dissolved Molecular Dynamics (fdMD), is based on the generation of simulation boxes of solvated small fragments that can be easily used to solvate around a target. It is worth to stress that the procedure of obtaining a solvated ligand box is done once and that box can be used later, in an easy way, to solvate any system and to prepare a database of ligand boxes to be used in a fragment-based drug discovery project. fdMD requires a low number of ligands, avoiding this way the risk of denaturation, and modifies the van der Waals interaction of a selected atom within the ligand using a pure repulsion potential to prevent ligand aggregation. Besides, we propose the use of six descriptors as a robust way to discard false binding sites. This work essays the performance of this methodology using four MD runs of 200 ns, which are extended up to 1 µs when needed. Page 16 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 17 After studying four test cases, which contain a total of seven protein systems with different types of binding sites, we concluded that our proposed method finds at least one experimental binding site, which usually corresponds to the most energetically favorable one according to MMGBSA and KDEEP results. Even in the case of ligands with several close binding sites, fdMD is able to recognize them despite the use of a repulsion potential for the ligands as several MD trajectories are used. Moreover, our results suggest that the effectiveness of the method is not dependent of the initial structure selected for the receptor, which is the usual case in a FBDD project. On the other hand, the number of replicas becomes crucial to discern between correct and false ligandprotein binding sites because it has a direct effect on the statistical average of the descriptors, being 4 MD runs reasonable. The average binding free energies estimated with MMGBSA or KDEEP approaches stand out as reliable descriptors to discard false binding sites, although we must stress that the computational cost of KDEEP is much smaller. Finally, our results suggest that MD runs of 200 ns length are a trustworthy starting point to deal with both exposed and unexposed binding sites remerging as a new alternative for the identification and optimization of potential fragments for the discovery of novel drugs. Some aspects mentioned in the discussion are on the table to improve the fdMD method. Future research will focus on the sampling of cryptic pockets or the identification of the correct binding modes by using enhanced-sampling techniques. 5. CONFLICTS OF INTEREST There are no conflicts of interest to declare 6. ACKNOWLEDGEMENTS This study was supported by The Agència de Gestió d’Ajuts Universitaris i de Recerca (AGAUR)-Generalitat de Catalunya (2017SGR1033) and the Maria de Maeztu MdM-2017-0767 program. José M. Granadino-Roldán thanks the Consejería de Economía y Conocimiento, Junta de Andalucía (FQM-337), and the Universidad de Jaén (Acción 1) for financial support. Page 17 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 18 Supporting Information Available: Detailed description of steps 1 and 3 of the fdMD flowchart. Tables S1 to S8 and Figures S1 to S3. The scripts used, together with a toy example, and a video showing the ligand for the f035 system oscillating between S3 and S6 pockets, can be downloaded at https://github.com/DrugDesignUBUJA/fdMD. Page 18 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 19 Table 1. Protein-ligand systems and number of ligands of each fdMD simulation. Liganda Protein Label PDB ID Ligandsb f031 5DPZ32 f031_5DR1 5DR132 9 f035 4Y3W32 f035_5DR1 5DR132 12 f207 4Y3T32 f207_5DR1 5DR132 12 f240 4YD632 8 TEST CASE I Endothiapepsin f240_5DR1 5DR132 12 Urokinase 1fv9 1FV967 26 Hsp90 2jjc 2JJC34 25 TEST CASE II BACE-1 2ohk 2OHK68 25 class1 20 TEST CASE III MCL-1 class2 4OQ537 26 dmso 8 FKBP dss 1D6O38 6 TEST CASE IV HPK4 aca 1PMK69 8 Page 19 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 20 a The dashed line circle indicates the C99/S99 atom (see below) b Number of ligands in the simulation box. Page 20 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 21 Table 2. Descriptors calculated for reactive trajectories of test case I. Pockets are defined in reference 32. MMGBSA and KDEEP energies in kcal mol-1. Systema Binding site Ave. MMGBSAb Ave. KDEEPc Best MMGBSA Best KDEEP N. react.d Best RTe Pocket S1´ -12.2 -6.0 -15.1 -6.4 4 180 f031 Pocket S6 -5.2 -3.5 -10.6 -7.2 2 75 Pocket S3 -6.1 -1.4 -24.4 -5.5 1 120 f035 Pocket S6 -2.9 -1.1 -11.5 -4.2 1 20 Dyadf -10.3 -4.3 -18.9 -6.1 3 180 f207 Pocket S6 -5.5 -2.5 -12.2 -5.2 2 80 Pocket S1 -16.9 -2.9 -18.3 -5.9 2 140 Natural Receptor f240 Pocket S6 -2.2 -1.1 -8.7 -4.4 1 40 Pocket S1´ -11.9 -5.2 -12.4 -6.2 4 180 f031_5DR1 Pocket S6 -2.7 -1.6 -10.7 -6.4 1 20 Pocket S3 -22.0 -4.3 -31.6 -6.3 3 170 f035_5DR1 Pocket S6 -3.9 -1.3 -15.7 -5.2 1 180 Dyadf -15.3 -5.3 -18.5 -5.9 4 170 f207_5DR1 Pocket S6 -2.9 -1.4 -11.4 -6.2 1 40 5DR1 Receptor f240_5DR1 Pocket S1 -15.9 -6.2 -17.3 -6.4 4 150 a See definition in Table 1.b Average MMGBSA binding energy (see Equation 1). c Average KDEEP binding energy. d Number of reactive trajectories interacting with that binding pocket. e In ns. f Catalytic dyad. In bold the best value for each descriptor. Page 21 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 22 Table 3. Descriptors calculated for reactive trajectories of set II. Pockets are defined in Table S6 and reference 57. MMGBSA and KDEEP energies in kcal mol-1. Simulation time Systema Binding site Ave. MMGBSAb Ave. KDEEPc Best MMGBSA Best KDEEP N. react.d Best RTe Experimental -11.6 -3.7 -17.8 -5.2 3 190 BS 1 -12.0 -2.7 -25.4 -5.7 2 190 BS 2 -3.1 -1.4 -12.3 -5.6 1 40 BS 3 -8.5 -2.4 -17.0 -4.9 2 138 BS 4 -3.0 -1.3 -12.1 -5.0 1 30 BS 5 -3.3 -1.3 -13.1 -5.0 1 152 BS 6 -4.0 -0.8 -16.0 -3.2 1 49 BS 7 -2.7 -1.4 -10.9 -5.7 1 54 200 ns BS 8 -2.8 -1.5 -11.2 -5.8 1 38 Experimental -15.6 -3.9 -25.4 -5.5 3 290 BS 1 -5.7 -1.4 -22.7 -5.6 1 290 300 ns 1fv9 BS 3 -8.4 -2.8 -16.6 -5.6 2 238 Experimental -6.9 -3.6 -11.0 -4.9 3 190 200 ns 2jjc BS 1 -4.6 -1.0 -18.2 -4.0 1 170 S1 -5.85 -2.6 -11.8 -5.2 2 161 BS 1 -2.48 -1.1 -9.9 -4.4 1 39 BS 2 -2.00 -0.9 -8.0 -3.6 1 20 BS 3 -3.60 -0.9 -14.4 -3.6 1 180 BS 4 -9.28 -3.6 -13.8 -5.6 3 141 BS 5 -1.68 -1.4 -6.7 -5.6 1 62 200 ns BS 6 -4.35 -2.3 -8.8 -5.5 2 47 S1 -5.3 -2.8 -11.0 -5.7 2 261 BS 3 -4.1 -1.1 -16.4 -4.5 1 280 300 ns BS 4 -9.6 -3.7 -12.8 -5.6 3 241 Experimental -2.5 -1.5 -9.9 -5.9 1 28 1 µs 2ohk BS 4 -3.3 -1.1 -13.2 -4.2 1 829 a See definition in Table 1.b Average MMGBSA binding energy (see Equation 1). c Average KDEEP binding energy. d Number of reactive trajectories interacting with that binding pocket. e In ns. In bold the best value for each descriptor. Page 22 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 23 Table 4. Descriptors calculated for reactive trajectories of set III. Pockets are defined in Table S6. MMGBSA and KDEEP energies in kcal mol-1. Simulation time Systema Binding site Ave. MMGBSAb Ave. KDEEPc Best MMGBSA Best KDEEP N. react.d Best RTe Experimental -21.1 -5.8 -22.3 -6.2 4 195 BS 1 -2.4 -1.1 -9.7 -4.5 1 25 BS 2 -7.8 -2.5 -18.7 -5.2 2 60 BS 3 -6.9 -2.1 -14.1 -4.2 2 50 BS 4 -16.4 -3.6 -17.4 -3.8 4 140 BS 5 -7.7 -1.9 -17.6 -4.1 2 60 BS 6 -6.0 -1.2 -24.0 -4.9 1 100 200 ns BS 7 -3.4 -1.4 -13.4 -5.5 1 50 Experimental -23.1 -5.6 -24.3 -6.4 4 295 BS 2 -3.1 -1.3 -12.4 -5.0 1 160 BS 4 -17.7 -3.6 -21.3 -4.0 4 260 300 ns class1 BS 6 -7.0 -1.0 -28.1 -3.8 1 200 Experimental -19.5 -4.6 -28.8 -6.3 3 175 BS 4 -11.9 -1.8 -24.9 -3.8 2 33 BS 5 -8.2 -1.6 -19.9 -3.2 2 135 BS 6 -11.6 -2.1 -31.7 -4.8 2 192 BS 8 -3.7 -1.1 -14.6 -4.3 1 60 BS 9 -3.7 -1.0 -14.6 -4.1 1 25 BS 10 -9.0 -2.3 -20.2 -4.7 2 110 BS 11 -5.2 -0.7 -20.6 -2.9 1 70 BS 12 -7.1 -1.3 -28.3 -5.0 1 110 BS 13 -4.2 -1.2 -16.6 -4.6 1 20 BS 14 -3.8 -0.9 -15.3 -3.7 1 122 200 ns BS 15 -3.1 -1.0 -12.4 -3.8 1 30 Experimental -18.5 -3.9 -29.8 -6.0 3 275 BS 4 -6.0 -1.0 -24.0 -3.9 1 133 BS 6 -7.3 -1.0 -29.0 -3.9 1 292 BS 10 -4.2 -1.0 -16.6 -4.1 1 178 BS 11 -5.1 -0.8 -20.2 -3.2 1 170 300 ns BS 14 -4.2 -0.9 -16.9 -3.6 1 222 1 µs class2 Experimental -15.9 -4.6 -23.2 -6.5 3 975 Page 23 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 24 a See definition in Table 1. b Average MMGBSA binding energy (see Equation 1). c Average KDEEP binding energy. d Number of reactive trajectories interacting with that binding pocket. e In ns. In bold the best value for each descriptor. Page 24 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 25 Table 5. Descriptors calculated for reactive trajectories of set IV. MMGBSA and KDEEP energies in kcal mol-1. Systema Binding site Ave. MMGBSAb Ave. KDEEPc Best MMGBSA Best KDEEP N. react.d Best RTe dmso dss Experimental -7.2 -6.0 -14.4 -6.0 2 155 aca Experimental -9.2 -5.3 -13.2 -5.9 4 23 a See definition in Table 1.b Average MMGBSA binding energy (see Equation 1). c Average KDEEP binding energy. d Number of reactive trajectories interacting with that binding pocket. e In ns. In bold the best value for each descriptor. Page 25 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 32 25. R. O. Dror, A. C. Pan, D. H. Arlow, D. W. Borhani, P. Maragakis, Y. Shan, H. Xu and D. E. Shaw, Proceedings of the National Academy of Sciences of the United States of America, 2011, 108, 13118-13123. 26. K. Takemura, C. Sato and A. Kitao, The Journal of Physical Chemistry B, 2018, 122, 7191-7200. 27. J. J. Perez, M. S. Tomas and J. Rubio-Martinez, Journal of Chemical Information and Modeling, 2016, 56, 1950-1962. 28. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, Journal of Chemical Physics, 1983, 79, 926-935. 29. D. Case, I. Ben-Shalom, S. Brozell, D. S. Cerutti, T. E. Cheatham III, V. W. D. Cruzeiro, T. Darden, R. E. Duke, D. Ghoreishi, H. Gohlke, A. W. Goetz, D. Green, R. Harris, N. Homeyer, S. Izadi, A. Kovalenko, T. Kurtzman, T. S. Lee, S. LeGrand, P. Li, C. Lin, R. Luo, B. Madej, K. M. Mermelstein, K. M. Merz, Y. Miao, G. Monard, H. Nguyen, H. T. Nguyen, I. Omelyan, A. Onufriev, D. R. Roe, A. Roitberg, C. Sagui, C. Simmerling, J. Smith, R. Salomon-Ferrer, J. Swails, R. C. Walker, J. Wang, R. M. Wolf, X. Wu, L. Xiao and P. A. Kollman, Amber 2018. 30. B. R. Miller, T. D. McGee, J. M. Swails, N. Homeyer, H. Gohlke and A. E. Roitberg, Journal of chemical theory and computation, 2012, 8, 3314-3321. 31. J. Jiménez, M. Škalič, G. Martínez-Rosell and G. De Fabritiis, Journal of Chemical Information and Modeling, 2018, 58, 287-296. 32. N. Radeva, S. G. Krimmer, M. Stieler, K. Fu, X. Wang, F. R. Ehrmann, A. Metz, F. U. Huschmann, M. S. Weiss, U. Mueller, J. Schiebel, A. Heine and G. Klebe, Journal of Medicinal Chemistry, 2016, 59, 7561-7575. 33. I. Schechter and A. Berger, Biochem Biophys Res Commun, 1967, 27, 157-162. 34. M. Congreve, G. Chessari, D. Tisi and A. J. Woodhead, Journal of Medicinal Chemistry, 2008, 51, 3661-3680. 35. M. K. Haider, H.-O. Bertrand and R. E. Hubbard, Journal of Chemical Information and Modeling, 2011, 51, 1092-1105. 36. A. Friberg, D. Vigil, B. Zhao, R. N. Daniels, J. P. Burke, P. M. GarciaBarrantes, D. Camper, B. A. Chauder, T. Lee, E. T. Olejniczak and S. W. Fesik, Journal of Medicinal Chemistry, 2013, 56, 15-30. 37. A. M. Petros, S. L. Swann, D. Song, K. Swinger, C. Park, H. Zhang, M. D. Wendt, A. R. Kunzer, A. J. Souers and C. Sun, Bioorganic & Medicinal Chemistry Letters, 2014, 24, 1484-1488. 38. P. Burkhard, P. Taylor and M. D. Walkinshaw, J Mol Biol, 2000, 295, 953-962. 39. L. Cheng, D. Pettersen, B. Ohlsson, P. Schell, M. Karle, E. Evertsson, S. Pahlen, M. Jonforsen, A. T. Plowright, J. Bostrom, T. Fex, A. Thelin, C. Hilgendorf, Y. Xue, G. Wahlund, W. Lindberg, L. O. Larsson and D. Gustafsson, ACS Med Chem Lett, 2014, 5, 538-543. 40. H. M. Berman, J. Westbrook, Z. Feng, G. Gilliland, T. N. Bhat, H. Weissig, I. N. Shindyalov and P. E. Bourne, Nucleic Acids Research, 2000, 28, 235-242. 41. G. Madhavi Sastry, M. Adzhigirey, T. Day, R. Annabhimoju and W. Sherman, Journal of Computer-Aided Molecular Design, 2013, 27, 221-234. 42. Schrödinger Suite 2016-2. Schrödinger, LLC, New York, 2016. 43. L. Hong, G. Koelsch, X. Lin, S. Wu, S. Terzyan, A. K. Ghosh, X. C. Zhang and J. Tang, Science, 2000, 290, 150. 44. M. P. Jacobson, D. L. Pincus, C. S. Rapp, T. J. Day, B. Honig, D. E. Shaw and R. A. Friesner, Proteins, 2004, 55, 351-367. Page 32 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 33 45. M. P. Jacobson, R. A. Friesner, Z. Xiang and B. Honig, Journal of Molecular Biology, 2002, 320, 597-608. 46. J. A. Maier, C. Martinez, K. Kasavajhala, L. Wickstrom, K. E. Hauser and C. Simmerling, Journal of Chemical Theory Computation, 2015, 11, 3696-3713. 47. A. W. Götz, M. J. Williamson, D. Xu, D. Poole, S. Le Grand and R. C. Walker, Journal of Chemical Theory and Computation, 2012, 8, 1542-1555. 48. R. Salomon-Ferrer, A. W. Gotz, D. Poole, S. Le Grand and R. C. Walker, Journal of Chemical Theory and Computation, 2013, 9, 3878-3888. 49. T. Darden, D. York and L. Pedersen, The Journal of Chemical Physics, 1993, 98, 10089-10092. 50. J.-P. Ryckaert, G. Ciccotti and H. J. C. Berendsen, Journal of Computational Physics, 1977, 23, 327-341. 51. M. J. S. Dewar, E. G. Zoebisch, E. F. Healy and J. J. P. Stewart, Journal of the American Chemical Society, 1985, 107, 3902-3909. 52. J. Wang, R. M. Wolf, J. W. Caldwell, P. A. Kollman and D. A. Case, Journal of Computational Chemistry, 2004, 25, 1157-1174. 53. A. Jakalian, B. L. Bush, D. B. Jack and C. I. Bayly, Journal of Computational Chemistry, 2000, 21, 132-146. 54. A. Jakalian, D. B. Jack and C. I. Bayly, Journal of Computational Chemistry, 2002, 23, 1623-1641. 55. C. W. Hopkins, S. Le Grand, R. C. Walker and A. E. Roitberg, Journal of Chemical Theory and Computation, 2015, 11, 1864-1874. 56. A. Onufriev, D. Bashford and D. A. Case, Proteins: Structure, Function & Bioinformatics, 2004, 55, 383-394. 57. P. Jain, P. K. Wadhwa, S. Rohilla and H. R. Jadhav, Bioorganic & Medicinal Chemistry Letters, 2016, 26, 33-37. 58. U. Neumann, M. Ufer, L. H. Jacobson, M.-L. Rouzade-Dominguez, G. Huledal, C. Kolly, R. M. Lüönd, R. Machauer, S. J. Veenstra, K. Hurth, H. Rueeger, M. Tintelnot-Blomley, M. Staufenbiel, D. R. Shimshek, L. Perrot, W. Frieauff, V. Dubost, H. Schiller, B. Vogg, K. Beltz, A. Avrameas, S. Kretz, N. Pezous, J.-M. Rondeau, N. Beckmann, A. Hartmann, S. Vormfelde, O. J. David, B. Galli, R. Ramos, A. Graf and C. Lopez Lopez, EMBO molecular medicine, 2018, 10, e9316. 59. K. Nakahara, K. Fuchino, K. Komano, N. Asada, G. Tadano, T. Hasegawa, T. Yamamoto, Y. Sako, M. Ogawa, C. Unemura, M. Hosono, H. Ito, G. Sakaguchi, S. Ando, S. Ohnishi, Y. Kido, T. Fukushima, D. Dhuyvetter, H. Borghys, H. J. M. Gijsen, Y. Yamano, Y. Iso and K.-i. Kusakabe, Journal of Medicinal Chemistry, 2018, 61, 5525-5546. 60. S. J. Veenstra, H. Rueeger, M. Voegtle, R. Lueoend, P. Holzer, K. Hurth, M. Tintelnot-Blomley, M. Frederiksen, J.-M. Rondeau, L. Jacobson, M. Staufenbiel, U. Neumann and R. Machauer, Bioorganic & Medicinal Chemistry Letters, 2018, 28, 2195-2200. 61. J. D. Low, M. D. Bartberger, K. Chen, Y. Cheng, M. R. Fielden, V. Gore, D. Hickman, Q. Liu, E. Allen Sickmier, H. M. Vargas, J. Werner, R. D. White, D. A. Whittington, S. Wood and A. E. Minatti, MedChemComm, 2017, 8, 11961206. 62. K. Fuchino, Y. Mitsuoka, M. Masui, N. Kurose, S. Yoshida, K. Komano, T. Yamamoto, M. Ogawa, C. Unemura, M. Hosono, H. Ito, G. Sakaguchi, S. Ando, S. Ohnishi, Y. Kido, T. Fukushima, H. Miyajima, S. Hiroyama, K. Koyabu, D. Page 33 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 34 Dhuyvetter, H. Borghys, H. J. M. Gijsen, Y. Yamano, Y. Iso and K.-i. Kusakabe, Journal of Medicinal Chemistry, 2018, 61, 5122-5137. 63. P. Johansson, K. Kaspersson, I. K. Gurrell, E. Bäck, S. Eketjäll, C. W. Scott, G. Cebers, P. Thorne, M. J. McKenzie, H. Beaton, P. Davey, K. Kolmodin, J. Holenz, M. E. Duggan, S. Budd Haeberlein and R. W. Bürli, Journal of Medicinal Chemistry, 2018, 61, 3491-3502. 64. J. D. Low, M. D. Bartberger, Y. Cheng, D. Whittington, Q. Xue, S. Wood, J. R. Allen and A. E. Minatti, Bioorganic & Medicinal Chemistry Letters, 2018, 28, 1111-1115. 65. Y. Tanaka, K. Aikawa, G. Nishida, M. Homma, S. Sogabe, S. Igaki, Y. Hayano, T. Sameshima, I. Miyahisa, T. Kawamoto, M. Tawada, Y. Imai, M. Inazuka, N. Cho, Y. Imaeda and T. Ishikawa, Journal of Medicinal Chemistry, 2013, 56, 9635-9645. 66. M. L. Stewart, E. Fire, A. E. Keating and L. D. Walensky, Nature Chemical Biology, 2010, 6, 595. 67. P. J. Hajduk, S. Boyd, D. Nettesheim, V. Nienaber, J. Severin, R. Smith, D. Davidson, T. Rockway and S. W. Fesik, Journal of Medicinal Chemistry, 2000, 43, 3862-3866. 68. C. W. Murray, O. Callaghan, G. Chessari, A. Cleasby, M. Congreve, M. Frederickson, M. J. Hartshorn, R. McMenamin, S. Patel and N. Wallis, Journal of Medicinal Chemistry, 2007, 50, 1116-1123. 69. K. Padmanabhan, T. P. Wu, K. G. Ravichandran and A. Tulinsky, Protein Sci, 1994, 3, 898-910. Page 34 of 35Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B 88x35mm (300 x 300 DPI) Page 35 of 35 Physical Chemistry Chemical Physics Physical Chemistry Chemical Physics Accepted Manuscript Published on 28 December 2020. Downloaded by Universitat Politecnica de Catalunya on 12/28/2020 10:40:28 AM. View Article Online DOI: 10.1039/D0CP05471B