Induced Effects of Sodium Ions on Dopaminergic G-Protein Coupled Receptors

G-protein coupled receptors, the largest family of proteins in the human genome, are involved in many complex signal transduction pathways, typically activated by orthosteric ligand binding and subject to allosteric modulation. Dopaminergic receptors, belonging to the class A family of G-protein coupled receptors, are known to be modulated by sodium ions from an allosteric binding site, although the details of sodium effects on the receptor have not yet been described. In an effort to understand these effects, we performed microsecond scale all-atom molecular dynamics simulations on the dopaminergic D2 receptor, finding that sodium ions enter the receptor from the extracellular side and bind at a deep allosteric site (Asp2.50). Remarkably, the presence of a sodium ion at this allosteric site induces a conformational change of the rotamer toggle switch Trp6.48 which locks in a conformation identical to the one found in the partially inactive state of the crystallized human β2 adrenergic receptor. This study provides detailed quantitative information about binding of sodium ions in the D2 receptor and reports a possibly important sodium-induced conformational change for modulation of D2 receptor function.

D 2 receptor (accession P14416) was retrieved from the Swiss-Prot database [1]. The generalized numbering scheme proposed by Ballesteros and Weinstein [2] was used for residues belonging to helix regions of the G-protein coupled receptors (GPCR). The high flexible N-terminal (sequence 1-29) and intracellular loop 3 (sequence 221-365) were removed since no adequate template exists and due to the fact that the putative binding site of sodium is located at Asp2.50 in the seven transmembrane (TM) region. The shortened sequence was aligned with the ClustalX software [3,4] using the PAM250 matrix and "gap open" and "gap elongation" penalties of 10 and 0.05, respectively. The resulting multiple sequence alignment was realigned with the crystal structure of the human ß 2 adrenergic G-protein coupled receptor (PDB entry 2RH1) [5,6] introducing secondary structure information derived from the crystal structure to avoid gaps within the seven helical segments. The alignment was then manually refined to ensure a perfect alignment of the highly conserved residues of GPCR superfamily, according to Baldwin et al. [7]. Extension of each helix was contemplated by taking into account the experimental length of the 2RH1 helices and the sequence conservation. The conserved disulfide bond between residues C3.25 at the beginning of TM3 and the cysteine in the middle of the extracellular loop 2 (a feature highly common among GPCRs) was also created and kept as a constraint in the geometric optimization. 3D models were then built using the MODELLER suite of programs [8] which yielded 15 candidate models. From these candidates, the best structures according to the MODELLER objective function and to visual inspection were selected. Models with interruptions or gaps in the TM regions, as identified by visual inspection, were discarded. The resultant structures of the receptors were optimized using the Amber99 force field [9] applying the molecular modeling program MOE (Molecular Operating Environment; Chemical Computing Group). PROCHECK software [10] was used to assess the stereochemical quality of the minimized structures resulting in good quality parameters with an excellent distribution of Psi and Phi angles in the Ramachandran plot (over 90% of the residues are in the most favored regions).

Generation of the dopaminergic D 2 receptor/POPC membrane system
In order to place the D 2 receptor model into the bilayer membrane, a hole was generated in a preequilibrated palmitoyloleoylphosphatidylcholine (POPC) bilayer membrane (membrane builder tool of the VMD version1.8.6) by removing POPC molecules. Lipids which were in close contact with the protein atoms (<1 Å distance from any protein atoms) were deleted, and water molecules were added to the system. An all-atom model of the D 2 receptor was generated using the CHARMM27 force field and parameter set with the TIP3 water model. The protonation state of titratable groups was predicted for a pH value at 7.4 based on PROPKA [11] using the implemented prediction tool of the MOE package. Seven aspartates located in the extracellular surface and TM region were built in the unprotonated form. The putative sodium binding site Asp2.50 in the TM domain were modeled with unprotonated carboxylic acid side chains (negatively charged) as the experimentally suggested presence of a sodium ion in the receptor would compensate the negative charge of Asp2.50. Furthermore, six glutamate residues which are located in extra-and intracellular regions were modeled in an unprotonated form; five histidine residues were kept in their neutral form; all arginine and lysine residues were modeled in protonated form. The protein model had a total net charge of +8e, which is balanced with 27 sodium and 35 chloride ions in the aqueous layer of 11094 water molecules containing a phospholipids bilayer with 170 POPC (palmitoyloleoylphosphatidylcholine) molecules. The coordinates for water and ions were generated using the solvate and autoionize modules of VMD1.8.6.

Equilibration protocol of the D 2 receptor/POPC membrane system
The initial equilibration of the system was carried out with NAMD2.6, using rigid bonds, time step 2 fs, cutoff 9 Å, PME grid size of 80x80x96, Langevin damping 0.1/ps. The system was first minimized for 1000 steps, then thermalized for 80 ps at 300 K reinitializing the velocities every 20 ps. A positional restrain of 1 kcal/mol/Å 2 was applied to all heavy atoms of the receptor. The system was then equilibrated at 1 atm for 100 ps, while the restraints were progressively reduced to zero. A final run of 20 ns was performed to equilibrate the pressure during which the potential energy and volume of the simulation box S4 were monitored to assure equilibration. After 15 ns both quantities remain stable. The final configuration box (78.76,78.76,94.21) Å was further equilibrated at constant temperature (NVT) using ACEMD for further 20 ns at 300 K, with a time step equal 2 fs, cutoff 9 Å, switch distance 7.5 Å, rigid bonds and a PME grid size of 80x80x96.

Production runs
For the production runs we have used ACEMD running on local workstations equipped with 4 Nvidia GPUs each and the distributed computing project GPUGRID. ACEMD is a new generation molecular dynamics software which runs exclusively on graphics processing units (GPUs) at the equivalent speed of hundred of standard processors [12,13] (see also reference [14] for a review of the impact of accelerator processors and graphics cards on computational biology). GPUGRID is a distributed computing project based on the use of graphics processors around the globe [15] and the ACEMD molecular dynamics software [12]. At the time of these runs GPUGRID comprised of over 4000 GPUs for a total peak performance of 1 petaflop, and a sustained performance of 150 Teraflops.
After the initial equilibration phase and for all the production runs the simulation are run in the NVT ensemble, with the same parameters used for the equilibration but a timestep of 4 fs thanks to the use of the hydrogen mass repartition scheme implemented in ACEMD [16]. This elegant scheme allows for longer timesteps up to 4 fs by using the property that individual atom masses do not appear explicitly in the equilibrium distribution, therefore changing them only affects the dynamical properties of the system (marginally) [16] but not the equilibrium distribution. The lipid cross section is at equilibrium 64.6 Å 2 /lipid as estimated from the equilibrated configuration and in good agreement with the experimental cross section for POPC [17]. As we run the long simulations in the NVT ensemble this cross section is maintained even during long runs and large sampling.
We have reported in this work three final runs according to Table 1 of the main manuscript: MD1, MD2, MD3.

Single, long trajectory (MD1).
A single run of 1.097 µs was performed on a workstation using ACEMD in parallel over three Nvidia 275 GPUs for a total time of 66 days of continuous computation. A movie of the entire simulation is available online at http://www.vimeo.com/7900216 and in Video S1.
Volumetric Maps (MD2). The second production run (MD2 in Table 1) consists of 100 independent runs of approximately 50 ns each for a total sampling of 4.7 µs which were successfully returned by the GPUGRID server. The trajectory dumping frequency was set to 50 ps to reduce data upload to the server.
The entire calculations required a couple of months although occupying only a fraction of the capability of GPUGRID.
For the analysis we have first discarded the first 20 ns of each run to avoid any bias of the initial condition, sampling the trajectory every 100 ps. The concatenated trajectory data was then wrapped for the periodic boundary conditions centered around the protein and then aligned to the same initial configuration. The tool VolMap from VMD was used to produce the volumetric map for ions and water (atom radius of 1.0 Å) using a grid size resolution of 0.5 Å and computing the ion and water occupancy per grid cell. The volumetric map of sodium ions (yellow isosurface, Figure 1A) was computed for a chemical potential µ = -4k B T which corresponds to an isosurface c = 0.002465 particle/ Å 3 based on: Δμ = -k B T log(c/c bulk ), with the bulk ions concentration being 150 mM NaCl (c bulk = 0.000090307 particle/ Å 3 ). The volumetric map of water (blue isosurface, Figure 1A) was computed for the chemical potential µ = 0.514k B T which corresponds to an isosurface c = 0.20000 particle/ Å 3 relative to the concentration of water at 300 K, 1 ATM (c water = 0.334438 particle/ Å 3 ).

Metadynamics free energy calculations (MD3). Metadynamics [18] is a biased dynamics technique
widely used to improve sampling for free energy calculations over a set of multidimensional reaction coordinates which would not be sampled exhaustively with normal unbiased simulations. It is implemented in the molecular dynamics code ACEMD using the PLUMED plugin interface [19].
Metadynamics is used to complement the unbiased simulation. According to the unbiased simulations, a S6 two-dimensional reaction coordinate is chosen to represent the slow degrees of freedom: The first coordinate is the z position of the bound sodium ion, while the second one is the dihedral angle χ 2 of tryptophan W6.48 formed between atoms (CA-CB-CG-CD1). The metadynamics parameters are set to a Gaussian hill height of 0.2 kcal/mol with a spread of 0.2 Å for the z coordinate and 0.1 rad for the dihedral angle. The deposition rate is one hill every 2 ps and a well-tempered bias factor of 8. Small changes on these reference parameters did not seem to affect the results. In order to estimate not only the free energy, but also the error associated with it, we performed 25 independent metadynamics runs of 14 ns each on GPUGRID. Each run started from an equilibrated configuration with the sodium ion outside of the deep binding pocket of the receptor. The simulations were binned over 100 bins between -20 to -1 Å for the z coordinate and between -π and π for the dihedral angle χ 2 . The average free energy W(z, χ 2 ) is reconstructed as the average over the estimate of the number states N(z, χ 2 ) = <exp(-W i (z, χ 2 )/k B T> i of each individual metadynamics run, where W i (z, χ 2 ) is the free energy determined by the individual metadynamics run i, and The error in the estimation of W(z, χ 2 ) is measured similarly by computing the standard deviation in the number of states S(z, χ 2 ) = <(exp(-Wi(z, χ 2 )/k B T-N) 2 > 1/2 . The total free energy error is then estimated as W err (z, χ 2 ) = -k B T log (N(z, χ 2 ) -S(z, χ 2 )) + k B T log (N(z, χ 2 )+S(z, χ 2 )), as the free energy difference between an error of a standard deviation in N.

System stability
The generated system is stable during 1.1 μs simulation (MD 1). First, no major structural changes in the TM region of the D 2 receptor occurred. As expected, larger structural rearrangements in the more flexible loop regions were observed (Figure 1) Second, comparison of the TM topology of the final D 2 receptor model relative to experimental X-ray crystal structures demonstrates that the primary 7TM topology is well maintained during simulation (RMSD < 2 Å in TM region, Table 1). Third, the binding pocket of the S7 D 2 receptor is filled with water molecules which mostly coincide with water molecules found in the X-ray crystal structure of the β 2 adrenergic receptor (see Figure 2).   Figure 2. Volumetric map of bulk water (blue isosurface) computed over 4.7 μs simulation data (MD2) in superimposition with the X-ray crystal water of the β 2 adrenergic receptor (PDB ID 2RH1).

Receptor states during the activation process
In the process of GPCR activation, the rotamer switch Trp6.48 undergoes a major conformational change [20] and rotates from a vertical position in which the nitrogen is directed to TM2 (gauche (g) conformation) into a horizontal position in which the nitrogen is directed to TM5 (trans (t) conformation) [21].During this process the receptor has to undergo intermediate conformational states