Molecular Dynamics Simulation at High Sodium ... - ACS Publications

Mar 2, 2010 - Concentration: Toward the Inactive Conformation of the. Human Adenosine A2A Receptor. Balázs Jójárt,*. ,†. Róbert Kiss,. ‡. Béla Viskolc...
0 downloads 0 Views 4MB Size
pubs.acs.org/JPCL

Molecular Dynamics Simulation at High Sodium Chloride Concentration: Toward the Inactive Conformation of the Human Adenosine A2A Receptor )

,^  j bert Kiss,‡ B €rgy M. Keseru Bal azs Jo art,*,† Ro ela Viskolcz,† Istv an Kolossv ary,§,# and Gyo ? †

)

Department of Chemical Informatics, Faculty of Education, University of Szeged, Boldogasszony sgt. 6, H-6725 Szeged, Hungary, ‡Chemical Engineering, School of Engineering and Physical Sciences, Heriot-Watt University, EH14 4AS Edinburgh, ert t er 4, H-1111 Budapest, U.K., §Department of Chemistry, Budapest University of Technology and Economics, Szt. Gell Hungary, Discovery Chemistry, Gedeon Richter Plc, POB 27, H-1475 Budapest, Hungary, and ^Department of General and Analytical Chemistry, Budapest University of Technology and Economics, Szt. Gell ert t er 4, H-1111 Budapest, Hungary

ABSTRACT The recently solved crystal structure of the human adenosine A2A receptor (hA2AR) shows the characteristics of a partially activated state. Experimental data suggests that high sodium chloride concentration shifts hA2AR to the inactive state. We found that molecular dynamics simulations at high sodium chloride concentration result in an inactive form of hA2AR reflected in the reformation of the “ionic lock” (Arg102(3.50)-Glu228(6.30)) as well as in the reduction of the RC-RC distance between the intracellular sides of transmembrane helices 3 and 6 (TM3 and TM6). Interestingly, no such stabilization effect was observed at physiological concentrations. Our results suggest that the effect of high sodium chloride concentration might be exploited to generate an inactive state of hA2AR, which is more favorable for identifying pharmacologically relevant antagonists or inverse agonists. SECTION Biophysical Chemistry

H

igh-quality structural information about G-protein coupled receptors (GPCRs) is highly demanded in drug discovery. Publication of the first crystal structures of bovine rhodopsin (BR) in the inactive state by Palcziewsky et al.1 was clearly a breakthrough in this field, allowing the development of GPCR homology models. Kobilka et al. reported the first structures of the human β2 adrenergic receptor (hβ2AR) cocrystallized with the partial inverse agonist carazolol.2 Shortly after the publication of the hβ2AR structures, the turkey β1 adrenergic receptor in complex with the antagonist cyanopindolol3 and the human adenosine A2A receptor (hA2AR) in complex with the antagonist ZM2413854 were crystallized. These novel GPCR structures have an invaluable impact on GPCR research, especially because they share significantly higher homology to several therapeutically privileged receptors than BR. Mobarec et al. analyzed the applicability of the novel GPCR structures for homology modeling.5 These authors concluded that the hβ2AR and hA2AR structures are suitable templates for 16% and 12% of all class A GPCRs, respectively, while BR is favorable in only 2% of the same class. The hA2AR structure is in a partially activated state indicated by its in vitro pharmacology.4 Structural characteristics such as the missing “ionic lock” (Arg102(3.50)Glu228(6.30) (Ballesteros-Weinstein numbering scheme6) and the elongated RC-RC distance between these residues (9.7 Å vs 8.0-9.1 Å in BR) also support this hypothesis. This makes the hA2AR crystal structure less suitable for designing

r 2010 American Chemical Society

antagonists, which might possess higher therapeutic potential than agonists. Experimental data suggest that binding of agonists and antagonists to adenosine receptors is modulated in a reciprocal manner by monovalent cations.4,7,8 In the case of hA2AR, it was found that high concentration (1.00 M) of sodium chloride shifts the active/inactive conformational equilibrium of hA2AR to the inactive state.4 We hypothesized that this effect could be exploited to develop an inactive hA2AR structure by molecular dynamics (MD) simulations at high sodium chloride concentration. In this paper, we report our results obtained from 2  3 independent, 70.08 ns long (Σt = 0.42 μs) MD simulations on hA2AR in an explicit membrane environment at 0.15 M (c1) and 1.00 M (c2) NaCl concentrations. During class “A” GPCR activation, the structural alterations confirmed by experimental measurements, are as follows: (1) the disruption of the ionic lock between residues Arg(3.50) and Glu(6.30), and (2) the separation of the intracellular (IC) sides of transmembrane helices 3 and 6 (TM3 and TM6).9 These preliminary experiments were verified by the newly solved crystal structures of opsin (metarhodopsin III state)

Received Date: December 15, 2009 Accepted Date: February 24, 2010 Published on Web Date: March 02, 2010

1008

DOI: 10.1021/jz900403a |J. Phys. Chem. Lett. 2010, 1, 1008–1013

pubs.acs.org/JPCL

Figure 2. Disrupted (a) and intact (b) ionic lock as shown at 30 ns of trajectories ID#2 (c1) and ID#6 (c2). Distance between the RC atoms and between possible bridge atoms of H-bond is depicted via dashed red and blue lines, respectively. In the receptor structure, the TM3 and TM6 domains are highlighted, and the secondary structure elements are indicated by different coloring.13

between the RC atoms of Arg(3.50) and Glu(6.30) and the minimum distance between the guanidine nitrogen atoms of Arg(3.50) and the carboxylate oxygen atoms of Glu(6.30) are indicative for the activate/inactivate state of hβ2AR.12 We therefore monitored the possible inactivation process of hA2AR by means of these parameters (Figure 1, Figure 2). Based on the suggested criterion of Shaw et al.,12 we considered a structure to be rather in the inactive form when the RC-RC distance (d(RC-RC)) was below 9.5 Å. The inactive state could be further stabilized by the presence of the ionic lock, which was monitored by the distance between the guanidine nitrogen atoms of Arg102(3.50) and the carboxylate oxygen atoms of Glu228(6.30) (d(min(O-N)), Figure 2). Accordingly, the individual snapshots of the trajectories were divided into four groups: case1: d(RC-RC) < 9.5 Å and d(min(O-N)) < 3.5 Å; case2: d(RC-RC) < 9.5 Å and d(min(O-N))>3.5 Å; case3: d(RC-RC)>9.5 Å and d(min(O-N))< 3.5 Å; case4: d(RC-RC) > 9.5 Å and d(min(O-N)) > 3.5 Å. More than 90% of the structures in the 1.00 M NaCl trajectories showed an RC-RC distance below 9.5 Å being indicative of the inactive state (see Table 1). Moreover, 71, 70, and 52% of the snapshots contained both a short RC-RC distance and intact ionic lock. The average d(RC-RC) values obtained in the presence of 1.00 M NaCl (8.3-8.7 Å) are in good agreement with distances found in BR crystal structures (8.0-9.1 Å).12 On the other hand, no such stabilization effect was detected in the case of 0.15 M NaCl concentration: the lifetime of the longest inactive state was approximately 400-500 ps only. The average RC-RC distances (9.8911.2 Å) varied among these simulations in agreement with the fact that the active and inactive forms of the GPCR apo structure are in equilibrium at normal salt concentration.14 In simulation ID#2, more than 99% of the snapshots showed

Figure 1. The alteration of the RC-RC distance (light red) and the minimum distance between H-bond forming atoms of Arg102(3.50) and Glu228(6.30) (light blue) in the course of the 0.15 M NaCl (ID#1, ID#2 and ID#3) and 1.00 M NaCl (ID#4, ID#5 and ID#6) simulations. Average values calculated over 98 ps are depicted with dark red and dark blue.

showing no interaction between the above-mentioned residues as well as an extended helical conformation of TM6 in the IC side.10,11 Recently, Shaw et al. showed that the distance

r 2010 American Chemical Society

1009

DOI: 10.1021/jz900403a |J. Phys. Chem. Lett. 2010, 1, 1008–1013

pubs.acs.org/JPCL

Table 1. Distribution of the Snapshots between the Four Groupsa

ID#1 ID#2

case1

case2

case3

case4

average d(RC-RC) ( SD /Å

0.06 0.00

27.87 0.82

0.01 0.00

72.07 99.18

9.90 ( 0.75 11.22 ( 0.50

ID#3

0.74

36.64

0.04

62.58

9.89 ( 0.88

ID#4

71.06

28.10

0.02

0.82

8.31 ( 0.46

ID#5

70.27

29.41

0.00

0.32

8.26 ( 0.46

ID#6

52.57

39.09

0.36

7.98

8.68 ( 0.58

a The following criteria were used to define groups: case1: d(RC-RC)< 9.5 Å and d(min(O-N)) < 3.5 Å; case2: d(RC-RC) < 9.5 Å and d(min(O-N)) > 3.5 Å; case3: d(RC-RC) > 9.5 Å and d(min(O-N)) < 3.5 Å; case4: d(RC-RC)>9.5 Å and d(min(O-N))>3.5 Å. The average RC-RC distances were also calculated for the last 65.28 ns of the trajectories (the first 4.8 ns was considered as equilibration).

higher RC-RC distance than 9.5 Å, and the ionic lock was disrupted during the whole simulation. The average d(RCRC) value in this trajectory (11.2 Å) is in the range found in the recently solved, partially activated β2AR and β1AR structures (11.0-11.2 Å) and far from those of both the inactive BR (8.0-9.1 Å) and the opsin structures (14.2 Å).12 We therefore suggest that snapshots from this trajectory represent a partially activated state of hA2AR. In the case of low salt concentration (see trajectories ID#1, ID#2, and ID#3 in Figure 1), we found that Arg102(3.50) and Glu228(6.30) were unable to form a stable salt bridge in either of the three trajectories. On the other hand, the ionic lock quickly formed in the presence of 1.00 M NaCl (see trajectories ID#4, ID#5, and ID#6 in Figure 1). It took only 3-4 ns in trajectories ID#4 and ID#5, and 14 ns in trajectory ID#6. As a comparison, in the study of Shaw et al. 30-130 ns long simulations of hβ2AR in low salt concentration were needed for the ionic lock to be formed 12. In our initial hA2AR structure, Glu228(6.30) forms a salt bridge with Arg220(6.22). It should be mentioned that the crystal structure did not contain Arg220(6.22) residue simply because its third intracellular (IC3) loop was replaced by T4 lysozyme (T4L) to make the crystallization possible. We found that Arg220(6.22) and the ionic lock forming Arg102(3.50) are competing for Glu228(6.30). At 0.15 M NaCl, this can be seen only sporadically (e.g., around 10 ns in simulation ID#1; see Figure 3), suggesting that the interaction between Glu228(6.30) and Arg220(6.22) is stable. On the other hand, this interaction frequently breaks at 1.00 M NaCl concentration, and, in parallel, the ionic lock with Arg102(3.50) is formed. In order to quantify this effect, we analyzed the distribution of Naþ-Glu228(6.30), Cl--Arg102(3.50) and Cl--Arg220(6.22) pairs (Figure 4). Interestingly, we found a clear difference depending on the NaCl concentration. Sodium ions were found in the proximity of Glu228(6.30) only in simulations at 1.00 M NaCl. Furthermore, chloride ions could be identified more frequently around Arg220(6.22) than Arg102(3.50) in simulations at 1.00 M NaCl. We suggest that both of these effects (i.e., more sodium and chloride ions in the vicinity of Glu228(6.30) and Arg220(6.22)) contribute in the weakening of the Glu228(6.30)-Arg220(6.22) interaction at high NaCl concentration.

r 2010 American Chemical Society

Figure 3. The alteration of the minimum distance between H-bond forming atoms of Arg102(3.50)-Glu228(6.30) (red) and Arg220(6.22)Glu228(6.30) (blue) in the course of the 0.15 M NaCl (ID#1, ID#2 and ID#3) and 1.00 M NaCl (ID#4, ID#5 and ID#6) simulations.

This is in agreement with the fact that Arg102(3.50) is more buried in the TM region than Arg220(6.22), which is located in the IC3 loop (see Figure S1 in the Supporting Information). Consequently, once the ionic lock with Arg102(3.50) has been formed, it can rarely be targeted by sodium and chloride ions.

1010

DOI: 10.1021/jz900403a |J. Phys. Chem. Lett. 2010, 1, 1008–1013

pubs.acs.org/JPCL

Figure 4. The distribution functions for sodium (Glu228(6.30) blue line), and chloride ions (Arg102(3.50) black line, Arg220(6.22) red line).: NOI: average number of ions in the course of the simulation.

On the other hand, the Arg220(6.22)-Glu228(6.30) interaction is more easily accessed. We have also calculated the solventaccessible surface area for Arg102(3.50) and Arg220(6.22) and found that Arg220(6.22) is significantly more accessible than Arg102(3.50) (see Figure S2 and Table S1 in the Supporting Information). Interestingly, position 6.22 is also occupied by an arginine in hβ2AR. Moreover, in the available crystal structures of the activated opsin (Protein Data Bank IDs: 3CAP and 3DQB), Glu247(6.30) forms interactions with Lys231(5.66) and Lys248(6.31) from IC3. This suggests that positively charged residues of IC3 in other class “A” GPCRs might compete for Glu(6.30) with Arg(3.50) and facilitate receptor activation. Our simulations on hA2AR suggest that this competition is regulated by the concentration of sodium chloride. It is also important to mention that binding of agonists and antagonists to several GPCRs linked to inhibition of adenylate cyclase is modulated in a reciprocal manner by monovalent cations,15-19 suggesting a similar mechanism of regulation. However, this can only be reinforced by further simulations on GPCRs. In conclusion, our MD simulation at high sodium chloride concentration resulted in an inactive state of hA2AR with an intact ionic lock and a reduced distance in the intracellular side of TM3 and TM6. We suggest that this structural model is more favorable for designing novel antagonists or inverse agonists for the hA2AR. Nevertheless, our results might also contribute to the understanding of interactions between ions and protein surfaces.19

r 2010 American Chemical Society

MATERIALS AND METHODS The hA2AR crystal structure (downloaded from the Protein Data Bank20,21) (PDB ID: 3EML), contains residues from Ile3 to Gln310, except the IC3 loop (Lys209-Ala221) and a short part of the second extracellular (EC2) loop (Gln148-Ser156). We incorporated the first two amino acids as well as the missing part of EC2 and IC3 to our model. Water molecules found to be closer than 4.5 Å to any receptor atoms in the crystal structure were also included to our system. The missing parts of the receptor were built using Molecular Operating Environment (version 2007.09).22 Optimization of the missing N-terminal EC2 and IC3 loops was performed using the simulated annealing protocol described previously.23 These parts were fully flexible during optimization, except the ω dihedral angles were restrained in the trans conformation; the remaining part of the receptor was restrained with a force constant of 10 kcal/ molÅ2. During the structure refinement, the CHARMM27 force field24 was utilized, and a 10 Å cutoff distance was used in the nonbonded interaction calculations. The IC3 loop of the lowest energy structure folded back unrealistically to the receptor surface, therefore the second lowest energy structure with an acceptable IC3 conformation was selected. The membrane environment with water molecules and ions was built with the Desmond GUI in Maestro.25 Two systems with different NaCl concentrations were built: 0.15 M (31 Naþ and 31 Cl- ions) (c1) and 1.00 M (234 Naþ and 234 Cl- ions) (c2). Ten chloride counterions were also added to both systems, 170 POPE molecules were used to build the explicit

1011

DOI: 10.1021/jz900403a |J. Phys. Chem. Lett. 2010, 1, 1008–1013

pubs.acs.org/JPCL

REFERENCES

membrane environment around the protein, and the number of water molecules was 11106 and 10960 for c1 and c2, respectively. The total number of the atoms in the two systems was 59691 (x = 78.83 Å, y = 74.71 Å, z = 100.43 Å) and 59659 (x = 76.41 Å, y = 72.41 Å, z = 106.83 Å) for c1 and c2, respectively. MD simulations were performed with the Desmond v20108 program package.26 Pre-equilibration of the two systems was performed in the following way: (1)10 kcal/ molÅ2 positional restraint was applied on hA2AR; (2) minimization with the steepest-descent and L-BFGS methods was performed (maximum step number: 2000, gradient tolerance: 1.0 kcal/molÅ; (3) gaps between the receptor surface and membrane were removed by means of 480 ps simulated annealing, where the system was heated from 250 to 310 K in the first 100 ps (Gaussian force routine was applied to prevent water influx to the gaps). After the preequilibration, the barrier and the restraints were removed, and minimization, simulated annealing (240 ps heating phase), and free MD simulation (70.08 ns) were performed. Three independent simulations were conducted for both systems using different initial velocities.27 Snapshots were saved every 4.8 ps, resulting in 14600 structures. The CHARMM27 parameters24 with the CMAP correction terms28 were applied for the proteins as well as for the membrane. Water molecules and ions were described by the TIP3P potential29 and by the Beglov parameters,30 respectively. The calculations were performed in the NPT ensemble (T = 310 K, p = 1.01325 bar) with anisotropic pressure scaling using the Berendsen thermoand barostat methods (τp = 1.0 ps, τT = 0.1 ps).31 The cutoff value of the near nonbonded interactions was set to 9 Å, and the far electrostatic interactions were calculated by the PME method32 with 80 80108 and 8072100 Fourier mesh points in the preequilibration stage and in the free MD runs, respectively. The time step was set to 2 fs, and bonded and short nonbonded interactions were calculated in every step, while long electrostatics was calculated in every third step. Molecular graphics were produced with the Visual Molecular Dynamics package.33

(1)

(2)

(3)

(4)

(5)

(6)

(7)

(8)

(9)

(10)

(11)

SUPPORTING INFORMATION AVAILABLE Details of the þ

-

solvent accessible surface area calculations and the Na -Cl radial distribution function. This material is available free of charge via the Internet at http://pubs.acs.org.

(12)

AUTHOR INFORMATION Corresponding Author: *To whom correspondence should be addressed. Address: Department of Chemical Informatics, Faculty of Education, University of Szeged, Boldogasszony sgt. 6, H-6725 Szeged, Hungary. E-mail: [email protected]; telephone: þ36-62-544-779; fax: þ3662-544-720.

(13) (14)

Present Addresses: #

(15)

Present address: D. E. Shaw Research, New York, NY 10036.

ACKNOWLEDGMENT We thank M. Labadi for excellent technical

(16)

support. We are also grateful to HPC Szeged for the computer facility.

r 2010 American Chemical Society

1012

Palczewski, K.; Kumasaka, T.; Hori, T.; Behnke, C. A.; Motoshima, H.; Fox, B. A.; Le Trong, I.; Teller, D. C.; Okada, T.; Stenkamp, R. E.; Yamamoto, M.; Miyano, M. Crystal Structure of Rhodopsin: A G Protein-Coupled Receptor. Science 2000, 289, 739–745. Cherezov, V.; Rosenbaum, D. M.; Hanson, M. A.; Rasmussen, S. G.; Thian, F. S.; Kobilka, T. S.; Choi, H. J.; Kuhn, P.; Weis, W. I.; Kobilka, B. K.; Stevens, R. C. High-Resolution Crystal Structure of an Engineered Human β2-Adrenergic G ProteinCoupled Receptor. Science 2007, 318, 1258–1265. Warne, T.; Serrano-Vega, M. J.; Baker, J. G.; Moukhametzianov, R.; Edwards, P. C.; Henderson, R.; Leslie, A. G.; Tate, C. G.; Schertler, G. F. Structure of a β1-Adrenergic G-ProteinCoupled Receptor. Nature 2008, 454, 486–491. Jaakola, V. P.; Griffith, M. T.; Hanson, M. A.; Cherezov, V.; Chien, E. Y.; Lane, J. R.; Ijzerman, A. P.; Stevens, R. C. The 2.6 Angstrom Crystal Structure of a Human A2A Adenosine Receptor Bound to an Antagonist. Science 2008, 322, 1211– 1217. Mobarec, J. C.; Sanchez, R.; Filizola, M. Modern Homology Modeling of G-Protein Coupled Receptors: Which Structural Template to Use? J. Med. Chem. 2009, 52, 5207-5216. Ballesteros, J. A; Weinstein, H. Integrated Methods for the Construction of Three-Dimensional Models and Computational Probing of Structure-Function Relations in G ProteinCoupled Receptors. Methods Neurosci. 1995, 25, 366–428. de Ligta, R. A. F.; Rivkeesb, S. A.; Lorenzenc, A.; Leursd, R.; IJzermana, A. P. A “Locked-On,” Constitutively Active Mutant of the Adenosine A1 Receptor. Eur. J. Pharmacol. 2005, 510, 1–8. Gao, Z. G.; Kim, S. K.; Gross, A. S.; Chen, A.; Blaustein, J. B.; Jacobson, K. A. Identification of Essential Residues Involved in the Allosteric Modulation of the Human A(3) Adenosine Receptor. Mol. Pharmacol. 2003, 63, 1021–1031. Farrens, D. L.; Altenbach, C.; Yang, K.; Hubbell, W. L.; Khorana, H. G. Requirement of Rigid-Body Motion of Transmembrane Helices for Light Activation of Rhodopsin. Science 1996, 274, 768–770. Park, J. H.; Scheerer, P.; Hofmann, K. P.; Choe, H. W.; Ernst, O. P. Crystal Structure of the Ligand-Free G-Protein-Coupled Receptor Opsin. Nature 2008, 454, 183–187. Scheerer, P.; Park, J. H.; Hildebrand, P. W.; Kim, Y. J.; Krauss, N.; Choe, H. W.; Hofmann, K. P.; Ernst, O. P. Crystal Structure of the Active G-Protein-Coupled Receptor Opsin in Complex with a C-Terminal Peptide Derived from the G-R Subunit of Transducin. Nature 2008, 455, 497–502. Dror, R. O.; Arlow, D. H.; Borhani, D. W.; Jensen, M. Ø.; Piana, S.; Shaw, D. E. Identification of Two Distinct Inactive Conformations of the β2-Adrenergic Receptor Reconciles Structural and Biochemical Observations. Proc. Natl. Acad. Sci. U.S.A. 2009, 106, 4689–4694. Frishman, D.; Argos, P. Knowledge-Based Secondary Structure Assignment. Proteins 1995, 23, 566–579. De Lean, A.; Stadel, J. M.; Lefkowitz, R. J. A Ternary Complex Model Explains the Agonist-Specific Binding Properties of the Adenylate Cyclase-Coupled β-Adrenergic Receptor. J. Biol. Chem. 1980, 255, 7108–7117. Limbird, L. E. Receptors Linked to Inhibition of Adenylate Cyclase: Additional Signaling Mechanisms. FASEB J. 1988, 2, 2686–2695. Pert, C. B.; Pasternak, G. W.; Snyder, S. H. Opiate Agonists and Antagonists Discriminated by Receptor Binding in Brain. Science 1973, 182, 1359–1361.

DOI: 10.1021/jz900403a |J. Phys. Chem. Lett. 2010, 1, 1008–1013

pubs.acs.org/JPCL

(17)

(18)

(19)

(20)

(21) (22) (23)

(24)

(25) (26)

(27)

(28)

(29)

(30)

(31)

(32)

(33)

Urwyler, S. Mono- and Divalent Cations Modulate the Affinities of Brain D1 and D2 Receptors for Dopamine by a Mechanism Independent of Receptor Coupling to Guanyl Nucleotide Binding Proteins. Naunyn-Schmiedeberg's Arch. Pharmacol. 1989, 339, 374–382. Nunnari, J. M.; Repaske, M. G.; Brandon, S.; Cragoe, E. J., Jr.; Limbird, L. E. Regulation of Porcine Brain R2-Adrenergic Receptors by Naþ, Hþ and Inhibitors of Naþ/Hþ Exchange. J. Biol. Chem. 1987, 262, 12387–12392. Jungwirth, P.; Winter, B. Ions at Aqueous Interfaces: FromWater Surface to Hydrated Proteins. Annu. Rev. Phys. Chem. 2008, 59, 343–66. Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.; Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank. Nuc. Acid. Res. 2000, 28, 235–242. Berman, H. M.; Henrick, K.; Nakamura, H. Announcing the Worldwide Protein Data Bank. Nat. Struct. Biol. 2003, 10, 980. Molecular Operating Environment (MOE 2007.09); C.C.G., Inc: Montreal, Quebec, Canada.  . Comparative Study of Eight Oxytocin J ojart, B.; Marki, A Antagonists by Simulated Annealing. J. Mol. Model. 2006, 12, 823–828. MacKerell, A. D.Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; et al. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins. J. Phys. Chem. B 1998, 102, 3586–3616. http://www.deshawresearch.com/resources_desmond.html (accessed September 1, 2008). Bowers, K. J.; Chow, E.; Xu, H.; Dror, R. O.; Eastwood, M. P.; Gregersen, B. A.; Klepeis, J. L.; Kolossv ary, I.; Moraes, M. A.; Sacerdoti, F. D.; Salmon, J. K.; Shan, Y.; Shaw, D. E. Scalable Algorithms for Molecular Dynamics Simulations on Commodity Clusters. Proceedings of the ACM/IEEE Conference on Supercomputing (SC06), Tampa, Florida, November 11-17, 2006. Caves, L. S. D.; Evanseck, J. D.; Karplus, M. Locally Accessible Conformations of Proteins: Multiple Molecular Dynamics Simulations of Crambin. Protein Sci. 1998, 7, 649–666. MacKerell, A. D.Jr.; Feig, M.; Brooks, C. L.III. Extending the Treatment of Backbone Energetics in Protein Force Fields: Limitations of Gas-Phase Quantum Mechanics in Reproducing Protein Conformational Distributions in Molecular Dynamics Simulations. J. Comput. Chem. 2004, 25, 1400–1415. Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926– 935. Beglov, D.; Roux, B. Finite Representation of an Infinite Bulk System: Solvent Boundary Potential for Computer Simulations. J. Chem. Phys. 1994, 100, 9050–9063. Berendsen, H. J. C.; Postma, J. P. M.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. Essman, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593. Humphrey, W.; Dalke, A.; Schulten, K. VMD - Visual Molecular Dynamics. J. Mol. Graphics 1996, 14, 33–38.

r 2010 American Chemical Society

1013

DOI: 10.1021/jz900403a |J. Phys. Chem. Lett. 2010, 1, 1008–1013