Probing Hydration Patterns in Class-A GPCRs via ... - ACS Publications

Article pubs.acs.org/JCTC

Probing Hydration Patterns in Class‑A GPCRs via Biased MD: The A2A Receptor Syeda Rehana Zia,† Roberto Gaspari,‡ Sergio Decherchi,†,§ and Walter Rocchia*,† †

CONCEPT Lab and ‡CompuNet, Italian Institute of Technology, Via Morego, 30, I-16163 Genova, Italy BiKi Technologies s.r.l., Via XX Settembre, 33/10, I-16121 Genova, Italy

§

S Supporting Information *

ABSTRACT: Herein, we present a new computational approach for analyzing hydration patterns in biomolecular systems. This protocol aims to efficiently identify regions where structural waters may be located and, in the case of protein−ligand binding, where displacing one or more water molecules could be advantageous in terms of affinity and/or residence time. We validated our approach on the adenosine A2A receptor, a target of significant pharmaceutical relevance. The results of the approach are enriched with an extensive analysis of hydration in A2A and other members of the A-family of GPCRs using available crystallographic evidence and reviewing existing literature. As per the protein−ligand complex case, we conducted a more detailed study of a series of triazine analogues inhibiting A2A. The proposed approach provides results in good agreement with existing data and offers interpretability and simple and fast applicability.

1. INTRODUCTION The interest in understanding how water exerts its essential role for life, particularly in supporting protein function, continues to grow. Increasingly often researchers are finding water molecules that are involved in maintaining protein structure and dynamics, in forming protein−ligand complexes, in allosteric regulation of receptors, in protein−protein interactions, and in other biological mechanisms, many of which are crucial for drug discovery.1−18 However, among the molecules comprising the hydration sheath, only a small subset can be termed “structural” in that they are crucial for a protein’s structure (and thus its activity), for instance, by forming a bridge between different protein regions that are distant in terms of sequence. Evidently, water molecules and hydration networks can both affect and be affected by the presence of a ligand in the binding site.3−6,19 Gaining insights into hydration can help drug designers engineer a compound that, upon binding, displaces waters in duly identified regions or, alternatively, establishes bonds with them. This ability can directly impact the formation of the complex in terms of thermodynamics, kinetics, or both. It has in fact been suggested that a ligand displacing “unhappy” water molecules upon binding may significantly affect its association kinetic rate.2,20 More ambitiously, displacing (or conversely stabilizing) structural waters can be instrumental in affecting the conformational equilibrium of the target protein’s structure. Finally, some studies have found an improvement in the performances of docking/virtual screening runs if a few carefully chosen water molecules were considered as being part of the binding site.21,22 © XXXX American Chemical Society

Along this line, a number of computational methods have focused on incorporating active site water molecules into structure-based drug design to confirm and enrich experimental approaches.16,23,24 They include JAWS,25,26 GCMC,27 WaterFLAP,28 WaterMap,12,16,29 SZMAP,23 STOW,30 a grid-based implementation of solvation theory (GIST),31 and many others.32−35,27,36 Most of these methods try to estimate a regional hydration free energy to quantify the “happiness” of water molecules inside or around a binding site. However, this estimation process is not straightforward; in particular, where hydration is low and where water molecules may become trapped, convergence may be hard to reach. As a consequence, taken individually, different computational methods may disagree in some cases. Nevertheless, a machine-learningbased approach, combining some of the techniques mentioned above, has managed to successfully correlate water displacement with binding affinity on more than 350 congeneric ligands.2,37 Class A G-protein-coupled receptors (GPCRs) are a relevant example of a protein family for which a significant number of conserved water molecules have been observed and discussed. By analyzing crystallographic structures across the family, researchers have identified these waters in a number of regions, organizing them in at least 15 different structural hydration sites (or water clusters).38 It appears that the main role of these sites is to create hydrogen bond networks that strengthen Received: May 10, 2016

A

DOI: 10.1021/acs.jctc.6b00475 J. Chem. Theory Comput. XXXX, XXX, XXX−XXX

Article

Journal of Chemical Theory and Computation

to label the water clusters in the GPCRs’ structures.42,38 Biased simulations have been performed on the same system in the apo form and complexed with triazine derivative inhibitors, which we named 4a, 4e, 4f, and 4g in accordance with Bortolato et al.’s nomenclature.2,37 Lengths in Å between braces represent distances. Residues are identified by the three-letter code followed by residue number when referring to the simulated systems, whereas Ballesteros numbering is used when referring to the analysis of conserved residues within the GPCR family A.43 Crystallographic waters are indicated as WXXXX, with XXXX being the residue number in the high-resolution A2A crystal structure with the PDB code 4EIY, unless otherwise stated. Throughout the text, occupancy maxima are named as “ΓxY”, where x identifies the specific run (i.e., “cli” for the i-th cluster analyzed by plain MD simulation of the apo structure, “bias” for the apo-biased runs, and “4a”, “4e”, “4f”, and “4g” for the holo-biased simulations). Each subscript is followed by the identification number Y of the maximum. For the sake of simplicity and when ambiguity is avoided, the superscript has been dropped in the figures and tables. 2.2. Reference for Water Conservation Analysis. To contextualize discussion of the obtained results, we rely on previously published structural analysis, where researchers have classified conserved residues and conserved water molecules.38 In particular, multiple-sequence-alignment-based comparison of five members of the A GPCR family has revealed 24 conserved amino acid residues located in the transmembrane (TM) helices (see also the Supporting Information for a multiple sequence alignment including the 4EIY sequence). Some of these are located in the highly conserved (D/E)RY motif (3TM), the WxPF/Y motif (6TM), and the NPxxY motif (7TM). Structural superimposition and analysis of these GPCRs led to the identification of 15 regions where water molecules are colocalized (named “water clusters”). We will refer to these regions using the naming convention proposed in Angel et al.38 These conserved waters seem to exert their structural role by making contacts with residues in the vicinity of the transmembrane helices.38 Conserved waters and structurally important residues contribute to the formation of so-called major hydrogen-bond networks (MHNs). Within the A-family GPCRs, MHNs connect several TM residues, including Asn1.50, Asp2.50, Ser3.39, Tyr5.58, Trp6.48 (in the conserved CWxP motif), Asn7.45, Ser7.46, Asn7.49, and Tyr7.53 (in the conserved NPxxY motif).44 Within the A family, several works have focused on the adenosine A2A receptor. One recent A2A study was based around a highresolution crystallographic structure (PDB ID: 4EIY).42 In contrast to Angel et al.,38 Liu et al.42 classified experimentally detected water molecules into extracellular (EC), central (CC), and intracellular (IC) clusters. To include these results in our analysis, we sought a correspondence between the different classifications. This is reported in the Supporting Information together with a more detailed description of water clusters and their postulated structural role. 2.3. System Setup. We considered the apo and four holo systems, namely, the 4a, 4e, 4f, and 4g complexes of A2A-staR, as named in Bortolato et al.2 The staR structure has the advantage of being stable in water. This allows for faster simulations and avoids the need for a large membrane model, which is not relevant for our purposes because we are not focusing here on water exit pathways. Incidentally, a visual inspection analysis of our trajectories did not provide evidence that water molecules exit from confined regions through paths

interactions between the different helices of the transmembrane receptor, and recent simulations have also proposed that the concerted motion of water networks can be key in regulating the receptor activity.39,40 Furthermore, for some members of this protein family, hydration shows a significant impact on ligand affinity and binding kinetic rate constants. In this respect, several experimental and simulative studies have been performed on the A2A receptor. For instance, substitutions affecting hydration sites were correlated with the binding affinity of a series of congeneric A2A inhibitors by means of the WaterMap tool.29 On the basis of the resulting evidence, the continuing increase in the available structural and functional information, and the tremendous biological and pharmaceutical importance, GPCRs can undoubtedly be considered an excellent test bed for methods that seek to characterize and predict hydration and its effects. In the present work, we propose a reductionist approach to identifying the role of individual water molecules in a protein’s structure and propensity for ligand binding. This approach divides loosely bound water molecules, which are natively short-lived, from tightly bound ones. Accordingly, we devised a molecular dynamics (MD)-based protocol that seeks to quickly identify most persistent waters in any predetermined region of a biomolecular system, but especially in a protein’s binding site. In this protocol, we force the desolvation of a desired region via a steered MD protocol acting on a purposely defined collective variable (CV). This bias produces an external force field that mimics a repulsive electrostatic field between the selected region of the protein (often a pocket) and the nearby water molecules. This bias is realized by suitably positioning on the protein fictitious charges that (coupled with the regular van der Waals forces exerted by the target) consistently and dynamically repel water oxygen atoms. The method is thus suitable, for example, for simulating protein cavity dehydration because it counteracts, by construction, atomic movements that would bring water molecules back into the cavity or against the protein surface. In the proposed method, short runs of steered MD are followed by a spatial density analysis to measure the local water persistency against the application of the bias. Complementary to the persistency estimation part is the identification of water molecules that may behave as persistent only due to steric reasons rather than to favorable hydrogen bond interactions. Identification of these cases is performed by stripping explicit water molecules from each MD frame and then computing the solvent-excluded surface.41 The degree of bulk accessibility calculated by averaging the resulting map along the simulation is, in our opinion, a good figure for estimating “trapping regions” and distinguishing them from sites where high water persistency is due to stabilizing interactions. We validated our approach by applying it to the adenosine A2A receptor both in apo form and complexed with a series of triazine derivative inhibitors. The method provided, at low computational cost, qualitative but still insightful information on structural waters and an easier grasp of the role of hydration in A2A inhibitor binding, obtaining good agreement with previous studies and with available experimental data.2

2. METHODS 2.1. Nomenclature, Notation, and Residue Numbering. The studied system is the staR (stabilized Receptor) adenosine-receptor A2A. To deliver the most comprehensive analysis, we used both Liu et al. and Angel et al.’s nomenclature B


Article


where the c constants are suitably normalized charges, r is the distance between the interacting particles x and y, and k is a spatial decay constant. This biasing potential gives rise to a force between the oxygen atom of any water molecule and any protein heavy atom of the prescribed protein region in a way that mimics a screened electrostatic repulsion. The screening, described by the exponential decay significantly mitigates the long-range nature of the Coulomb (1/r) term, avoiding the need of the numerical treatment of long-ranged forces and focusing the effects of the bias on a neighborhood of the predetermined repulsive region. The pairwise additivity is wellsuited to the force fields currently used in MD simulations. By appropriately choosing the “charged” residue’s location, it is possible to obtain field lines that naturally, consistently, and dynamically point out a region as concave and tortuous as a binding site. Although there may be many diverse ways to devise such a repulsive bias, the behavior of electrostatic field lines, obtained suitably combining w(r)-like terms, seems to be ideal for accelerating desolvation. A comparison with the simplest “coordination” variable is provided as Supporting Information. This form is implemented in Plumed2 and can be applied to any MD engine supported by Plumed2 itself.65 The initial choice of the amplitude of the c constants is not relevant: the CV is calculated along a reference short simulation with predetermined default parameters, and then the protocol steers the CV until it reaches the desired fraction of the initial value. As per our preferred protocol, let Group A be the set of heavy atoms of the binding pocket, Na+ ion, and of the ligand (if present), and let Group B be the set of all water oxygen atoms. A neighbor list is used to improve the performance. As with other nonbonded interactions, this list is updated during runtime. The accelerated desolvation is obtained by reducing the value of the electrostatic-like interaction energy among the groups. This triggers a collective repulsive effect between Groups A and B. We devised two desolvation schedules: a slow one in which the CV is decreased along several plateaus of 2 ns, which was mainly used for validating and tuning the protocol, and a fast schedule (see Figure 1) that envisions only two plateaus. In the first, the CV is left at the initial reference value (w0). In the second, the value of the CV is reduced to half of w0 in a short time frame. The reference value w0 of the collective variable can be computed from the system’s initial configuration, which usually consists of the hydrated and equilibrated system. To do this, we used the Plumed2 driver on the last equilibration step.65 Overall, the fast protocol comprises three phases: (i) the reference value w0 is maintained for a prescribed time, (ii) its value is halved using a decreasing ramp, and (iii) this halved value is kept constant until the simulation ends. In all the following results, the k value was parametrized to obtain a decay distance of 10 Å, and the c constants per each atom of the same group were equal in value and sign. More complex scenarios could be envisioned, but this appeared to be sufficient to achieve the desired results. 2.5. Choice of the Repulsive Region and Persistence Estimation. In the A2A system, we computed the collective variable value by considering, in Group A, all the heavy atoms of the pocket as detected by the NanoShaper software tool66,67 together with the Na+ ion and, in Group B, all the oxygen atoms of water (within a neighborhood of 7 Å). The binding site residues were thus Asp52, Val55, Ala59, Ala63, Ile66, Ile80, Ala81, Cys82, Phe83, Val84, Leu85, Ala88, Gln89, Ser91, Ile92, Leu167, Phe168, Glu169, Met177, Asn181, Trp246, Leu249,

that would intersect the lipid bilayer in the membraneembedded receptor model. Crystal structures are available for the 4g (PDB ID: 3UZA) and 4e (PDB ID: 3UZC) complexes. Because no crystal structure of the apo A2A-staR is yet available, its structural model was taken from the 3UZC structure after ligand removal. The unresolved missing part of the ECL2 (extracellular loop 2, residues 150−157) was inserted in the apo and holo models using the Modeller software and the 4EIY structure as a template.45,46 To determine the protonation state for the titratable residues, we used the H++ server,47−49 which estimates the residue’s pKa using a Poisson−Boltzmann-based method.50 We calculated the pKa for different values of the protein-relative dielectric constant (ε = 2, 6, and 10) and at a bulk ionic strength of 0.15 M. The continuum solvent-relative dielectric constant was always set to 80. For His278, this method yielded a pKa value larger than 12 for all parameter combinations. We thus assumed the doubly protonated state. Similarly, we used the H++ server predictions to choose between ε and δ histidine protonations. All other titratable residues were kept in their standard protonation state at pH 7. The structural model used for MD simulations includes four disulfide bridges linking Cys166-Cys77, Cys146-Cys74, Cys159-Cys71, and Cys259-Cys262. The model included the structurally important sodium ion, which is present in 4EIY but undetected in the 3UZC/3UZA structures. The initial position of the ion was retrieved after structural alignment of the backbone carbons shared by the 4EIY and 3UZC structures. The initial positions of the 4a and 4f ligands were obtained by structural alignment on the 4e ligand scaffold. The protein force field used for the MD simulations was Amber99SB-ILDN with TIP3P waters. Ligands were parametrized using GAFF51 and RESP52 for partial charges (basis set 6-31g++ at the DFT-B3LYP level of theory) via the R.E.D. server.53 The ligand-protein system was solvated in a triclinic (dodecahedral) unit cell and neutralized by adding 10 Cl− ions. We used the Gromacs 4.6.1 simulation software for all simulations.54 In the starting configuration, the minimum distance between the protein and the box edges was set to 12 Å. All MD equilibration steps were performed by constant annealing from 0 to 297 K in 200 ps and by further 20 ns simulation at 297 K in the NPT ensemble using the velocity rescaling thermostat55 and Parrinello−Rahman barostat56−58 with coupling time constants of 0.1 and 2 ps, respectively. Target pressure was set to 1 atm. During equilibration, the protein, ligand, and Na+ ion were restrained by a weak positional harmonic restraint with a force constant of 10,000 kJ mol−1 nm−2 for all of the systems. The intensity of this restraint is approximately one-fourth of that used in other grid-based approaches, such as in GIST,31 but nevertheless, it limits the fluctuations to acceptable values (see section 3.3.5). The simulation time step was set to 2 fs. 2.4. Dehydration Protocol. The proposed method aims to quickly estimate water persistence in a given region by using a dehydrating bias. The bias acts on a collective variable that mimics the interaction energy of charged particles in an electrolytic solution modeled according to the Debye−Huckel theory.59−64 It can be written as the sum of several terms of the form w(r ) = cxcy

e−kr r C


Article


His250, Asn253, Thr256, His264, Leu267, Met270, Tyr271, Leu272, Ala273, Ile274, Ser277, His278, and Asn280. When the ligand was present, the corresponding heavy atoms were also included in the first group. During the biased simulations, the protein, ligand, and Na+ ion were subjected to a harmonic positional restraint with a spring constant of 10,000 kJ mol−1 nm−2. This allowed us to consistently define a grid along the trajectory and to avoid unwanted distortions of the protein induced by the bias. This also restricts the validity of the results to the starting pose, which is however consistent with our aims. We then obtained information on water molecule persistence under bias application by calculating water density on a grid that encompassed the region of interest averaged along the plateau. In the specific case dealt here, A2A, the analysis included the binding site and the region around the sodium ion. From the density map, we could have an indication of reluctance, in relative terms, of water molecules to leave any given grid cube, and we could identify the local occupancy maxima. We set the side of the grid to 0.5 Å to provide sufficient resolution. 2.6. Bulk Accessibility Estimation. Water accessibility to the bulk was calculated by stripping explicit water molecules from each frame of the trajectory and calculating the solvent-

Figure 1. Fast desolvation protocol: the collective variable value is halved with respect to its starting value. Upon reduction of the CV, the site can become desolvated by different extents depending on the resilience of the surrounding water interaction network. Several tests with different values of duration of the desolvation part and of the second plateau indicate that 150 and 1000 ps, respectively, can be a reasonable choice (see Supporting Information).

Figure 2. Class-A GPCR conserved waters, organized in clusters as in Angel et al.,38 are mapped to the high-resolution 4EIY structure and compared with occupancy from our apo plain simulation. The average simulated protein structure (cyan) is overlaid on the crystal (pink); crystallographic waters are shown in tan color, and simulation results are represented as water oxygen maxima. (A) Cluster 6, located near the ligand binding pocket, connects Ile80, Phe83, and Val84 in helix III. The 4EIY member of this cluster is W2515 (extracellular cluster) and is crucial to stabilize the kink in helix III. It also mediates interaction with Ala59 and Ala63 of helix II. (B) The Na+ coordination in the “central cluster”, which corresponds to water cluster 3 is shown. (C) One water molecule was found near residues related to water cluster 2, namely, Val45, Ala49, and Tyr288. (D) Water cluster 9 is adjacent to the ligand binding pocket and stabilizes the proline kink of the highly conserved WxP6.50F/Y motif of TM6. A persistent water maximum is observed near W2534. (E) W2532’s position in water cluster 13 is found near a persistent maximum. It interacts with Asn284 of the highly conserved NPxxY motif in helix VII. D


Article


Å of the ion, was consistently 2.7, which is in accordance with the three crystallographic waters that coordinate the ion. Near the coordinating waters, namely, W2529, W2552, and W2579, we found two density maxima, ΓNa8 and ΓNa10, 0.58 and 1.53 Å from W2552 and W2579, respectively. We also found a maximum, ΓNa4, 0.62 Å from W2528. Further analysis of the solvation around the sodium ion reveals the presence of water making persistent hydrogen bonds with three residues, namely, with OD1 of Asp52, OG of Ser91, and OD1 of Asn280. Moreover, the simulation largely agrees, in terms of water− sodium coordination, with Fenalti et al.’s recently resolved crystal structure of the delta-opioid receptor (1.8 Å resolution),69 where the Na+ ion is coordinated by two waters and by the oxygen atoms of two GPCR-conserved residues, Asp952.50 and Ser1353.39, and one nonconserved residue, Asn1313.35, thus maintaining a total coordination number of 5. We also analyzed the water network around Asn24, Asp52, and Asn284 (the contact residues of water cluster 3). We found a pronounced water maximum, Γcl35 (aka Γcl2100), 0.80 Å from W2530. Γcl35 is 3.46 Å from backbone N of Ala49, 3.81 Å from ND2 of Asn284, 2.82 Å from OD1 of Asn24, and 2.98 Å from the OD2 of Asp52. Another water maximum, Γcl32, is 1.44 Å from W2533 in the 4EIY central cluster. Γcl32 is 3.84 Å from the ND2 of Asn24, 3.53 Å from the OD2 of Asp52, 3.23 Å from the backbone oxygen of Asn280, 3.48 Å from ND2, 2.79 Å from the backbone nitrogen of Asn284, and 3.59 Å from OD1 of Asn24. These results are shown in panel B of Figure 2. 3.1.2. Cluster 6. We found a maximum, Γcl625 (Figure 2, panel A), 1.18 Å from W2515, near the hydrophobic residues from Ile80-Val84 of helix III in the 4EIY crystal structure (contact residues of water cluster 6). During the simulation, the persistent water at this position interacts with the backbone oxygen of Ile80 (average distance of 2.97 Å), Ala81 (3.66 Å), and is also located near the backbone nitrogen of Phe83 (3.37 Å), Val84 (2.97 Å), Cys82 (3.89 Å), and Ala81 (4.01 Å). Γcl619 is 4.12 Å from Γcl625, 1.11 Å from W2523, and in possible interaction with O of Ile80 (3.58 Å) and with the backbone N of Ala63 (3.98 Å). A further maximum, Γcl613, is 1.27 Å from W2584. Γcl613 is within interaction distance (2.69 Å) of Γcl69, which in turn is 0.80 Å from W2522 and 4.14 Å from backbone O of Ala59. In the high-resolution crystal structure of 4EIY, we observed a water network that stabilized interactions between helix III and helix II (e.g., via Ala59 and Ala63). 3.1.3. Cluster 9. A clear maximum Γcl925 corresponds to water cluster 9. In 4EIY, only one member of this cluster can be found, W2534, at a distance of 1.67 Å. Panel D of Figure 2 represents the high-density water maximum in this region. This maximum is adjacent to the ligand binding pocket near the residues of Ala273, Cys245, and Leu249. This water is known to stabilize the proline kink of the highly conserved WxP6.50F/Y motif in the sixth TM. However, in our simulation, this maximum is slightly displaced from the Pro248 and Leu249 but is found close to the backbone oxygen of Ala273 and Cys245 residues (at distances of 2.71 and 2.8 Å, respectively). 3.1.4. Cluster 13. We found a well-defined maximum, Γcl137, 0.74 Å from W2532 (Figure 2, panel E). This water belongs to water cluster 13, which stabilizes the conserved NPxxY motif in helix VII. It is close to the OD1 and ND2 of Asn284 (at distances of 2.35 and 3.67 Å, respectively) and also 3.68 Å from the backbone N of Phe242, thus connecting helices VI and VII. A more extended analysis and comparison with existing data in the literature is reported as Supporting Information.

excluded surface (SES), according to the Lee and Richards definition,41 using a probe radius of 1.4 Å with NanoShaper.66 For each frame, it is easy to derive an accessibility map made on the same grid used for persistence estimation. This allows timeaveraging of both quantities along the trajectory. 2.7. Comparison with Hydration Free Energy-Based Methods. As noted by Nguyen et al.,68 correlation between binding affinity and regional hydration free energy in the area covered by the bound ligand can be observed with difficulty. This can be due to sensitivity to the system’s geometry or because the low mobility of some water molecules hinders the convergence of the estimator. In this approach, we try to see whether creating a bias that lowers, or even reverses, the dehydration barrier, can provide useful information on the hydration of a biomolecular system. In this sense, we accelerate the “unbinding” of water molecules from the selected region and use observed exit events to infer local water persistency. Although the proposed method does not require simulating the entire number of dehydration processes that would be needed to get quantitative estimates, we still need to assess whether the information provided by a single run is consistent. To do this, we performed a reversibility analysis on one complex. We calculated a water density map along 1 ns of equilibration. Then, we applied the fast dehydration protocol followed by 900 ps where the bias was kept on. Then, the CV was again brought back to the initial value in 50 ps. At this point, two possibilities have been explored: continuing the simulation for 1 ns after switching off the bias or continuing the simulation for 1 ns while forcing the CV at the initial value. These results are shown in section 3.3.5.

3. RESULTS 3.1. Apo-Unbiased Simulation. By analyzing the water occupancy map derived from our 100 ns plain molecular dynamics simulation, we identified several persistent water molecules in regions corresponding to previously identified and published conserved water clusters38,42 plus the sodium ion surrounding region (see Figure 2). Where possible, we also compared our results with experimental evidence provided by water molecules present in the 4EIY and/or 3EML crystals. Clearly, this was possible only out of the region of the binding site occupied by the ligand because ours is an apo simulation. Results of the following sections, represented in terms of the observed maxima, are gathered in six pdb files, included in the Supporting Information. 3.1.1. Clusters 2 and 3 and Na+ Coordination. Overall, we observed a significantly high probability of finding water molecules in the region around the residues Val45, Ala49, and Tyr288, which has been related to water cluster 2. This suggests a water-mediated structural interaction between helix II and helix VII (Figure 2, panel C). Here, we observed two significant maxima, Γcl2100 and Γcl2101. Γcl2100, which was also observed in the adjacent cluster 3, exhibiting a strong persistence and corresponding to the crystal water W2530 (0.80 Å). Γcl2101 has a lower occupancy of approximately 60% and is within hydrogen bond distance of both Γcl2100 (2.69 Å) and the OH group of Tyr288 in helix VII (2.61 Å) and Γcl2120, located 3.56 Å from the backbone oxygen of Val45. The plain MD shows that the region surrounding the Na+ ion is characterized by a highly dynamic behavior with considerable water interchange. Compatible with the crystal, the average hydration number of the Na+ ion along the simulation, i.e., the average number of water molecules within 4 E


Article


Figure 3. Dehydration plot upon bias application in apo structure simulation. The target value of the CV (in green) is brought from the initial value to 0; however, the actual value of the CV does not closely follow the target and stops at approximately 30% of the initial value. Correspondingly, the number of nearby water molecules, in light blue, decreases but does not fall to zero.

Table 1. Summary of Main Correspondences between Experimental Data and Apo Simulations water cluster ID

4EIY res. 1.8 Å

3EML res. 2.6 Å

Apo (plain) Figure 2

Apobias (desolv.) Figure 4

reference

dist. wrt 4EIY

dist. wrt 4EIY

dist. wrt 4EIY

Γcl2101(2.05 Å) panel C Γcl35(0.80 Å) Γcl32(1.44 Å) panel B Γcl625 (1.18 Å) panel A Γcl925 (1.67 Å) panel D Γcl137 (0.74 Å) panel E

Γbias131(0.57 Å) panel A Γbias10(1.39 Å)

2

W2512

W509(0.28 Å)

3

W2530 W2533 (CC)

W504(0.29 Å)

6

W2515 (EC)

W501(1.44 Å)

9

W2534

W502(2.5 Å)

13

W2532

W512(1.54 Å)

panel B

Γbias9(0.71 Å) panel C

We found water maxima Γbias12 and Γbias13 at distances of 0.40 and 0.98 Å from W2552 and W2579, respectively. The maximum Γbias12, located near W2552, is 2.82 Å from the Na+ ion, whereas the other maximum, Γbias13, is 2.04 Å from Na+ (Figure 4, panel B). The coordination of the Na+ ion with water W2529 of 4EIY was not maintained after the application of the bias. We found a water maximum Γbias10 in the region corresponding to water cluster 3 (Figure 4, panel B). This maximum is 1.39 Å from the crystallographic water W2530, 3.15 Å from the OD2 of Asp52, and 2.36 Å from the OD1 of Asn24. The density maximum, which is observed in the unbiased simulation near crystallographic water W2533, is lost upon dehydration. We also found a maximum, Γbias131, which remarkably overlaps W2512 (cluster 2) at a distance of 0.57 Å (Figure 4, panel A). The contacting residues of water cluster 2 are found close to this persistent water at distances of 2.75, 3.63, and 2.67 Å from the backbone oxygen of Val45, the backbone nitrogen of Ala49, and the OH of Tyr288, respectively. Another maximum we found is Γbias9 (Figure 4, panel C), which is 0.71 Å from W2532. This water belongs to cluster 13 and stabilizes the conserved NPxxY motif in helix VII. It is located 2.48 and 3.64 Å from the OD1 and ND2 of Asn284, respectively, and 3.55 Å from the backbone N of

3.2. Apo-Biased Simulation. In this simulation, we applied the slow dehydration protocol (as explained in the Methods) to the A2A binding pocket (and including the sodium ion in the repulsion group, see Figure 1SI in the Supporting Information). The total amount of simulation time was 28 ns, divided into 14 plateaus, each lasting 2 ns, during which the collective variable value was kept constant. Along the 5th and 11th plateaus, which correspond to the intermediate and the strongest biasing forces, respectively, we computed a new grid for the water persistence calculation. As can be seen in Figure 3, before applying the bias, there were 27.5 ± 1.5 water molecules in the pocket identifying the binding site. As expected, this number decreased at each step due to the repulsive bias, but in a nonlinear way, indicating that the system was reluctant to release waters in proportion to the increment of the bias. As a result, during the final dehydration step, an average of 3.6 water molecules remained in the surroundings of the binding pocket. These molecules, mentioned also in Table 1, belong to the conserved water clusters described in the literature, supporting their potential structural role. As with the unbiased simulation, we analyzed the water molecules in the central region of the apo receptor, where the Na+ ion is located. We observed two water molecules coordinating the Na+ ion during biased dehydration near Γbias12 and Γbias13. F


Article


Figure 4. Apo-biased simulation. The conserved water clusters and the persistent waters during simulation are compared. The protein structure during simulation (cyan) is overlaid on the 4EIY crystal (pink), crystallographic waters are shown in tan color. Local occupancy maxima are represented as sphere color coded by occupancy value. (A) Persistent water maximum during simulation is found near W2512 in water cluster 2. (B) Persistent waters coordinating with the Na+ ion. This region corresponds to the “central cluster” of the 4EIY crystal and is close to water cluster 3. (C) A pronounced maximum was found very close to the crystallographic water W2532 in cluster 13 and stabilizes the conserved NPxxY motif in helix VII.

Phe242, indicating its connecting role between helices VI and VII. Table 1 summarizes the main correspondences between the water molecules in the 4EIY crystal and our simulations, organized in the different clusters of Angel et al.38 Tables 1SI and 2SI in the Supporting Information instead report occupancy data for each maximum of apo-unbiased and -biased simulations. 3.3. Biased Dehydration of A2A in Complex with Triazine Derivatives. When simulating the holo structures, we adopted the fast desolvation protocol as described in the Methods. The repulsive atomic set in our dehydration bias consisted of the heavy atoms of the binding site residues, those of the ligand at hand, and the Na+ ion. The first plateau of 50 ps was followed by a decreasing ramp of 50 ps so that the bias finally reached half of its initial value. Then, an overly long second plateau of 1900 ps was performed to estimate the time needed for the system to relax to the equilibrium at the new CV value for a total simulation time of 2 ns for each complex. Sensitivity analysis, shown in the Supporting Information, concerning the durations of these steps indicates that a desolvation ramp of 150 ps followed by a second plateau of 1 ns is enough to reach convergence. To validate this fast protocol, we selected four members of the series of 1,2,4-triazine derivatives, which inhibit the A2A receptor (see Figure 5). 5,6-diphenyl-1,2,4-triazine-3-amine (4a) served as a scaffold from which the original series of 12 triazine derivatives was obtained.37 They present significant differences in residence times and have been analyzed in Bortolato et al.2 The triazine core is connected to two phenyl rings A and B with variable substituents. Ring A is connected to C6 and ring B is connected to C5 of the triazine core. The selected compounds 4e, 4f, and 4g were developed by chemical modification of the parent (4a) and aimed to optimize the

Figure 5. Triazine-derivative inhibitors series considered in this work.

ligand potency. The modifications consisted of decorating the A ring while leaving the rest of the structure unchanged. They allowed the identification of two lead compounds (4e and 4f). 4e and 4f have common phenolic hydroxyls at their para positions on the A ring. Additionally, 4f has 3,5-disubstituted methyl groups whereas 4e has a chlorine substitution at one of the meta positions of its phenolic hydroxyl ring. 4e and 4g were crystallized in complex with A2A (PDB IDs: 3UZC and 3UZA, respectively). The residence time of the parent compound was extremely short, less than 1 s, whereas that of 4e was 990 s. Replacing the chlorine group in 4e with a 3,5-dimethyl moiety of the phenolic hydroxyl ring in 4f slightly decreased the residence time from 990 to 735s.2 Bortolato et al. have identified two subpockets in the binding site in terms of being G


Article

Journal of Chemical Theory and Computation able to host, or entrap, water molecules.2 Water pocket A (see Figure 2SI in the Supporting Information) is surrounded by Ala59, Phe62, Ala63, Ile66, Ile80, Ala81, Phe83, Val84, Phe168, and His278. Water pocket B is surrounded by Leu85, Ala88, Gln89, Ile92, Ile135, Gly136, Pro139, Tyr176, Met177, Phe180, Asn181, Cys185, Phe242, Trp246, and His250. Below, we present our analysis of each ligand together with a corresponding 3D representation, including the average occupancy map, which is recalculated for each ligand. Hereafter, we consider buried waters to be those that are located in a solvent-excluded region for more than the 80% of the frames of the second plateau. Table 3SI contains the occupancy and accessibility values for all the observed local occupancy maxima. As detailed below, on the basis of this definition and on the results of our simulations, not all water molecules located in pockets A and B can be considered as buried, as expected. 3.3.1. Ligand 4a. As shown in Figure 6, in this conformation, nitrogen N3 of the ligand interacts with OD1

(3.27 Å), and OG1 of Thr256 (2.95 Å). A water molecule in Γa7 is 2.88 Å from the OD1 and 4.06 Å from the backbone N of Asn253. In contrast to the previously mentioned maxima, Γa10, Γa12, and Γa17 are basically never solvent accessible. In particular, Γa10 is located in pocket A facing the phenyl ring. A water here can only mediate the interaction between the side chain NE2 of the doubly protonated His278 and the backbone oxygen of Ala59 (respective distances of 2.45 and 3.15 Å). Similarly, the number of potential interactions of Γa12 is limited (2.64 Å from ND2 of Asn253, 3.39 Å from the backbone oxygen of Met177). The remaining region surrounding Γa12 is basically hydrophobic, suggesting that the affinity might benefit from a modification that displaces it from that location. Finally, a water molecule in Γa17, which is located in pocket B, could interact with NE2 of Gln89, with OD1, and with the backbone oxygen of Asn181, being 2.69, 3.20, and 3.23 Å from them, respectively. Overall, the short residence time and high rate of unbinding of this ligand can be accounted for by the fact that almost only the triazine moiety, which is pointing out of the binding site, is stabilized by water and water-mediated interactions. 3.3.2. Ligand 4e. In the biased simulation of the 4e complex, on the basis of the 3UZC structure, we found some tightly bound water molecules. Moreover, we observed that ligand 4e maintains a well-characterized water network throughout the biased simulation. In particular, the hydroxyl group present in the para position of the phenyl ring is involved in a watermediated interaction network that connects the ligand and several residues of the receptor binding site (Figure 7). In water

Figure 6. Hydration in the 4a complex.

of Asn253 and OE2 of Glu169. In turn, N4 of the ligand’s amino-triazine scaffold interacts with ND2 of Asn253. Overall hydration analysis shows that, unlike for the other compounds described below, almost half of the water molecules located near the ligand in the equilibrated system left their position upon the application of the bias (data not shown). As a result of our protocol, we hypothesize that N3 of the ligand mediates the interaction with the two water molecules identified by the Γa6 and Γa7 maxima located 3.09 and 3.16 Å from N3, respectively. Moreover, the surroundings of the ligand exhibit interesting hydration patterns, which are described below. Waters located at Γa6 and Γa7 are in a favorable position (av distance of 2.24 Å) to make hydrogen bonding. Further interactions can be identified between Γa7 and Γa4 (2.6 Å), between Γa6 and OE2 of Glu169 (2.99 Å), and between Γa4 and Γa6 (3.84 Å). Furthermore, a water molecule in Γa4 can interact with the backbone oxygen of Ile252 (3.92 Å), the backbone nitrogen of Ala265 (3.07 Å), the ND1 of the (doubly protonated) His264

Figure 7. Hydration in the 4e complex.

pocket A, we found two density maxima, Γe6 and Γe11, located in partially and fully buried positions, respectively, that were still favorable for establishing interactions with both the ligand and the protein. Γe6 is 2.68 Å from the hydroxyl group of 4e, 2.46 Å from the backbone oxygen of Ala59, and 3.5 Å from the backbone nitrogen of Ala63. In turn, Γe11 is within bonding distance (3 Å) of the backbone oxygen of Ile80, and of the H


Article


and 2.70 Å, respectively. Γf14 is also connected to Γf16 (2.83 Å). In turn, Γf16 is within bonding distance of the backbone oxygen of Leu85 (3.96 Å), the side chain OH of Tyr176 (3.47 Å), and OD1 of Asn181 (2.54 Å). Three maxima, Γf6, Γf10, and Γf13, are located near the N3 of the ligand’s amino-triazine core with 21, 0, and 50% accessibility values, respectively. Γf6 is 3.34 Å from the N3 of the ligand. It mediates the interaction with OD1 of Asn253 (3.34 Å) and with the backbone oxygen of Leu249 (2.71 Å). Γf6 is also 3.26 Å from the N4 of the ligand’s amino-triazine core. Another maximum, Γf10, is 3.64 Å from the N3 of the amino-triazine core. Γf10 is located between OE1 of Glu169, ND2 of Asn253, and ND1 of His264 at distances of 2.36, 3.14, and 3.52 Å, respectively. Γf10 is within bonding distance of another pronounced maximum, Γf13 (2.69 Å), which is in turn connected with the backbone oxygen atoms of Glu169 and Val172 and the backbone nitrogen of Met174 at distances of 3.06, 2.71, and 3.45 Å, respectively. Two more water maxima, Γf3 and Γf4, have also been observed near the ligand. Γf3 is 2.24 Å from Γf4, which is in turn connected to OE2 of Glu169 at a distance of 3.07 Å. For 4f, we observed a persistent water maximum in pocket A near this ligand’s hydroxyl moiety. The most significant maximum, Γf8, is connected to the oxygen atom of the OH moiety of the 4f ligand at a distance of 3.69 Å. Γf8 is connected to the NE2 of the His278 at a distance of 2.60 Å, backbone oxygen of Ala59 and nitrogen of Ala63 at distances of 2.93 and 4.12 Å, respectively. On the basis of these water networks, we suggest that for 4f (similarly to 4e), the hydroxyl moieties in the para position on their phenolic A ring make a crucial contribution to the experimentally verified binding affinity. Here too we found good agreement overall with Bortolato et al.’s results,2 and a possible reason justifying the fact that, in the analyzed series, 4f’s residence time was slightly shorter than that of 4e. This can be due to the fact that the methyl substituents in 4f interfere with the water network by occupying the position of the water molecule bound to Ile80. 3.3.4. Ligand 4g. In the biased simulation of the 4g complex, on the basis of the 3UZA structure, we observed the following water-mediated interaction between the protein and the ligand. The N3 atom of the ligand’s triazine moiety can interact with water located at the Γg9 occupancy maximum (which has an accessibility around 64%), which is 3.14 Å away. Γg9 is also within bonding distance of the OE2 of Glu169 (2.81 Å) and the OD1 (3.87 Å) and ND2 (2.89 Å) of Asn253. Additionally, water in Γg9 can interact with that in Γg10 (3.54 Å). In turn, a water molecule located in Γg10 can interact with NE2 of His264 (2.81 Å). A water molecule located at Γg10 can also interact with OE2 of Glu169 (2.62 Å). Moreover, OE1 of Glu169 is 3.34 Å from Γg5, a fully solvent accessible maximum that can be connected with another accessible maximum, Γg4 (2.69 Å). Figure 9 depicts the most likely hydration pattern around ligand 4g. In pocket A, we found three maxima, Γg7, Γg8, and Γg11, with average accessibilities of 98, 38, and 100%, respectively (see Table 3SI in the Supporting Information). The nitrogen of 4g’s dimethylpyridine moiety is 3.57 Å from Γg7. Γg7 is 2.49 Å from NE2 of His278 and 2.69 Å from Γg8. Γg8 is 2.52 Å from the backbone oxygen of Ala59 and 2.87 Å from the maximum Γg11. In turn, Γg11 is connected to the backbone oxygen of Ala59, Ile80, and Ala81 and to the backbone nitrogen of Val84, at distances of 3.62, 2.92, 3.05, and 3.71 Å, respectively. On the other side, two buried maxima, Γg19 and Γg16, are located in water pocket B near the Asn181 and Leu85 residues,

backbone nitrogen atoms of Phe83 (3.49 Å) and Val84 (3.41 Å). The unsubstituted phenolic ring B of the ligand is adjacent to water pocket B. We also observed five fully solvent-accessible and persistent water molecules that correspond to the density maxima Γe12, Γe13, Γe14, Γe15, and Γe17 in this hydration site. Γe12 is located near the phenyl ring and is more external, whereas Γe14 and Γe15 are in the interior of the cavity. Γe12 allows a bond with the backbone oxygen atoms of Ala88 (3.85 Å) and Cys185 (3.61 Å). Moreover, Γe12 is within bonding distance of the backbone (3.87 Å) and NE2 (3.94 Å) nitrogen atoms of Gln89. It is 3.5 Å from another well-defined water density maximum, Γe14. Γe14 is located in the innermost part of pocket B and is 2.69 Å from Γe15. Γe14 is 3.11 Å from the backbone oxygen of Leu85, and 3.06 and 4 Å from NE2 and the backbone nitrogen of Gln89, respectively. The other maximum, Γe15, in turn allows a hydrogen bond with OD1 of Asn181 (3.19 Å) and with the OH of Tyr176 (3.59 Å). Via its nitrogen N3, the ligand can interact with a water molecule located in the persistent and solvent-accessible maximum Γe9 (2.92 Å). In turn, a water molecule located in Γe9 can interact with OE2 of Glu169 (2.56 Å), with ND2 of Asn253 (3.04 Å), and also with a 73.9% buried water molecule located at Γe7 (2.69 Å). Finally, Γe7 is within bonding distance of the OE2 of Glu169 (2.07 Å) and the NE2 of the doubly protonated His264 (2.97 Å). The OE1 atom of Glu169 is also 3 Å from Γe7. Γe7 is 3.75 and 3.96 Å from OD1 and ND2 of Asn253, respectively. 3.3.3. Ligand 4f. The water network around the ligand is shown in Figure 8. In water pocket B, we found three maxima

Figure 8. Hydration in the 4f complex.

having accessibility around 50%: Γf12, Γf14, and Γf16. Γf12 and Γf14 are located in front of the unsubstituted phenyl ring, whereas Γf16 is located in the interior of the pocket. Moreover, Γf12 and Γf14, which are near the phenyl ring, are within bonding distance (2.69 Å). Γf12 is also 3.28 Å from ND2 of Asn181. Γf14 is connected with the backbone oxygen of Leu85, OD1 of Asn181, and NE2 of Gln89 at distances of 2.92, 3.88, I


Article


literature,3 which also describe the chain of these three waters in pocket A as highly stable.2,3 Further analysis concerning identified positions and possible interactions is provided as Supporting Information. 3.3.5. Effect of Restraining and Reversibility Analysis. RMSD analysis of the protein upon application of the fast dehydration protocol with a restraint of 10,000 kJ mol−1 nm−2, which is weaker than that used in similar approaches,31 provides an average RMSD value for the protein of 0.14 Å with a standard deviation of 0.01 Å. Reversibility analysis was performed as described in section 2.7 on the 4e complex. As shown in Figure 3SI, 14 maxima (with occupancy >80%) were found during the initial plain trajectory. After a fast dehydration, four positions remained unaffected, four were displaced, and six disappeared. When the bias was removed, three of the displaced maxima spontaneously regained the original positions as well as four of the ones that previously disappeared, for a total of 11 recovered maxima (>78%). In the last phase, when the CV was forced to take the original average value, 13 maxima were recovered (∼80% or higher), and one was displaced. Results are collected in Table 2 and reported graphically in Figures 3SI and 4SI in the Supporting Information.

4. DISCUSSION AND CONCLUSIONS In this work, we characterized the hydration patterns of the apo and 4 holo forms of the A2A receptor using a novel moleculardynamics-based approach, which uses a bias to facilitate dehydration around a given region of a molecular system and uses the reaction of water molecules as an immediate measure of their persistence. The peculiarity of the method is rooted in the definition of an electrostatic-like bias that aims at forcing a physically sound desolvation of the site of interest. The form of the bias, which is however not unique, takes into account the inherently collective character of solvation and desolvation processes and the possible geometric complexity of the region to be analyzed. Other existing computational methods try to provide a quantitative estimate, in terms of free energy, of the impact of the displacement of individual water molecules. Achieving this task is perhaps not essential in the drug design context and can

Figure 9. Hydration in the 4g complex.

which is in line with Bortolato et al.’s results.2 Five waters in this position have been described by Bortolato et al.,2 who classified three of those waters as “unstable” or “unhappy”. This might be because they are far from 4g’s phenyl moiety and the unsubstituted phenyl ring lacks any favorable interaction with them. Our simulation found these waters to be buried in the pocket. Figure 9 shows the water maxima in the surroundings of 4g. The interaction of waters with the ligand and the receptor helps rationalize 4g’s off-rate, which is approximately 1 order of magnitude faster than that of the congeneric compound 4e. These results are also in good agreement with the stable water clusters recently identified by water fluid dynamics (WFD) maps and with other sources in the Table 2. Reversibility Analysis Resultsa

a

The maximum marked with a star (*) was still observed but only with a persistency of approximately 53%. J


Article


decorating ring B of the parent compound and pointing toward water-pocket B (namely, 4h, 4i, 4j, and 4l). The proposed analysis also suggests a possible general strategy for building improved agonists/antagonists for GPCRs. The exploration of interhelical spaces (or near-helical spaces) seems promising for further expanding a ligand scaffold to displace waters or connect to extended water networks. Our results support the opinion that optimizing a ligand in light of the impact on hydration is subtler but can be as important as optimizing the direct interaction of the target protein. In the context of ligand structure−kinetics relationships, this approach can be paired to other computational methods that try to rank congeneric ligands based on residence time.70 Overall, the proposed method seems to have good potential, offering a quick dynamical approach to water behavior and an intuitive way to suggest directions where to grow a lead compound. Moreover, it does not involve training on a data set to tune its parameters, which are minimal. This approach is one example of how MD can be of use in the drug design context71 and is in total accordance with the conclusions of a recently published review of scientists from Roche Pharmaceutical Research,72 who advocate the adoption in practical drug discovery of qualitative approaches that try to describe a process in terms as simple as possible and to provide hints to the medicinal chemistry teams.

be far from trivial, for instance, when the analyzed region contains partially buried water molecules for which estimating statistical mechanical figures in an MD-based study is possibly hindered by lack of convergence. In contrast, the aim of the suggested approach is to simply identify regions where the water network is easier to disrupt, assuming that they are the most advantageous for growing the candidate in the hit-to-lead and lead optimization phases. These regions emerge as the first that become depleted upon application of the bias. This simplification allows an implementation that is quite fast with respect to other MD-based methods, is numerically stable, and provides results with easy interpretation, see for instance the water depletion movie in the Supporting Information. In fact, the originally envisioned simulated length of the dehydration protocol, of 2 ns turned out to be longer than needed, leading to a shorter protocol of 1.2 ns. Particular attention is paid to the presence of water molecules that become buried upon ligand binding, which might be persistent only because of steric reasons. These cases are analyzed by considering the interactions that these waters are making and could be integrated by a comparison with the corresponding apo-biased simulation. If persistent and buried water is located in a region that was depleted in the apo-biased simulation, it is reasonable to assume that it is persistent but “unhappy”, especially in the case of water molecules displaying significant interactions with hydrophobic moieties. Interestingly, buried and “unhappy” water molecules have been pointed to as particularly responsible for shorter ligand residence times.2 Some analysis along this line is presented in Table 4SI of the Supporting Information. The validation of the protocol, made on the apo A2A structure, was conducted by probing waters located in several sites that have been identified in other class A GPCRs. Interestingly, by gradually increasing the dehydrating bias, we could single out a few water molecules that have a recognized structural role (Table 1). Furthermore, the application of the validated protocol on congeneric inhibitors of the A2A receptor provided results in substantial agreement with recent literature correlating the presence of “happy” and “unhappy” waters to the residence times of triazine-derivative A2A inhibitors.2,3 Summarizing the analysis on the studied complexes, in the 4a case, we observed only one highly persistent but buried water in each of the pockets identified by Bortolato. In the 4e ligand complex, we observed two persistent and mostly buried water maxima in pocket A and five highly persistent and accessible water molecules in pocket B. In this respect, our study integrates the work of Bortolato et al., suggesting that the bulk accessibility of the two pockets is also a dynamical feature. The trapped waters in pocket A mediate favorable interactions between the hydroxyl moiety of the 4e ligand and the pocket residues. In 4f, we observed three persistent but partially solvent-accessible water maxima in pocket B and one persistent and buried maximum, Γf8, in pocket A. In 4g, we observed three persistent maxima, with two of them being fully accessible, in pocket A, and two buried and highly persistent maxima in pocket B. These observations are collected in Table 3SI of the Supporting Information. The modifications of the compounds considered in this work involve ring A; however, our analysis indicates the presence of persistent and accessible waters in pocket B. This is in agreement with the reduction of affinity and residence time of ligands (9, 11, 23, and 24 s, respectively) obtained by

■

ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jctc.6b00475. Extended Analysis; Table 1SI, maxima observed in the Apo-unbiased simulation; Table 2SI, maxima observed in the Apo biased simulations; Figure 1SI, selected residues for the application of the repulsive bias during dehydration of A2A; Figure 2SI, identification of water pockets A and B; Table 3SI, main maxima observed in the 4 holo-biased simulations; Table 4SI, summary of the findings observed in the 4 holo-biased simulations; Figure 3SI, results of reversibility analysis I; Figure 4SI, results of reversibility analysis II; Figure 5SI, statistical robustness and sensitivity analysis; Table 5SI, sensitivity analysis for the duration of the dehydration ramp (PDF) PDB file containing apo A2A structure and occupancy maxima obtained along plain MD simulation (PDB) PDB file containing apo A2A structure and occupancy maxima after maximum dehydration (lowest plateau) (PDB) Movie illustrating the dehydration process for the 4e holo case giving an intuitive picture of where less stable waters are located (MPG) PDB file containing the 4a holo structure and the occupancy maxima after the application of the fast dehydration protocol (PDB) PDB file containing the 4e holo structure and the occupancy maxima after the application of the fast dehydration protocol (PDB) PDB file containing the 4f holo structure and the occupancy maxima after the application of the fast dehydration protocol (PDB) K


Article


■

PDB file containing the 4g holo structure and the occupancy maxima after the application of the fast dehydration protocol (PDB) Water occupancy maxima in the pdb files are given a residue name that corresponds to the water cluster number as identified by Angel et. al.

(10) Mancera, R. L. Molecular Modeling of Hydration in Drug Design. Curr. Opin. Drug Discovery Devel. 2007, 10 (3), 275−280. (11) Wong, S. E.; Lightstone, F. C. Accounting for Water Molecules in Drug Design. Expert Opin. Drug Discovery 2011, 6 (1), 65−74. (12) Young, T.; Abel, R.; Kim, B.; Berne, B. J.; Friesner, R. A. Motifs for Molecular Recognition Exploiting Hydrophobic Enclosure in Protein-Ligand Binding. Proc. Natl. Acad. Sci. U. S. A. 2007, 104 (3), 808−813. (13) Bissantz, C.; Kuhn, B.; Stahl, M. A Medicinal Chemist’s Guide to Molecular Interactions. J. Med. Chem. 2010, 53 (14), 5061−5084. (14) Baron, R.; Setny, P.; McCammon, J. A. Water in Cavity-Ligand Recognition. J. Am. Chem. Soc. 2010, 132 (34), 12091−12097. (15) Riniker, S.; Barandun, L. J.; Diederich, F.; Krämer, O.; Steffen, A.; van Gunsteren, W. F. Free Enthalpies of Replacing Water Molecules in Protein Binding Pockets. J. Comput.-Aided Mol. Des. 2012, 26 (12), 1293−1309. (16) Abel, R.; Young, T.; Farid, R.; Berne, B. J.; Friesner, R. A. Role of the Active-Site Solvent in the Thermodynamics of Factor Xa Ligand Binding. J. Am. Chem. Soc. 2008, 130 (9), 2817−2831. (17) Abel, R.; Wang, L.; Friesner, R. A.; Berne, B. J. A DisplacedSolvent Functional Analysis of Model Hydrophobic Enclosures. J. Chem. Theory Comput. 2010, 6 (9), 2924−2934. (18) Hummer, G. Molecular Binding: Under Water’s Influence. Nat. Chem. 2010, 2 (11), 906−907. (19) Siebert, X.; Hummer, G. Biochemistry 2002, 41, 2956. (20) Pearlstein, R. A.; Hu, Q.-Y.; Zhou, J.; Yowe, D.; Levell, J.; Dale, B.; Kaushik, V. K.; Daniels, D.; Hanrahan, S.; Sherman, W.; Abel, R. New Hypotheses about the Structure-Function of Proprotein Convertase Subtilisin/kexin Type 9: Analysis of the Epidermal Growth Factor-like Repeat A Docking Site Using WaterMap. Proteins: Struct., Funct., Genet. 2010, 78 (12), n/a−n/a. (21) Sun, H.; Zhao, L.; Peng, S.; Huang, N. Incorporating Replacement Free Energy of Binding-Site Waters in Molecular Docking. Proteins: Struct., Funct., Genet. 2014, 82 (9), 1765−1776. (22) Lemmon, G.; Meiler, J. Towards Ligand Docking Including Explicit Interface Water Molecules. PLoS One 2013, 8 (6), e67536. (23) SZMAP 1.2.0.7: OpenEye Scientific Software, Santa Fe, NM. http://www.eyesopen.com, 2013; http://docs.eyesopen.com/szmap/ citation.html (accessed Mar 24, 2015). (24) Baroni, M.; Cruciani, G.; Sciabola, S.; Perruccio, F.; Mason, J. S. A Common Reference Framework for Analyzing/comparing Proteins and Ligands. Fingerprints for Ligands and Proteins (FLAP): Theory and Application. J. Chem. Inf. Model. 2007, 47 (2), 279−294. (25) Michel, J.; Tirado-Rives, J.; Jorgensen, W. L. Prediction of the Water Content in Protein Binding Sites. J. Phys. Chem. B 2009, 113 (40), 13337−13346. (26) Michel, J.; Tirado-Rives, J.; Jorgensen, W. L. Energetics of Displacing Water Molecules from Protein Binding Sites: Consequences for Ligand Optimization. J. Am. Chem. Soc. 2009, 131 (42), 15403−15411. (27) Ross, G. A.; Bodnarchuk, M. S.; Essex, J. W. Water Sites, Networks, And Free Energies with Grand Canonical Monte Carlo. J. Am. Chem. Soc. 2015, 137 (47), 14930−14943. (28) Sciabola, S.; Stanton, R. V.; Mills, J. E.; Flocco, M. M.; Baroni, M.; Cruciani, G.; Perruccio, F.; Mason, J. S. High-Throughput Virtual Screening of Proteins Using GRID Molecular Interaction Fields. J. Chem. Inf. Model. 2010, 50 (1), 155−169. (29) Higgs, C.; Beuming, T.; Sherman, W. Hydration Site Thermodynamics Explain SARs for Triazolylpurines Analogues Binding to the A2A Receptor. ACS Med. Chem. Lett. 2010, 1 (4), 160−164. (30) Li, Z.; Lazaridis, T. Computing the Thermodynamic Contributions of Interfacial Water. Methods Mol. Biol. 2012, 819, 393−404. (31) Nguyen, C. N.; Cruz, A.; Gilson, M. K.; Kurtzman, T. Thermodynamics of Water in an Enzyme Active Site: Grid-Based Hydration Analysis of Coagulation Factor Xa. J. Chem. Theory Comput. 2014, 10 (7), 2769−2780.

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Author Contributions

S.R.Z. performed the calculations and the analysis; R.G. and S.D. performed some of the analyses. S.D. implemented the method. W.R. invented the method and designed the experiments. The manuscript was written through contributions of all authors. All authors have given approval to the final version of the manuscript. Notes

The authors declare the following competing financial interest(s): Dr. Decherchi is with BiKi Technologies, a company providing software solutions for Drug Discovery. Some of the techniques described in this work are part of these solutions.

■

ACKNOWLEDGMENTS We acknowledge PRACE for awarding us access to the computing resource at FERMI based at CINECA, Italy. We also thank Dr. Andrea Spitaleri and Prof. Andrea Cavalli for useful discussions.

■

REFERENCES

(1) De Beer, S. B. a; Vermeulen, N. P. E.; Oostenbrink, C.; Oostenbrick, C. The Role of Water Molecules in Computational Drug Design. Curr. Top. Med. Chem. 2010, 10 (1), 55−66. (2) Bortolato, A.; Tehan, B. G.; Bodnarchuk, M. S.; Essex, J. W.; Mason, J. S. Water Network Perturbation in Ligand Binding: Adenosine A(2A) Antagonists as a Case Study. J. Chem. Inf. Model. 2013, 53 (7), 1700−1713. (3) Sabbadin, D.; Ciancetta, A.; Moro, S. Perturbation of Fluid Dynamics Properties of Water Molecules during G Protein-Coupled Receptor-Ligand Recognition: The Human A2A Adenosine Receptor as a Key Study. J. Chem. Inf. Model. 2014, 54 (10), 2846−2855. (4) Poornima, C. S.; Dean, P. M. Hydration in Drug Design. 1. Multiple Hydrogen-Bonding Features of Water Molecules in Mediating Protein-Ligand Interactions. J. Comput.-Aided Mol. Des. 1995, 9 (6), 500−512. (5) Poornima, C. S.; Dean, P. M. Hydration in Drug Design. 3. Conserved Water Molecules at the Ligand-Binding Sites of Homologous Proteins. J. Comput.-Aided Mol. Des. 1995, 9 (6), 521− 531. (6) Ladbury, J. E. Just Add Water! The Effect of Water on the Specificity of Protein-Ligand Binding Sites and Its Potential Application to Drug Design. Chem. Biol. 1996, 3 (12), 973−980. (7) Li, Z.; Lazaridis, T. Thermodynamic Contributions of the Ordered Water Molecule in HIV-1 Protease. J. Am. Chem. Soc. 2003, 125 (22), 6636−6637. (8) García-Sosa, A. T.; Firth-Clark, S.; Mancera, R. L. Including Tightly-Bound Water Molecules in de Novo Drug Design. Exemplification through the in Silico Generation of poly(ADPRibose)polymerase Ligands. J. Chem. Inf. Model. 2005, 45 (3), 624− 633. (9) Li, Z.; Lazaridis, T. Thermodynamics of Buried Water Clusters at a Protein-Ligand Binding Interface. J. Phys. Chem. B 2006, 110 (3), 1464−1475. L


Article

Journal of Chemical Theory and Computation (32) Czapiewski, D.; Zielkiewicz, J. Structural Properties of Hydration Shell around Various Conformations of Simple Polypeptides. J. Phys. Chem. B 2010, 114 (13), 4536−4550. (33) Haider, K.; Huggins, D. J. Combining Solvent Thermodynamic Profiles with Functionality Maps of the Hsp90 Binding Site to Predict the Displacement of Water Molecules. J. Chem. Inf. Model. 2013, 53 (10), 2571−2586. (34) Sun, X.; Ågren, H.; Tu, Y. Functional Water Molecules in Rhodopsin Activation. J. Phys. Chem. B 2014, 118 (37), 10863−10873. (35) Setny, P. Prediction of Water Binding to Protein Hydration Sites with a Discrete, Semiexplicit Solvent Model. J. Chem. Theory Comput. 2015, 11 (12), 5961−5972. (36) Gerogiokas, G.; Southey, M. W. Y.; Mazanetz, M. P.; Hefeitz, A.; Bodkin, M.; Law, R. J.; Michel, J. Correction: Evaluation of Water Displacement Energetics in Protein Binding Sites with Grid Cell Theory. Phys. Chem. Chem. Phys. 2015, 17 (24), 16213−16213. (37) Congreve, M.; Andrews, S. P.; Doré, A. S.; Hollenstein, K.; Hurrell, E.; Langmead, C. J.; Mason, J. S.; Ng, I. W.; Tehan, B.; Zhukov, A.; Weir, M.; Marshall, F. H. Discovery of 1,2,4-Triazine Derivatives as Adenosine A(2A) Antagonists Using Structure Based Drug Design. J. Med. Chem. 2012, 55 (5), 1898−1903. (38) Angel, T. E.; Chance, M. R.; Palczewski, K. Conserved Waters Mediate Structural and Functional Activation of Family A (rhodopsinLike) G Protein-Coupled Receptors. Proc. Natl. Acad. Sci. U. S. A. 2009, 106 (21), 8555−8560. (39) Yuan, S.; Filipek, S.; Palczewski, K.; Vogel, H. Activation of GProtein-Coupled Receptors Correlates with the Formation of a Continuous Internal Water Pathway. Nat. Commun. 2014, 5, 4733. (40) Yuan, S.; Hu, Z.; Filipek, S.; Vogel, H. W246 6.48 Opens a Gate for a Continuous Intrinsic Water Pathway during Activation of the Adenosine A 2A Receptor. Angew. Chem. 2015, 127 (2), 566−569. (41) Lee, B.; Richards, F. M. The Interpretation of Protein Structures: Estimation of Static Accessibility. J. Mol. Biol. 1971, 55 (3), 379−IN4. (42) Liu, W.; Chun, E.; Thompson, A. A.; Chubukov, P.; Xu, F.; Katritch, V.; Han, G. W.; Roth, C. B.; Heitman, L. H.; IJzerman, A. P.; Cherezov, V.; Stevens, R. C. Structural Basis for Allosteric Regulation of GPCRs by Sodium Ions. Science 2012, 337, 232−236. (43) Ballesteros, J. A.; Weinstein, H. Integrated Methods for the Construction of Three-Dimensional Models and Computational Probing of Structure-Function Relations in G Protein-Coupled Receptors. Methods Neurosci. 1995, 25 (C), 366−428. (44) Zhang, X. C.; Sun, K.; Zhang, L.; Li, X.; Cao, C. GPCR Activation: Protonation and Membrane Potential. Protein Cell 2013, 4, 747−760. (45) Sali, A.; Blundell, T. L. Comparative Protein Modelling by Satisfaction of Spatial Restraints. J. Mol. Biol. 1993, 234 (3), 779−815. (46) Eswar, N.; Webb, B.; Marti-Renom, M. A.; Madhusudhan, M. S.; Eramian, D.; Shen, M.-Y.; Pieper, U.; Sali, A. Comparative Protein Structure Modeling Using Modeller. Curr. Protoc. Bioinformatics 2006, 5.6.1. (47) Gordon, J. C.; Myers, J. B.; Folta, T.; Shoja, V.; Heath, L. S.; Onufriev, A. H++: A Server for Estimating pKas and Adding Missing Hydrogens to Macromolecules. Nucleic Acids Res. 2005, 33 (Web Server issue), W368−W371. (48) Myers, J.; Grothaus, G.; Narayanan, S.; Onufriev, A. A Simple Clustering Algorithm Can Be Accurate Enough for Use in Calculations of pKs in Macromolecules. Proteins: Struct., Funct., Genet. 2006, 63 (4), 928−938. (49) Anandakrishnan, R.; Aguilar, B.; Onufriev, A. V. H++ 3.0: Automating pK Prediction and the Preparation of Biomolecular Structures for Atomistic Molecular Modeling and Simulations. Nucleic Acids Res. 2012, 40 (Web Server issue), W537−W541. (50) Fogolari, F.; Zuccato, P.; Esposito, G.; Viglino, P. Biomolecular Electrostatics with the Linearized Poisson-Boltzmann Equation. Biophys. J. 1999, 76 (1Pt 1), 1−16. (51) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25 (9), 1157−1174.

(52) Wang, J.; Cieplak, P.; Kollman, P. A. How Well Does a Restrained Electrostatic Potential (RESP) Model Perform in Calculating Conformational Energies of Organic and Biological Molecules? J. Comput. Chem. 2000, 21 (12), 1049−1074. (53) Vanquelef, E.; Simon, S.; Marquant, G.; Garcia, E.; Klimerak, G.; Delepine, J. C.; Cieplak, P.; Dupradeau, F.-Y. R.E.D. Server: A Web Service for Deriving RESP and ESP Charges and Building Force Field Libraries for New Molecules and Molecular Fragments. Nucleic Acids Res. 2011, 39 (Web Server issue), W511−W517. (54) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4 (3), 435−447. (55) Bussi, G.; Donadio, D.; Parrinello, M. Canonical Sampling through Velocity Rescaling. J. Chem. Phys. 2007, 126 (1), 014101. (56) Parrinello, M.; Rahman, A. Crystal Structure and Pair Potentials: A Molecular-Dynamics Study. Phys. Rev. Lett. 1980, 45 (14), 1196− 1199. (57) Parrinello, M. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52 (12), 7182. (58) Parrinello, M. Strain Fluctuations and Elastic Constants. J. Chem. Phys. 1982, 76 (5), 2662. (59) Yukawa, H. On the Interaction of Elementary Particles. I *. Prog. Theor. Phys. Suppl. 1955, 1, 1−10. (60) Fogolari, F.; Brigo, A.; Molinari, H. The Poisson-Boltzmann Equation for Biomolecular Electrostatics: A Tool for Structural Biology. J. Mol. Recognit. 2002, 15 (6), 377−392. (61) Tuinier, R. Approximate Solutions to the Poisson−Boltzmann Equation in Spherical and Cylindrical Geometry. J. Colloid Interface Sci. 2003, 258 (1), 45−49. (62) Rocchia, W. Poisson-Boltzmann Equation Boundary Conditions for Biological Applications. Math. Comput. Model. 2005, 41 (10), 1109−1118. (63) Debye, P.; Huckel, E. The Theory of Electrolytes. I. Lowering of Freezing Point and Related Phenomena. Phys. Zeitschrift 1923, 24, 185−206. (64) Debye, P.; Huckel, E. Theory of Electrolytes-Part II: Law of the Limit of Electrolytic Conduction. Phys. Zeitschrift 1923, 24, 305−325. (65) Tribello, G. A.; Bonomi, M.; Branduardi, D.; Camilloni, C.; Bussi, G. PLUMED 2: New Feathers for an Old Bird. Comput. Phys. Commun. 2014, 185 (2), 604−613. (66) Decherchi, S.; Rocchia, W. A General and Robust Ray-CastingBased Algorithm for Triangulating Surfaces at the Nanoscale. PLoS One 2013, 8 (4), e59744. (67) Decherchi, S.; Rocchia, W. Building and Analyzing Molecular Surfaces: A Tutorial on NanoShaper. In Computational Electrostatics for Biological Applications; Springer International Publishing: Switzerland, 2015. (68) Nguyen, C. N.; Young, T. K.; Gilson, M. K. Grid Inhomogeneous Solvation Theory: Hydration Structure and Thermodynamics of the Miniature Receptor cucurbit[7]uril. J. Chem. Phys. 2012, 137 (4), 044101. (69) Fenalti, G.; Giguere, P. M.; Katritch, V.; Huang, X.-P.; Thompson, A. a; Cherezov, V.; Roth, B. L.; Stevens, R. C. Molecular Control of Δ-Opioid Receptor Signalling. Nature 2014, 506 (7487), 191−196. (70) Mollica, L.; Decherchi, S.; Zia, S. R.; Gaspari, R.; Cavalli, A.; Rocchia, W. Kinetics of Protein-Ligand Unbinding via Smoothed Potential Molecular Dynamics Simulations. Sci. Rep. 2015, 5, 11539. (71) De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. The Role of Molecular Dynamics and Related Methods in Drug Discovery. J. Med. Chem. 2016, 59, 4035. (72) Kuhn, B.; Guba, W.; Hert, J.; Banner, D. W.; Bissantz, C.; Ceccarelli, S. M.; Haap, W.; Körner, M.; Kuglstatter, A.; Lerner, C. D.; Mattei, P.; Neidhart, W.; Pinard, E.; Rudolph, M. G.; Schulz-Gasch, T.; Woltering, T. J.; Stahl, M. A Real-World Perspective on Molecular Design. J. Med. Chem. 2016, 59, 4087.

M


Probing Hydration Patterns in Class-A GPCRs via ... - ACS Publications

Recommend Documents