Ligand design by targeting a binding site water

Solvent reorganization is a major driving force of protein–ligand association, but the contribution of binding site waters to ligand affinity is poorly understood. We investigated how altered interactions with a water network can influence ligand binding to a receptor. A series of ligands of the A2A adenosine receptor, which either interacted with or displaced an ordered binding site water, were studied experimentally and by molecular dynamics simulations. An analog of the endogenous ligand that was unable to hydrogen bond to the ordered water lost affinity and this activity cliff was captured by molecular dynamics simulations. Two compounds designed to displace the ordered water from the binding site were then synthesized and evaluated experimentally, leading to the discovery of an A2A agonist with nanomolar activity. Calculation of the thermodynamic profiles resulting from introducing substituents that interacted with or displaced the ordered water showed that the gain of binding affinity was enthalpy driven. Detailed analysis of the energetics and binding site hydration networks revealed that the enthalpy change was governed by contributions that are commonly neglected in structure-based drug optimization. In particular, simulations suggested that displacement of water from a binding site to the bulk solvent can lead to large energy contributions. Our findings provide insights into the molecular driving forces of protein–ligand binding and strategies for rational drug design.

DDH and DDS values were determined using MD/FEP from the slope and intercept of a linear regression of DDG/T versus 1/T (Fig. 3a). The DDG values are averages calculated based on three independent sets of MD/FEP calculations at 300 K. Uncertainties for DDH and -TDDS at 300 K were calculated as the SEM from eight values obtained from leave-one-out analysis of the van't Hoff plot. Potential energy values are averages ± SEM from 10 independent simulations (1 ns each). b Relative ligand potential energy representing ligand-receptor, ligand-water, and ligand-ligand interaction energies. c The potential energy from within the surroundings (DDUS) was calculated from the difference between the total enthalpy (DDH) and the ligand enthalpy (DDUL).  Each receptor-ligand complex was then prepared based on a snapshot of the equilibrated system. The binding modes of compounds 2-4 were modelled based on the A2AARadenosine complex and waters overlapping with the ligands were removed. The MD/FEP calculations were performed with Q 9 using the OPLS force field implemented for this program 10 . The simulations were performed under spherical boundary conditions with a sphere radius of 21 Å centered on the ligand (unless noted otherwise). The OPLSAA_2005 force field from the ffld_server program (Schrödinger, LLC, New York, NY, 2017) combined with 1.14*CM1A-LBCC partial charges were used to parameterize compounds 1-4 11 . Atoms outside the sphere were excluded from non-bonded interactions. Ionizable residues close to the sphere edge were set to their neutral form and atoms within 3 Å of the sphere edge were restrained to their initial coordinates. Water molecules at the sphere border were subject to radial and polarization restraints according to the surfaceconstrained all atom solvent (SCAAS) model 12 . Solvent bonds and angles were constrained S10 using the SHAKE algorithm 13 . A cutoff of 10 Å was used for non-bonded interactions except for ligand atoms, for which no cutoff was applied. Long-range electrostatic interactions were treated using the local reaction field approximation 14  heavy atoms were gradually released and the system was heated to the target temperature.
This was followed by a 250 ps production phase and energies were collected every 50 fs.
Each ligand was also prepared for simulations in aqueous solution using a water droplet of the same size. In this case, the equilibration and production phase times were 350 ps and 100 ps, respectively. In these simulations, a harmonic restraint was applied to a central ligand atom to prevent it from approaching the sphere edge. Free energies were calculated using the Zwanzig equation and a thermodynamic cycle, as described previously 16 .

S11
Hydration site analysis were carried out using GROMACS 5 and SSTMap 17 using the same force field as for the MD/FEP calculations. Three replicas of each complex were simulated with solute heavy atoms tightly restrained to their initial coordinates. The systems were first heated to 300 K in the NVT ensemble, which was followed by a 1 ns equilibration step in the NPT ensemble. Subsequently, three independent NVT simulations of 10 ns of each complex were performed and 10,000 snapshots were collected. Clustering in SSTMap was used to identify ordered waters, followed by calculation of enthalpies and entropies for these hydration sites using the approach by Young et al. [17][18][19] .
Molecular Docking. Docking calculations with AutoDock Vina were performed using default settings in a 15 Å box centered on the ligand 20 . The GLIDE docking calculations were carried out using the OPLS3e force field and GLIDE-SP scoring function 21 . In the GLIDE calculations, the ring puckering of adenosine was fixed in order to obtain the same binding mode as that observed the adenosine-bound A2AAR crystal structure 1 .
Biological Assays. Radioligand binding assays were performed as previously described 22 using membrane preparations from HEK293 cells stably expressing the human A1AR, Data are expressed as mean ± standard error.
Functional assays were carried out in HEK293 cells stably expressing a single hAR subtype. The cells were cultured in DMEM supplemented with 10% fetal bovine serum,  Alternative synthesis route for compound 3. As described in Scheme S2, the amino group of 7-methyl-1H-imidazo[4,5-c]pyridin-4-amine 5 was protected with benzoyl chloride in the presence of TMS-Cl to afford the corresponding protected derivative 12 in 43% yield (slightly higher than with M4S). Compound 12 was subjected to a Vorbrüggen coupling with 8 using TMSOTf to give the corresponding undesired β-N 5 -ribosylated derivative/riboside (13, kinetic product). Initially, we expected that the coupling product would be the more stable β-N 1 -ribosylated derivative (16, thermodynamic product), which was difficult to differentiate from 13 by its mass and 1 H-NMR spectrum in chloroform-d.
Consequently, the riboside derivative 13 was treated with sat. NH3 in MeOH for 20 h, and the product 15 was isolated as a white solid. However, subsequent analysis revealed that the product was neither the N 1 -nor the N 3 -ribosylated derivative, but the β-N 5 -ribosylated

S15
Framski et al. 24 . The assigned structure was confirmed by NMR. The absence of the N 1addition product under these conditions was probably due to the shielding of the N 1position by the pyridine 7-methyl group.
When the solvent was changed from acetonitrile to toluene, the more stable thermodynamic Several attempts were made to brominate the methyl group in compound 16 using the conditions as used for compound 9, but no desired mono-bromo compound could be isolated, although a product peak was observed by mass spectrometry. We did not optimize the reaction conditions for the bromination of 16. Had we been able to isolate the desired bromo derivative, compound 4 could have been prepared by this route. Based on these results, it is emphasized that the Vorbrüggen reaction regioselectivity with protected forms of 5 depends upon the solvent, temperature and the protecting group/substituent on N 4atom. The overall yield of compound 3 using N-benzoyl protection was 7.6%, and we feel that overall yield could be improved by optimizing reaction conditions. However, the more general synthetic approach proved to be using the M4S protecting group.