Subscriber access provided by United Arab Emirates University | Libraries Deanship
B: Biophysics; Physical Chemistry of Biological Systems and Biomolecules 4
Charge Discrimination in P2X Receptors Occurs in Two Consecutive Stages Gustavo Pierdominici-Sottile, Vanesa Racigh, Agustín Ormazábal, and Juliana Palma J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.8b10249 • Publication Date (Web): 09 Jan 2019 Downloaded from http://pubs.acs.org on January 9, 2019
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Charge Discrimination in P2X4 Receptors Occurs in Two Consecutive Stages Gustavo Pierdominici-Sottile, Vanesa Racigh, Agust´ın Ormaz´abal, and Juliana Palma∗ Departamento de Ciencia y Tecnolog´ıa, Universidad Nacional de Quilmes, S´aenz Pe˜ na 352, Bernal B1876BXD, Argentina. Consejo Nacional de Investigaciones Cient´ıficas y T´ecnicas. E-mail:
[email protected] Phone: +54 (11) 4365 7100 ext 5657
1
ACS Paragon Plus Environment
The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Abstract P2X receptors are a group of trimeric cationic channels that activate by ATP. They perform critical roles in the membranes of mammalian cells and their improper functioning is associated with numerous diseases. Despite the vast amount of research devoted to them, several aspects of their operation are currently unclear. Among these remain the causes of their charge selectivity. We present the results of molecular dynamics simulations which shed light on this issue for the case of P2X4 channels. We examined in detail the behavior of Na+ and Cl− ions inside the receptor. The examination reveals that charge discrimination occurs in two stages. First, cations bear precedence over anions to enter the extracellular vestibule. Next, cations at the extracellular vestibule are more likely to cross the pore than anions in an equivalent position. In this manner, a thorough but straightforward analysis of computational simulations suggests a stepwise mechanism, without a unique determinant factor.
Introduction P2X receptors are a family of cationic channels triggered by the binding of ATP. 1–3 They are widely distributed among the tissues of mammals and display a vast range of physiological roles. 4–15 It has been found that their malfunction is associated with neurological disorders, cardiovascular problems, and cancer. 10,16–18 Because of this, they have become promising targets for the development of new drugs. 19 The number of investigations seeking to clarify different aspects of their functioning is huge. 20–27 Several reviews describe and put them into context. 19,28–30 In 2009 the X-ray structure of zebra fish P2X4 (zfP2X4 ) in the closed conformation was elucidated. 31 This represented a significant step forward in our understanding of these channels. 32–34 A few years later the crystal structure of zfP2X4 in the ATP-bound open form was also resolved. 35 The comparison between the two structures revealed the conformational changes that take place in the extracellular and the transmembrane (TM) domains, upon 2
ACS Paragon Plus Environment
Page 2 of 32
Page 3 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
the binding of ATP. However, the behavior of the cytoplasmic side remained elusive since this part was lacking in the reported structures. Only recently has this gap of information been bridged over, after the X-ray structures of human P2X3 (hP2X3 ) in the closed, open, and desensitized conformations were resolved. 36 This open structure showed, for the first time, an important fraction of the intracellular side of the channel. Most P2X receptors are homotrimeric and possess threefold symmetry around the transmembrane axis. The chains have a dolphin-like shape. Their N and C terminals lie in the cytoplasmic side. These terminals are linked to the extracellular domain via two transmembrane alpha-helix motifs, TM1 and TM2. The extracellular region is relatively large. The three ATP binding sites are located at the outside of its upper part. Opposite to them, inside the protein, there is a small chamber known as the upper vestibule. Below, there is a sizable cavity formed by two chambers: the central vestibule and the extracellular vestibule. The extracellular vestibule lies above the cell membrane and is flanked by three wide lateral portals, typically referred to as fenestrations. These fenestrations are the primary path of access to the channel for ions and molecules. A scheme of the protein, portraying all these features, is presented in Fig. 1. There are seven subtypes of P2X receptors. They are identified with numbers from 1 to 7 (P2X1 , P2X2 , ... , P2X7 ). 37 Subtypes 1, 2, 4 and 7 exhibit high selectivity for the permeation of monovalent cations over anions (PCl /PNa ≤ 0.1) 38–40 while subtype 5 presents the lowest selectivity (PCl /PNa ∼ 0.5). 39,41 The reasons for this behaviour are still unknown. The crystal structure of zfP2X4 shows that there is a net negative charge inside the extracellular vestibule, just above the entrance to the pore, due to the presence of residues D59 and D61. It has been proposed that they favor the preferential entrance of Na+ and Ca+2 over Cl− . 31,42 However, further experiments have proved that the presence of these charges does not prevent the entrance of anions. Scanning Cysteine Accessibility Method assays have demonstrated that negatively charged molecules can reach, not only the locations of D59 and D61 but positions which are ∼13-15 ˚ A inwards the pore. 43,44 One explanation for this
3
ACS Paragon Plus Environment
The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
might be that neither the lateral portals nor the extracellular vestibule perform a critical role in charge selection. Alternatively, discrimination could take place in an inner part of the channel. Indeed, it has been found that mutations lying deeper into the pore modify the PCl /PNa ratio in P2X2 and P2X7 . 40,45 However, the residues spotted in this way are not highly conserved among the members of the P2X family. 29 Therefore, it is doubtful that they could fulfill a key role in charge selection. 29 Finally, one might also think that there is not a unique factor that controls selectivity, but rather a number of contributing elements. In such case, the observed permeability would result from the interplay between all these elements, turning the situation more complex to analyze. The computational studies on P2X receptors are scarcer than the experimental ones. Normal mode analysis (NMA) and molecular dynamics (MD) simulations were employed to analyze the alternative poses that ATP and similar molecules can adopt when they are in the binding pocket. 46,47 In addition, an NMA study described the coupling between the binding of ATP and the opening of the pore. 32 Electrostatic free energy computations were performed to examine alternative pathways for the access of ions. 42 The presentation of the crystal structure of hP2X3 was accompanied by MD simulations which showed the cytoplasmic fenestrations were a plausible route for ion egress. 36 Finally, the dynamics of the closed conformation was analyzed by decomposing its fluctuations into intra-chain deformations and inter-chain movements. 48 In this article, we present the results of a computational study focused on the behavior of Na+ and Cl− ions in interaction with the P2X4 receptor. We conducted MD simulations of the zfP2X4 channel, embedded in a POPC membrane and solvated in 0.15M NaCl solution. The dynamics of the ions was analyzed and their most important interactions were identified. The closed and open conformations of the receptor were studied in this way. The analysis of the trajectories revealed significant differences between the behavior of Na+ and Cl− . Umbrella sampling (US) simulations were carried out to force the ions to move across the channel. In this way, we could estimate the corresponding potentials of mean force. On the
4
ACS Paragon Plus Environment
Page 4 of 32
Page 5 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
whole, the analysis of the results allows rationalizing previous experimental findings, sheds light onto the behavior of cations and anions inside P2X4 receptors and provides relevant information to understand their preference for the passage of cations.
Theoretical methods The selectivity of ionic channels can be studied by determining the potential of mean force (PMF) for the path followed by the ions that go from the external solution to the internal one passing through the pore. This strategy, however, is difficult to implement in the current case. Previous theoretical and experimental studies have determined that ions enter the extracellular vestibule passing through the lateral fenestrations. 42,43 Then, they move across the pore to reach the cytoplasmic side. Hence, the process is composed of ionic movements in various directions. Because of this characteristic, we have been unable to find a single reaction coordinate capable of describing the entire process. Fortunately, we identified a way to circumvent the problem. We noticed that the distribution of the ions between the external solution and the extracellular vestibule got equilibrated within an affordable computational time. According to this, we performed standard simulations to let the ions to equilibrate between the external solution and the vestibule. From these calculations, we estimated the probability of having Na+ or Cl− within the receptor, in an appropriate position to initiate the movement towards the intracellular solution. In addition, we used umbrella sampling (US) simulations to move the ions across the pore, starting from the extracellular vestibule. Then we computed the corresponding PMFs. This allowed us to assess the probability that the ions pass through the pore. Finally, gathering both pieces of information, we were able to estimate the probability that the ions travel from the external solution to the internal one, passing through the pore. In the next sections we present the methods employed to run the simulations and analyze the results. In particular, we describe the protocol followed to prepare the models of zfP2X4
5
ACS Paragon Plus Environment
The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
embedded in a POPC bilayer. Also, we provide the numerical details needed to reproduce the calculations. Finally, we present some equations required for the analysis of the results, and discuss the computations performed to estimate their statistical uncertainty and consistency.
Models construction The initial structure of the closed conformation was taken from the Protein Data Bank, entry 3I5D. 31 The protein model employed in the current calculations was built from this crystal structure, by applying the adjustments described in Ref. 48. The dimensions of the lipid bilayer was also the same as in Ref. 48. The preparation of the model for the open conformation was more complicated. It required two crystal structures. One of them was the open form of zfP2X4 , PDB entry 4DW1. 35 The other one was the open form of hP2X3 , PDB entry 5SVK. 36 The latter is the only crystal structure of an open P2X receptor which contains a significant fraction of the cytoplasmic domain. Thus, we built a more complete model for open zfP2X4 using the cytoplasmic part of hP2X3 as a template. The next paragraphs describe in detail how we assembled the model. The numbering of the residues and the names of the atoms used hereafter are those found in the 4DW1 and 5SV3 structures. The sequences of zfP2X4 and hP2X3 were aligned with program Clustal Omega 49 as provided by the Uniprot Web Page, http://www.uniprot.org/align/. The results of the alignment are presented in the Supporting Information Section. To model the N-terminal region and most of the TM1 domain of zfP2X4 , we superposed the Cα atoms of L32, L33 and I34, belonging to TM1 of hP2X3 , with the corresponding residues of zfP2X4 : A40, L41 and V42. After that, the N-terminal part of zfP2X4 up to V42 was deleted and the backbone atoms of the corresponding hP2X3 residues were added in their place. In this situation, the first residue of the N terminal region corresponds to F15 of zfP2X4 (F7 of hP2X3 ). Similarly, to model the C-terminal region and most of the TM2 domain of zfP2X4 , the Cα atoms of V326, A327 and A328, belonging to TM2 of hP2X3 , were superposed with G342, A343 6
ACS Paragon Plus Environment
Page 6 of 32
Page 7 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
and G344, of zfP2X4 . Then, the C-terminal region of zfP2X4 up to G342 was deleted and the backbone atoms of the corresponding hP2X3 residues were added in their place. This generates the conformation of TM2 up to D380 of zfP2X4 (E363 of hP2X3 ). After that, in those locations in which the two sequences differ, we mutated the residues of hP2X3 to recover the native sequence of zfP2X4 . Besides, positions N78, N187 and H252 that had been modified to obtain the 4DW1 structure, 35 were mutated back to the native sequence. We added residues D9 to C14 to the N terminal part and residues F381 to K389 to the C-terminal part, using program I-TASSER. 50 The Leap module of the AMBER16 package 51 concluded the preparation of the model by incorporating all the missing atoms of the edited PDB file. Therefore, while the crystal structure of open zfP2X4 consists in 324 residues, from R36 to I359, the model employed in the current calculations consisted of 381 residues, from D9 to K389. The 4DW1 structure presents intersubunit crevices that should not be present in membrane-embedded receptors. 52 Previous MD simulations have demonstrated that lipid molecules can enter the channel through these crevices. When building our model for P2X4 , we purposely chose to superpose the structure of P2X4 with that of P2X3 using residues located in the upper part of the transmembrane domains. By this action, the TM domains of P2X4 got mostly aligned with those of P2X3 . This drastically reduced the size of the intersubunit crevices. Nonetheless, in all the MD simulations of the open structure we checked that there were no lipids inside the channel. The protein model was soaked into a square lipid bilayer of 138.0 ˚ A2 with 256 POPC molecules in each layer. The graphical user interface of the CHARMM membrane builder 53,54 was used for this purpose. The protein was displaced along the Z axis so that the top of the membrane was aligned with residue Y52.
Standard MD simulations Standard MD simulations were performed for both, the closed and open forms of the receptor. The following simulation protocol was applied for the two models. 7
ACS Paragon Plus Environment
The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Initially, the system was fed into the Leap module of AMBER16 where it was first neutralized. Then Na+ and Cl− ions were incorporated to produce a 0.15M concentration. Water molecules were added to form an octahedral cell. Special care was devoted to introduce waters molecules into the extracellular vestibule and transmembrane channel. Otherwise either the simulations become unstable or the pore shrinks. The Amber99SB force field 55 was employed for the protein and water molecules while the Lipid14 force field was chosen for POPC. 56 The SANDER and PMEMD modules of the AMBER16 package were employed to run the simulations. 51 This initial structure was first minimized. Next, it was heated from 0 K to 100 K at NVT conditions during 120 ps. The aim of this stage is to achieve a temperature high enough to assure the proper functioning of the barostat. At that point, we changed from NVT to NPT conditions to allow the density to relax. We continued the heating for 1 ns until the system got to 303 K. An harmonic restraint of 1.5 kcal/mol˚ A2 was applied on the Cα atoms of the protein and the oxygen atoms of the water molecules in the two heating periods. This was followed by four consecutive 5.0 ns simulations at 303 K in which the restraints were ˚2 ). gradually diminished (0.5, 0.1, 0.05 and 0.01 kcal/mol A A final equilibration stage of 2.5 ns was performed for the open structure while 7.0 ns were required for the closed structure. These times were needed to equilibrate the spacial distribution of the ions. Initially, the models had no ions inside the receptor. The simulations showed this situation is unrealistic since some Na+ moved into the central and extracellular vestibules even in the heating stages. We considered that the spacial distribution of the ions was equilibrated when the number the ions inside the receptor fluctuated about a definite average value. All the results presented below were obtained using snapshots taken from the equilibrated system. The production runs were performed with a restrain of 0.01 kcal/mol˚ A2 on the Cα atoms of residues I100 to I260. These residues are located in the upper vestibule, far apart from the entrance to the pore. During a few long simulations with no restrictions, we noted that charged residues lying in the lower part of the extracellular vestibule got close
8
ACS Paragon Plus Environment
Page 8 of 32
Page 9 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
to the charged heads of the lipids. When this occurred, the channel tilted as a consequence of the strong Coulombic interaction, and a small part of the fenestration got inside the membrane. We considered this behaviour an artifact of the model. The tiny restrains on the Cα atoms of residues I100 to I260 prevented this from occurring. The SHAKE algorithm was implemented to constrain the bonds involving hydrogen atoms. Accordingly, the time step could be set to 2.0 fs. A cutoff radius of 10.0 ˚ A was imposed to the electrostatic interactions and the Particle Mesh Ewald method was employed to calculate the long-range Coulombic force. 57,58 Standard MD simulations were started from the final structure of the equilibration stage. For the closed conformation we run eight independent simulations of 50 ns (adding up a total of 0.4 µs) collecting snapshots every 50 ps. For the open conformation we run twelve simulations of 50 ns (a total sampling time of 0.6 µs), collecting snapshots every 250 ps. In all cases the initial velocities were randomly chosen from a Maxwellian distribution at 303 K. The standard MD simulations were employed to visualize the behavior of the ions and also to perform a statistical analysis. We estimated the probability of finding Na+ or Cl− inside the protein, at different locations along the pore. We considered an ion was inside the protein when it was within a cylinder of radius rc = 16.0 ˚ A, whose axis was aligned to the pore axis. From now on we will refer to this axis as the Z-axis. The left pannel of Fig. 1 shows a picture of this cylinder. The probability density to find Na+ or Cl− inside the protein, at different locations along the Z-axis, were calculated as follows. Let’s Nj (Z ∗ , Z ∗ + ∆Z; r < rc ) be the number of snapshots in which ion j is found inside the cylinder of radius rc (see Fig. 1), having its Z-coordinate between Z ∗ and Z ∗ + ∆Z. Then, the probability that ion j fulfills the two conditions is, Pj (Z ∗ , Z ∗ + ∆Z; r < rc ) =
Nj (Z ∗ , Z ∗ + ∆Z; r < rc ) , Nsnap
where Nsnap is the total number of snapshots. In order to consider the contributions from all the ions of a given kind, we sum the probability of each ion and divide by the number of 9
ACS Paragon Plus Environment
The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 32
Figure 1: Scheme of the cylinders used to estimate the probability densities. Left panel: cylinder used to compute the probability distributions of the ions along the pore axis. Right panel: scheme of the cylindrical slab employed to compute the radial probability distributions. ions of that kind, Nions ,
∗
∗
P (Z , Z + ∆Z; r < rc ) =
N ions X j=1
=
Pj (Z ∗ , Z ∗ + ∆Z; r < rc ) , Nions
N (Z ∗ , Z ∗ + ∆Z; r < rc ) , Nions .Nsnap
where N (Z ∗ , Z ∗ + ∆Z; r < rc ) is the total number of snapshots in which any ion of the kind under consideration is inside the cylinder and has its Z coordinate in the interval Z ∗ : Z ∗ + ∆Z. If this interval is small, the probabilities so computed depend linearly on the size of ∆Z. Therefore, it is more convenient to evaluate probability densities, ρ(Z ∗ , r < rc ), from ρ(Z ∗ , r < rc ) =
P (Z ∗ , Z ∗ + ∆Z; r < rc ) ∆Z
(1)
since these quantities are independent of the interval size. It has to be noted that the probability densities calculated with Eq. 1 do not integrate to 1.0 when the whole range of the Z coordinate is considered. They only do that if rc is set to a very large value so that 10
ACS Paragon Plus Environment
Page 11 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
the whole volume of the simulation box is taken into account.
Umbrella sampling simulations Umbrella Sampling computations were started from snapshots corresponding to the equilibrated open structure, with Na+ ions at the top of the central vestibule. These Na+ were not restricted at all. In spite of this, they never moved to the extracellular vestibule. Therefore, they never directly interacted with the ion subjected to the bias potential. Instead, they hovered around their initial location which is ∼25-35 ˚ A away from the entrance of the pore. The reaction coordinate was defined as the value of the Z-coordinate. The origin was set at the center of mass of the Cα atoms of D59 (see Fig. 1). These residues are sited at the extracellular vestibule, above the entrance to the pore. The Z-axis was pointing downwards. Accordingly, positive Z-values correspond to the cytoplasmic and membrane regions while negative values correspond to the extracellular side. The reaction coordinate was sampled from 0.0 ˚ A to -8.0 ˚ A and from 0.0 ˚ A to 40.0 ˚ A, with a spacing of 0.7 ˚ A between the centers of adjacent windows. The last conformation of each window was used as the starting point for the next one. A harmonic restraint of 7.5 kcal/mol˚ A2 was applied to force the ion to wander around the center of each window. Individual US-simulations lasted for 10.0 ns. The actual values of the reaction coordinate were recorded every 4.0 fs. The first 2.0 ns of each simulation were not considered to compute the PMF. A cylindrical restriction was applied to the ion subjected to the bias potential in US-windows located at the extracellular vestibule and cytoplasmic domains. This was done to facilitate the convergence of the sampling. 59 The radius of the restraining potential was set to 30.0 ˚ A in the extracellular side and 40.0 ˚ A in the cytoplasmic side. This large value was required to include the whole cytoplasmic cap inside the cylinder. We checked that variations of both radii of ±5.0 ˚ A had no impact on the results. Other ions of the solution were not allowed to move into the cylinders to get a one-ion PMF. 60 We implemented two alternative procedures to calculate the PMFs from the data collected 11
ACS Paragon Plus Environment
The Journal of Physical Chemistry 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
in the US simulations. They are the Weighted Histogram Analysis Method (WHAM) 61 and the Dynamic Histogram Analysis Method (DHAM). 62 For WHAM we utilized the program developed by Grossfield’s and co-workers. 63 For DHAM we developed our own FORTRAN code. It employs subroutine DGEEV of the LAPACK library to diagonalize the matrix of transition probabilities. In both cases, the number of bins employed in the calculation was set as four times the number of US-windows. The two procedures afforded the same PMFs within the statistical uncertainty. This agreement indicates that the bias potential applied in the US-simulations was appropriate. To evaluate the convergence of the results we divided the data of each US-window into four equivalent sets and we computed a PMF from the data of each set. Figure S2 of the supporting information presents the profiles so obtained. The similarity between the alternative profiles indicates an adequate convergence of the results. The PMF presented below was calculated from the complete set of data. Its statistical uncertainty was estimated by the standard deviation of the curves computed with the four sets. Finally, to assess the consistency between the raw data provided by the US-simulations and the probabilities calculated from WHAM and DHAM, we computed the relative entropy at each window, using Eq. 2 of Ref. 64. Entropies below 0.0025 were obtained for all windows, proving the consistency of the results (see Fig. S3 of the Supporting information section).
Results Ions in the extracellular domain In this section we present and discuss the results afforded by the standard MD simulations. We observed that, when there are no other ions inside, several Na+ enter the receptor and move towards the upper part of the central vestibule. They attain a stable condition there, remaining for at least ∼25 ns. Some ions stay there for the complete length of the simulation. This process occurs for both, the closed and open structures, and takes place at the initial 12
ACS Paragon Plus Environment
Page 12 of 32
Page 13 of 32 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
stages of the simulations when the systems are being equilibrated. We never observed that an ion moved from the central to the upper vestibule. Eventually, an equilibrium situation is reached in which there are between four and five Na+ inside the extracellular domain. Some of them stay in the central vestibule and the rest is in the extracellular vestibule. Ions in the extracellular vestibule eventually exchange with those of the external solution or the central vestibule. Under these circumstances, Cl− can also enter the extracellular vestibule. Table 1 presents the probabilities for the alternative occupation numbers of the upper part of the receptor, for Na+ and Cl− . Movies 1 and 2 of the Supporting Information section display the behavior of Na+ and Cl− in standard MD simulations of the closed and open conformations, respectively. Table 1: Percentages for the alternative occupation numbers in the upper part of the receptor. Ions in the central vestibule, extracellular vestibule and upper half of the pore were taken into account. Nions Closed Na+ Cl− Open Na+ Cl−
0 0.0 96.5
1
2
3
4
5
6
0.0 6.5 13.9 52.8 26.4 0.4 3.5 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 69.4 27.5 2.5
0.8 57.9 41.0 0.3 0.6 0.0 0.0 0.0
The quantities in Table 1 demonstrate the preference of the interior of the extracellular domain for cations. This result agrees with the proposal of Kawate et. al., 31,42 who suggested that the negative charges of this region help to maintain inside a higher concentration of cations than anions. However, as stated in the a previous paragraph, several Na+ that enter the receptor are not available to move toward the intracellular side. Instead, they find a much more stable location at the upper part of the central vestibule. Consequently, it is not clear that the larger number of Na+ found inside can be used to justify the preference of the receptor for the transfer of cations. The situation becomes clearer after analyzing the distributions of Na+ and Cl− along the Z-axis. These distributions are presented in the upper panel of Fig. 2. They were computed with Eq. 1, employing the data collected in 13
ACS Paragon Plus Environment
The Journal of Physical Chemistry
ρ(Z; r