Middle-way flexible docking: Pose prediction using mixed-resolution Monte Carlo in estrogen receptor α

There is a vast gulf between the two primary strategies for simulating protein-ligand interactions. Docking methods significantly limit or eliminate protein flexibility to gain great speed at the price of uncontrolled inaccuracy, whereas fully flexible atomistic molecular dynamics simulations are expensive and often suffer from limited sampling. We have developed a flexible docking approach geared especially for highly flexible or poorly resolved targets based on mixed-resolution Monte Carlo (MRMC), which is intended to offer a balance among speed, protein flexibility, and sampling power. The binding region of the protein is treated with a standard atomistic force field, while the remainder of the protein is modeled at the residue level with a Gō model that permits protein flexibility while saving computational cost. Implicit solvation is used. Here we assess three facets of the MRMC approach with implications for other docking studies: (i) the role of receptor flexibility in cross-docking pose prediction; (ii) the use of non-equilibrium candidate Monte Carlo (NCMC) and (iii) the use of pose-clustering in scoring. We examine 61 co-crystallized ligands of estrogen receptor α, an important cancer target known for its flexibility. We also compare the performance of the MRMC approach with Autodock smina. Adding protein flexibility, not surprisingly, leads to significantly lower total energies and stronger interactions between protein and ligand, but notably we document the important role of backbone flexibility in the improvement. The improved backbone flexibility also leads to improved performance relative to smina. Somewhat unexpectedly, our implementation of NCMC leads to only modestly improved sampling of ligand poses. Overall, the addition of protein flexibility improves the performance of docking, as measured by energy-ranked poses, but we do not find significant improvements based on cluster information or the use of NCMC. We discuss possible improvements for the model including alternative coarse-grained force fields, improvements to the treatment of solvation, and adding additional types of NCMC moves.


Introduction
Computational structure-based drug design can play an important role in drug 2 development, as exemplified in the development of inhibitors of HIV protease, which 3 have a major impact on treatment for AIDS patients. [1,2] Because of their potential 4 to reduce the cost and time associated with drug development, a multitude of methods 5 have been developed to screen potential drug candidates virtually and prioritize 6 possible structures for synthesis. [3][4][5][6] Some of the most popular docking methods 7 represent the potential energy due to the receptor using grids and optimize the ligand 8 conformation with respect to this potential. These include Autodock Vina, [7] and the 9 related smina [8], Schrodinger Glide, [9][10][11] CDOCKER, [12] and others. [13] Once a 10 grid is constructed, however, the protein conformation represented by that grid is 11 fixed, which is a serious approximation, since a number of structural studies have 12 shown that "hidden" protein conformations and protein flexibility play important roles 13 in protein-ligand binding; methods that consider protein flexibility and multiple 14 structures improve docking performance compared to those that make use of only one 15 structure. [14][15][16][17][18][19][20][21] 16 A number of approaches have been developed to incorporate flexibility into 17 docking. [22] One approach is to try to incorporate protein flexibility into grid-based 18 approaches by collecting multiple confomations of the receptor and docking to all of 19 them, a practice commonly known as ensemble docking. [14,15,18,21,23,24] Another 20 grid-based strategy involves leaving certain amino acid side chains out of the grid and 21 optimizing their conformation alongside that of the ligand during the docking 22 procedure. CDOCKER has been modified in this way for example, [25] and Autodock 23 Vina and smina are also capable of this. [7,8] Other examples of grid-based force 24 fields, some of which are aimed at the related protein-protein docking problem, include 25 the ATTRACT force field [26,27] and LightDock. [13] However, the amount of 26 flexibility that can be incorporated by both these methods is limited. In particular, 27 allowing only a few side chains to be flexible means these few side chains must be 28 carefully chosen and that important protein motions involving the backbone are not 29 represented. Likewise, ensemble docking requires a careful choice of conformations to 30 be used and only allows for a limited degree of backbone flexibility. Since these 31 conformations come from simulations or structures of the protein without the 32 corresponding ligand, ensemble docking cannot take into account the mutual induced 33 fit that may occur when a ligand binds to a protein. 34 The RosettaLigand approach to ligand docking [28][29][30] makes use of the Rosetta 35 knowledge-based force field and, in principle, allows for full receptor flexibility. Like 36 other knowledge-based force fields, it relies on the assumption that the system under 37 discussion is similar to known protein-ligand complexes. Furthermore, the 38 RosettaLigand docking approach used a complex protocol involving several rounds of 39 minimization, which loses information about the relative entropy of the energy minima 40 that are found. Although in theory the protocol should be able to allow full receptor Multiscale simulation techniques show promise for reducing computational time 59 while maintaining full flexibility of the protein and physical accuracy. They seem 60 especially appropriate for the docking problem because the regions of the protein far 61 from the ligand can be simulated with a less physically accurate, but less 62 computationally expensive, coarse-grained model. Several efforts have been made to 63 implement this idea in simulations, although not always in the context of 64 protein-ligand simulations. An early approach that bears some similarities our own 65 involves dividing the system into three regions: an atomistic region, a coarse-grained 66 region using a Gō-like potential, and an intermediate region. [38] Our group has also 67 done some preliminary work on mixed resolution models, combining a residuel level 68 Gō model with the OPLS-AA force field and conducting some preliminary tests on 69 self-docking to the estrogen receptor. [39] The popular MARTINI force field for water 70 has also been combined with an atomistic force field for proteins; [40] however, 71 balancing the strength of electrostatic forces in the two force fields proved difficult. 72 Others have tried to combine an atomistic region with an elastic network model. [41]. 73 Feig and co-workers have also combined their PRIMO force field with the 74 CHARMM36 atomistic force field, [42] and obtained similar results to fully atomistic 75 or fully coarse-grained simulations, but they note issues with weakened hydrophobic 76 packing interactions, and the amount of speedup they obtained is generally modest 77 because PRIMO is relatively detailed for a coarse-grained model. A few workers have 78 also attempted to adapt grid-based docking methods to work in a multiscale manner; 79 some of these approaches are aimed at the related protein-protein docking 80 problem. [13,26,27]

81
Here we extend our previous mixed-resolution Monte Carlo software [39] and use it 82 to systematically study key aspects of docking in a challenging system. In the 83 software, the majority of the protein is modeled using a residue-level coarse-grained 84 model currently based on the Gō model [43][44][45][46] while the ligand and binding site are 85 modeled using an atomistic force field. Interactions between the two regions are 86 treated in a fully atomistic manner. Full flexibility of the ligand and receptor is 87 maintained, allowing the modeling of mutual induced fit, including "breathing" 88 motions away from the binding site. [47,48] Importantly, the coordinates of all atoms 89 are tracked throughout the simulation, which increases computational cost but enables 90 maintaining proper backbone geometry via standard non-bonded terms even in the 91 coarse region. The computational cost of the method can be adjusted by changing the 92 size of the atomistic region; for the region used here, the mixed-resolution model 93 reduces the computational cost by a factor of 2-4 compared to a fully atomistic 94 treatment of the protein using the same implementation. In addition, coarse-grained 95 models produce potential energy surfaces that are smoother than corresponding 96 atomistic surfaces, because of averaging over omitted degrees of freedom. [49] While 97 we use a Gō-like model for the coarse-grained region in the protein in this work, the 98 mixed-resolution concept can be extended to incorporate other coarse-grained models.

99
The use of Monte Carlo simulation with full flexibility in principle also allows for the 100 inclusion of entropy effects.

101
In this work, we have replaced the discontinous Gō model functional form used in 102 previous work [45,46] with a continous Lennard-Jones functional form. We have also 103 replaced the OPLS-AA force field used in our previous work with the AMBER 99SB 104 force field [50] so that ligands can be automatically parameterized using the 105 compatible GAFF forcefield and antechamber tool. [51,52]  which the co-crystal structure for the tested compound is not used.

145
The first section of the paper describes the mixed resolution potential and the 146 Monte Carlo and NCMC protocols that were used with it for docking. For this study, 147 ligands with known bound crystal structures were used, so we can compare the docked 148 conformations to experimental crystal structures. We study the docking protocol with 149 and without NCMC and with different levels of protein flexibility, and also tested  constitute the atomistic region and are treated using an all-atom force field. The remainder of the protein is treated as the coarse-grained region is represented by particles located at each α carbon (purple spheres) with native attractions between them (yellow lines). Rigid structures of each amino acid (thin structure) and moved along with the coarse-grained particles and used to calculate interactions between the coarse-grained and atomistic regions. Helix 12 is indicated in orange in all three panels.

163
Simulations were done with a mixed-resolution potential, in which the majority of the 164 protein is treated with a Gō model [43][44][45][46]65] and the atomistic region around the 165 binding site is treated using the AMBER 99SB force field. [50] The overall potential 166 has the following form: The coarse-grained portion of the potential U CG is in turn given by where the sum is taken over all pairs of residues in which both residues belong to the 169 coarse-grained region. The coarse-grained potential also includes a backbone terms are not included). The 12-10 form used here replaces a square well potential 173 used in previous work [45,46] and has shown good performance in other 174 settings. [65][66][67] The Gō model well depth ε and hard core radius r HC are listed in 175 table 1. r 0 ij is the native distance between atoms i and j, determined from the original 176 structure.

177
The all-atom potential U AA is the potential energy of the all-atom region, 178 according to the AMBER 99SB force field. [50] The U CG/AA term represents the 179 interaction between the two regions, which is also computed atomistically using the 180 AMBER 99SB force field. The all-atom region was defined to contain the ligand and 181 all residues with at least one non-hydrogen atom within 3Å of the ligand in any of the 182 66 crystal structures of ligand-ER complexes that were used as references; in addition, 183 residues 533-548, which comprise helix 12 and the neighboring loop, were also included 184 in the all-atom region. Ligands were parameterized with the antechamber tool [52] 185 using the GAFF force field. [51] 186 In order to incorporate solvation effects in a computationally efficient manner, the 187 electrostatic term of the AMBER 99SB was modified to use a solvent-exposure 188 dependent distance dependent dielectric originally developed by Garden and Zhorov,189 which was previously found to work well in docking simulations with an AMBER force 190 field. [53] In this model, the electrostatic interaction between atoms i and j is given by 191 where ε 0 and ε 1 are low and high dielectric constants, respectively, and s kl is the    In order to enhance the sampling of ligand poses, docking runs were also undertaken 209 using a modified version of the nonequilibrium candidate Monte Carlo (NCMC) 210 algorithm. [54,55] In this method, the system is subjected to NCMC "moves" which 211 are in fact short nonequilibrium simulations (using standard MC) during which the 212 potential energy function is occasionally perturbed such that sampling for particular 213 degrees of freedom may be enhanced.

214
In the method used here, the potential energy is modified by scaling the van der 215 Waals and electrostatic components of the protein-ligand interaction energy by the 216 factors λ VDW and λ elec : Each NCMC move consisted of 800 individual Monte Carlo moves, which were divided 218 into four 100-move phases in which the coupling parameters λ VDW and λ elec were 219 changed according to the schedule shown in Fig. 2. First, the charges on the ligand 220 were removed by systematically driving λ elec towards 0. Second, λ VDW was also driven 221 towards 0, such that the ligand was completely uncoupled from the protein and free to 222 rotate, translate, or change conformation without interference. Then, in the second 223 two phases, first λ VDW and then λ elec weree gradually transitioned back to 1, so that where U is the scaled potential given above. This leads to canonical sampling suitable 229 for the given set of λ values.

230
At the end of each NCMC move, the nonequilibrium work w performed on the 231 system (which accounts for the changes in λ values) was calculated using If the NCMC move was rejected, the conformation of the system prior to the sequence 238 of MC steps making up the given NCMC move was restored. Because of this, and the 239 low acceptance rate of NCMC moves (about 2%), in pure NCMC simulations the 240 relaxation of the system toward low energy conformations was extremely slow.

241
Consequently, in order to promote more rapid relaxation, the NCMC moves were corresponding reference structures are found in tables S1 Table and S2 Table.   hierarchical clustering [74] was then applied to this distance matrix to construct a 299 "phylogenetic tree" using ligand heavy atom RMSD taking into account chemically 300 equivalent groups. The clusters were then defined from this tree, using a cutoff defined 301 as the 10th percentile of all distances between structures for each drug. This ensured 302 that the size of the clusters in configuration space would be appropriately scaled to the 303 overall distribution of the poses for each drug.

306
We first examined the impact of flexibility on pose generation. Fig. 5 shows  drugs. This average RMSD is then plotted as a function of N . It is clear that a much 338 greater proportion of these best poses are close to the corresponding crystal structures 339 when the protein is allowed to be fully flexible than when the protein is rigid or only 340 sidechains move, and that the average best RMSD achieved is lower as well when more 341 flexibility is allowed. This is especially true for the antagonists, which are generally 342 more flexible than the agonists and whose docking is consequently more challenging.

343
This difference is probably because the inactive structure of ER α is more open and 344 consequently the protein as a whole, particularly helix 12, is more flexible.

345
Consequently, the additional flexibility provided by the MRMC method may be more 346 important for the inactive state.

347
In principle, allowing protein flexibility should enable us also to predict the changes 348 in protein structure that occur upon docking to different ligands. We constructed 349 Ramachandran and Janin plots [75] for the amino acid residues in the atomistic region 350 for the final conformations from our docking runs. Fig. 7 shows these plots for two   probability (given by eq. 9) as a function of the overall displacement or rotation angle, 377 were both computed and plotted as shown in Fig. 8. The corresponding distribution 378 of move sizes and average acceptance probability were also plotted for comparison.

379
The plots show that, compared to regular MC, the distribution of moves generated 380 by NCMC favors smaller moves in both translation and rotation. In addition, while 381 acceptance rates for small moves are comparable for NCMC and MC, the acceptance 382 rates for larger moves are many orders of magnitude smaller for NCMC. This is 383 particularly (and unexpectedly) true for ligand rotations, where acceptance rates for 384 all but the smallest rotations are much smaller for NCMC than for MC. The combined 385 effect of these two trends is that in NCMC, a greater number of small translations and 386 rotations are generated and accepted, compared to MC. The reason for these results requires large nonequilibrium work and consequently leads to rejections. We note that 390 these results are specific to our implementation, as described below in the Discussion.      all three of these methods. The differences between them are small, but it appears 421 that method 3 above performs slightly better than the others in that the 422 conformations identified by this method are closer to the crystal structures.

423
Note that we did not employ a more quantitative entropy estimation process 424 because the sampling did not appear to be sufficient -i.e., there were very few jumps 425 between poses in a given MC run ( Fig. 9) implying true Boltzmann sampling was not 426 achieved. Also, using the total energy in place of the protein-ligand interaction energy 427 and combining it with similar clustering approaches gave similar results.

428
Multiple Poses for WAY-169916 429 The crystal structure of the ER partial agonist WAY-169916 (PDB code 3OS9, ligand 430 ID KN1) [17] bound to the inactive conformation of ER α shows two separate poses 431 for the ligand, with occupancies of approximately 70% and 30%. This provided an 432 opportunity to test whether our protocol can find multiple bound poses for a ligand.

433
The RMSD to each of these poses was calculated separately for our docking 434 simulations of this drug, and a scatterplot is shown in Fig. 11. The ensemble of   Table 3 shows a comparison of computation speeds for fully atomistic,

463
In this paper, three separate tactics for improving protein-ligand docking were tested. was also evaluated by studying the relationship between acceptance rate and move size 472 for NCMC moves. We found that allowing for full protein flexibility using the 473 mixed-resolution potential significantly improved the docking results, and the use of 474 NCMC produces a further modest improvement. However, clustering did not appear 475 to offer any significant advantages over simply ranking the poses by protein-ligand 476 interaction energy. Overall, we obtained a correct pose within 2Å for about half the 477 ligands, so there is significant room for improvement. Below, we discuss ways that 478 several aspects of the approach could be improved. We note that systematic 479 investigation of these different aspects of the docking problem is facilitated by having  However, with other target proteins this may not be adequate. Additional flexibility 500 could be incorporated by making use of a double well Gō potential [45] or by replacing 501 the Gō potential by another potential that is less dependent on native interactions.

502
The popular MARTINI force field [76,77] has been used in a mixed resolution 503 configuration [40] but frequently leads to distorted structures for soluble proteins 504 unless reinforced by an elastic network model. [78] Other potential force fields that 505 could in principle be used include OPEP [79] and UNRES. [80] We have also 506 developed a tunable coarse-grained force field based on constructing interaction energy 507 tables and applying variable amounts of smoothing to them [81]; one of the goals for 508 this force field was to use it in a mixed resolution setting, but substantial additional 509 effort would be required to implement that combination.

510
Compared to the coarse-grained potential, the atomistic potential used here would 511 seem to offer less room for improvement. The AMBER 99SB force field is a 512 well-tested, commonly used force field, although it could conceivably be replaced with 513 another, newer force field. A more significant area for improvment concerns the 514 treatment of solvation in the simulations. The SEDDD method is based on a distance 515 dependent dielectric with a dielectric constant that includes some solvent exposure.

516
The linear dependence of the dielectric constant on interatomic distance is not 517 physically correct, however, since the dielectric constant should approach that of water 518 as the distance between two atoms increases. There is also no explicit term 519 representing the hydrophobic effect. A generalized Born model [82] would be more 520 physically realistic but also more computationally expensive. Explicit solvent (using 521 water molecules restrained around the atomistic region) would be even better, but 522 would require even more computational cost, and steric clashes between the ligand 523 and water molecules would make it more difficult to sample ligand configurations via 524 Monte Carlo. Another relatively inexpensive solvation model is the Sheffield solvation 525 model, [83] which might provide improved results.

526
The choice of atomistic region also plays an important role in establishing the 527 tradeoff between physical accuracy and computation time. We selected the atomistic 528 region used here to be as small as possible while including all those residues in direct 529 contact with the ligand in any of the reference structures. We also included helix 12 dihedrals. [84] We may be observing the same effect here, where couplings between the 557 ligand position and orientation and the side chain degrees of freedom of nearby amino 558 acids reduce the efficiency of NCMC. Alternatively, as suggested by Chodera [personal 559 communication] our NCMC protocol may require better targeted ligand MC moves 560 such as "smart darting," [85] although our initial tests of this idea did not show a 561 substantial improvement.

562
There are also straightforward means to improve sampling. Two simple ways are to 563 run longer docking simulations or perform more docking runs for each drug. The time needed to dock a drug. In addition, using a longer or otherwise redesigned 568 schedule for λ VDW might improve the NCMC acceptance rate, which was relatively 569 low (approximately 2%) in the work reported here, and thereby improve the sampling. 570 Another way to improve the NCMC acceptance rate might be to use a soft core 571 potential for the van der Waals interactions. This might reduce the potential energy 572 changes associated with steric clashes between the ligand and receptor as λ VDW is 573 increased from 0 to 1, which contribute to large values of the nonequilibrium work and 574 consequently to poor acceptance rates. conformations should also give information on the relative free energies of different 582 binding conformations, which can then be used to calculate binding free 583 energies. [86][87][88] Motivated by this reasoning, we sought to apply clustering algorithms 584 to the ensemble of configurations we obtained from our docking simulations and use 585 this information to aid in identifying the most representative structures. We found, 586 however, that clustering information did not improve the ranking of poses in practice.

587
The main reason why the clustering was not useful may simply have been that the 588 protocol used here did not allow for adequate sampling, as described above. Another 589 flaw may have been the choice of clustering algorithm or metric used. The complete 590 linkage clustering algorithm used here is relatively crude, being primarily intended for 591 the construction of phylogenetic trees using distances between protein or DNA 592 sequences. [74] Despite this, it was selected because many other algorithms, such as 593 the K-means algorithm, rely on averaging coordinates from distinct configurations, an 594 operation of unclear physical meaning. In addition, complete linkage clustering 595 guarantees that any two configurations classified in the same cluster will have an 596 RMSD less than the selected cutoff. The Cheatham group has tested a number of 597 clustering algorithms on MD trajectories; [89] while they recommend average-linkage 598 hierarchical clustering for circumstances in which the number of clusters is not known 599 in advance (as here) they also point out that the performance of a clustering algorithm 600 is influenced by the choice of atoms used for pairwise comparison and that hierarchical 601 clustering is sensitive to outliers. pipelines. This study is a step toward developing such a highly adaptable platform. 616 We recognize that further improvements to sampling and entropy-based pose 617 evaluation will be necessary to make a middle-way tool more valuable for the 618 drug-design enterprise.