J . Phys. Chem. 1985, 89, 2870-2876
2870
Na+ and K+ Ion Transport through a Solvated Gramlcldln A Transmembrane Channel: Molecular Dynamics Studles Using Parallel Processors Kwang S. Kim, H. L. Nguyen, P. K. Swaminathan, and E. Clementi* IBM Corporation, Kingston, New York 12401, and National Foundation for Cancer Research, Bethesda, Maryland 20814 (Received: November 16, 1984; In Final Form: March 4, 1985)
The energetics, structure, and dynamics of two different cations (Na+ and K'), water, and gramicidin A (GA) systems are studied by molecular dynamics methods with parallel computationaltechniques. The parallel computation techniques employed in this study are described. In our simulations, we study the dynamical behavior of 8 1 water molecules and a single ion of each kind, placed at one of two different locations along a channel. The ion and water move under the dynamic force fields from each other as well as the static force field of the channel, which is compared of the GA framework constrained to remain the Urry dimer conformation. All force fields are of ab initio origin. The trajectories illustrate the characteristic time scales of typical events at the chosen locations. The velocity and force autocorrelation functions from the simulations (the first computed for any GA model) are discussed.
I. Introduction The problem of ion transport in biological channels such as gramicidin A (GA) is central to understanding many current important issues in chemistry and biology. Although experimental and theoretical progress is being made toward elucidating the structure and dynamics of ion transport through GA channels, we have still not reached a detailed understanding of atomic and molecular motion in the G A channels. For recent experimental reviews, see Urry,' Ovchinnikov,2 Eisenman and Horn,3 and Anderson.* GA is uniquely suited to theoretical studies of the underlying molecular basis for ion transport due to its relatively simple and well-defined structure and a wealth of experimental data. Several conformational studies investigating the local minima of the GA structure in the channel have been r e p ~ r t e d ,but ~ , ~these energetic studies were performed with only selected torsional degrees of freedom of GA in the absence of ions and water molecules. Also, though other electrostatic calculations7 have reported activation barriers of 20-30 kJ/mol, these models (treating the membrane as a low-dielectric continuous medium) are probably too simple to yield predictive information at the molecular level. Among the earlier treatments are molecular dynamics (MD) studies on models chosen for GA itself! yet these simulations imposed restrictions such as neglecting water molecules and constraining the ion to move only along the channel axis. More recently, Wilson and his collaborators9 used M D techniques to study the interaction of GA with cations and eight water molecules. Their model used vacuum rather than bulk water outside the channel, thereby neglecting solvation contributions from bulk water. Also of concern are the uncertainties in the assumed potentials which, as pointed out by Wilson et al.,9 constitute a shortcoming of that study. Other recent theoretical studies of energy profiles of Na+ interacting with GA have been reported by Pullman and Etchebest,lo wherein temperature, dynamics, statistical effects, and solvation influence have bee ignored, and the binding energy of Na+ and GA was computed by using parameters which are chosen to reproduce the results of ab initio calculations on small systems. In a recent series of papers we reported on ab initio 6-12-1 atom-atom pair interaction potentials of water-GA," Na+-GA,12 and K+-GA,I3 and the Monte Carlo (MC) simulations for water-ion-GA These M C simulations were undertaken to investigate (1) the number of molecules which can be placed inside the GA channel, (2) the energy profiles of a Na+-water-GA system, (3) the structure of the water molecules inside the channel in the presence of Na+ and K+ ions, (4) solvation energetics of ions inside and outside the channel, ( 5 ) the most probable path of ions in the channel, and ( 6 ) the effects of an *Address correspondence to this author at IBM
0022-3654/85/2089-2870$01 SO10
externally applied electric potential along the channel. Libration effects of the carbonyls and other dynamical structural changes of GA, as well as the role of phospholipid^,^^ are currently being investigated in our laboratory. Considering the large number of underlying elementary processes in the ion transport through GA, there appears to be a need to sort out their respective roles. It is natural to proceed through a sequence of models such that the behavior of the real channel is approached gradually in stages as realism is increased. Ideally, in theoretical studies, one should also be able to attain some predicting power; thus, from our sequence of models, we should be able to ascertain the consequences of changes in an important set of control variables. There is much speculation on what the control variables may be in GA, including (1) the role of GA in opening up the dielectric barrier of the membrane, (2) the role (1) (a) Urry, D. W. In "Membranes and Transport", Vol. 2, Martonosi, A. N., Ed.; Plenum Press: New York, 1982; pp 285-294. (b) Wry, D. W. Top. Current Chem., in press. (2) (a) Arseniev, A. S.;Bystrov, V. F.; Ivanov, V. T.; Ovchinnikov, Y. A. Fed. Eur. Biochem. SOC.1984,165,51-56. (bl Ovchinnikov, Y.A.; Ivanov, V. T. In 'Conformation in Biology", Srinivasan, R., Sarma, R. H., Eds.; Adenine Press: New York, 1983; pp 155-174. (3) Eisenman, G.; Horn, R. J. Membr. Biol. 1983, 76, 197-225. (4) (a) Andersen, 0. S . Annu. Rev. Physiol. 1984,46, 531-548. (b) Levitt, D.G . Biophys. J . 1982, 37, 575-587. (5) Venkatachalam, C. M.; Urry, D. W. J. Comput. Chem. 1983, 4 , 461-469. (6) (a) Fraga, S.; Nilar, S. H. M. Can. J. Biochem. Cell Biol. 1983, 61, 856-859. (b) Venkatachalam, C. M.; Urry, D. W. J. Comput. Chem. 1984, 5, 64-71. (c) Ramachandran, G. N.; Ramachandran R. Ind. J . Biochem. Biophys. 1972, 9, 1-11. (d) Popov, E. M.; Lipkind, G. M. Mol. Biol. (Moscow) 1979,13,363-376. ( e ) Venkatram-Prasad, B. V.; Chandrasekaran, R. Int. J. Peptide Protein Res. 1977, 10, 129-138. (7) (a) Parseigan, V. A. Ann. N.Y. Acad. Sci. 1975, 264, 161-174. (b) Levitt, D. G. Biophys. J . 1978,22,209-219. (c) Jordan, P. C. Biophys. Chem. 1981, 13, 203-212. (8) (a) Brickmann, J.; Fischer, W. Biophys. Chem. 1983,17,245-258. (b) Fischer, W.; Brickmann, J.; Lauger P. Biophys. Chem. 1981,13, 105-116. (c) Fischer, W.; Brickmann, J. In "Bioelectrochemistry" and Bioenergetics: Proceedings of the International Meeting, Physical Chemistry of Transmembrane Ion Motions"; Paris, 1982. (d) Schroder, H.; Brickmann, J.; Fisher, W. Mol. P h p . 1983,11, 1-11. (9) Mackay, D. J.; Berens, P. H.; Wilson, K. R.; Hagler, A. T. Biophys. J . , in press. (10) (a) Pullman, A.; Etchebest, C. Fed. Eur. Biochem. SOC.1983,163, 199-202. (b) Etchebest, C.; Pullman, A. Fed. Eur. Biochem. SOC.1984, 170, 191-195. (c) Etchebest, C.; Ranganathan, S.;Pullman, A. Fed. Eur. Biochem. SOC.1984. 173. 301-306. (1 1) Fornili, S . L.. Vercauteren, D. P.; Clementi, E. Biochim. Biophys. Acta 1984, 771, 151-164. (12) Kim, K.S.;Vercauteren, D. P.; Welti, M.; Fornili, S. .; Clementi, E. IBM Technical Report Pok.42, April 20, 1984. (13) Kim, K. S.;Vercauteren, D. P.; Welti, M.; Chin, S.; Clementi, E. Biophys. J . 1985, 47, 327-335. (14) Kim, K. S.; Clementi, E. J. Am. Chem. SOC.,submitted for publication. (15) Clementi, E.; Corongiu, G.; Detrich, J.; Chin, S.; Domingo, L. I n t . J . Quantum Chem. 1984, 18, 601-617.
0 1985 American Chemical Society
Na+ and K+ Ion Transport of water in bulk form as well as in the form of a chain within the channel solvating the ion, (3) the role of fast collective translations and/or rptations/tumbling motions of the water molecules, (4) the effect of membrane (as modeled by phospholipids) on the solvent and possibly the direct coupling to the GA dimer, ( 5 ) the role of the bending of the head parts (the center of the dimer), the motion of hydroxyl oxygens of the tail parts, and the three C=O groups librating at each one of the two mouths of the channel, a motion ‘that reflects the flexibility of the complexly networked hydrogen-bonded helical structure and indeed the possibility that thgse motions may be influenced by coupling to the phospholipid motions, (6) the role of the entire channel due to intra- and intermolecular interactions of GA, like a Snake swallowing food, and (7) the effects of electrical potential and ion concentration on the ion transport. When the mechanisms mentioned above are considered, many possibilities arise. It is clear that, by proceeding via reasonable control variables within the larger umbrella of the concept of a “linked set of models”, we might examine the variables by bringing them in an orderly manner, going from the simple to the complex. A major part of our effort will be to define the variables more precisely. At the same time this will provide clues as to the more elaborate intramolecular potential calculation needed for GA to bring in its relevant motion. The way in which to add complexity to a model becomes apparent from the results of a simpler level. This research project, although clearly most demanding in computational facilities, is feasible because of the parallel computation s e t ~ p . ’ ~InJ ~addition, it is a challenging project in model development. Here we present as a first step an MD simulation of the motion of the ion and water molecules, with the restriction that all the GA atoms are clamped a t the configuration provided by Urry’s atomic coordinate^.^ With this choice, only the static and instantaneous solvent response due to the GA atoms on the ion motion are included. The structural and full dynamical solvent effects due to water molecules on the cation motion are also included. Static and dynamic solvent effects from the phospholipds are neglected in this study due to the uncertainty of the location of the phospholipid atoms. Rather, we impose a “softwall” effect outside the GA channel, as described in section 11. We note here that GA structural changes are strongly anharmonic;” thus, if a calculation is performed without properly considering the anharmonic terms, it might even lead to spurious results. Although the GA structural changes due to intramolecular motions can be significant in the real system, there is little reason to believe that these effects would be dominant because the carbonyls may not librate too much inside the channel except a t the binding sites. This argument can be supported by the strong hydrogen bonding present between carbonyl oxygens and amide hydrogens in GA13I4 and also by the strongly anharmonic potentials for carbonyl librations. Although ref 9 shows the carbonyls librate strongly inside the channel, the potentials used did not consider the strong anharmonicity of the potentials due to the change of the conformations of amide bondings and carbonyl librations. On the other hand, although the selection of one GA configuration necessarily imposes a limitation on our findings, it does provide information on the extent of significance of the neglected features in our model. 11. Theoretical Approach For the simulations treated in this paper, the intermolecular water interaction is assumed to be pairwise additive and modeled by a rigid MCY potential.18 This ab initio pair potential is based on a quantum mechanical calculation of a configuration interaction study on the water dimer. The ion-water interaction potential (16) (a) Corongiu, G.; Detrich, J. H. IEM J . Res. Deu. in press. (b)
Clementi, E.;Corongiu, G.; Detrich, J. H.;Khanmohammadbaigi, H.; Chin, S.;Domingo, L.; Laaksonen, A:; Nguyen, H.L. In “Structure and Motions: Membranes, Nucleic Acids, Proteins”, Clementi, E., Corongiu,G., Sarma, M. H., Sarma, R. S., Eds.; Adenine Press: New York, 1985; pp 49-85. (17) Kim, K. S.; Clementi, E. unpublished data. (18) Matsuoka, 0.; Clementi, E.; Ycahimine, M. J. Chem. Phys. 1976,64, 1351-1361.
The Journal of Physical Chemistry, Vol. 89, No. 13, 1985 2871 is also taken to be pairwise additive; for these we used 6-12-1 potentials fitted to a simplified three-site model” using the S C F results of Kistenmacher et at.19 The intermolecular interactions between GA and cations and between GA and water are approximated by the painvise sum of the atomatom pair interactions of the well-known 6-12-1 form.11-13*z0 Our concern is, first, to find out energy transfer patterns between the water and the ion (especially when the Gramicidin framework is assumed rigid). Since the energy-transfer mechanisms of GA to the ion can arise via water or directly (and indeed will depend on where the ion is in the channel), our knowledge of direct water-ion patterns in the GA static field will tell us later what part is due to GA motions (such as the intramolecular problem of GA modes or, more remotely, via GA but from phospholipids). In this way we can also judge the role of selected modes of the surroundings of the ion by turning them on/off systematically. In our MD simulation an equilibrated Monte Carlo configuration is used as a starting configuration. Here 81 water molecules are placed in a cylindrical volume whose length and radius are 48 and 5.5 A, respectively. Translation symmetry constraint along the z axis of the cylinder is imposed for the water molecules. The x and y components of the center-of-mass (CM) velocities of the water molecules near and outside the surface of the channel cylinder are influenced by the “soft” forces. This is a special kinematic condition treating the outside of the cylinder as continuous media. This boundary condition means that any water molecule near or leaving the cylinder volume is pushed inside the cylinder. Following this constraint, we employed a waterlike cylindrical softwall interaction potential of the form V(r) = C8(r - r0)*,where r is the radial distance for each water molecule, and c8 and ro are constants chosen as 60 (kJ/mol)/(A)s and 5 A, respectively. This constant partially simulates water-water repulsive interaction. The total interaction energy of the studied system consists of the water-water, ion-water, water-GA, ion-GA, and watersoft wall interactions. The Verlet leapfrog algorithm*l was used to solve the equations of motion both for the translation and the rotation. Quaternion algebraZZwas used to describe the molecular orientations of rigid water molecules. An integration time step of 0.0002 ps was used. Our approach to studying the “solventlike” influences of GA on the ion, by systematically introducing motion starting from a motionless GA to begin with, as done here, resembles the overall philosophy adopted in the clamping molecular dynamics method.23 This method has a mathematical justification and is formally analogous to the gas-phase Born-Oppenheimer approximation used in quantum structural calculations. For a given set of constrained solvent coordinates, in performing molecular dynamics simulations, one must study the motion of the solute and solvent to determine the dynamic response of the solvent to solute motion. By successively relaxing the set of the constrained solvent coordinates, the partial and final full solvent dynamic responses are constructed from MD trajectory calculations. The objective of this method is obvious: specific partitioning of the energy flows by successively relaxing the elements of the set of the constrained coordinates. Furthermore, these calculations will reveal other interesting informations, for example, the dependence of the dynamic and static solvent response on the ion positions along its reaction path in the GA channel.
111. Parallel Computation Parallel computational methods are different from conventional sequential methods in both hardware and software. We used a (19) Kistenmacher, H.; Popkie, H.; Clementi, E. J. Chem. Phys. 1973,59, 5842-5848. (20) Corongiu, G.; Clementi, E. Gazz. Chim. Itol. 1978, 108, 273-306. (21) (a) Verlet, L1. Phys. Reu. 1967, 165, 21. (b) Fincham, D. ‘Information Quarterly-CCPS”; No. 2, Sept 1981, p 6. (22) (a) Evans, D. J.; Murad, S. Mol. Phys. 1977, 34, 327-332. (b)
Goldstein, H. ‘Classical Quantum Mechanics”; Addison-Wesley: Reading, MA, 1980, 2nd ed. (23) Adelman, S.A. J. Chem. Phys. 1984,81, 2276.
2872
Kim et al.
The Journal of Physical Chemistry, Vol. 89, No. 13, 1985
high-performance computing system consisting of a distributed system of IBM host computers (one 4381 and two 4341’s) with several attached processors (FPS- 164). The system is essentially a multiinstruction stream, multidata stream (MIMD) system,24 in the form of a distributed network of nodes. Each node of the network is a full processor capable of executing complete programs. It is important that a task be distributed equally among the processors. Thus, a slower processor should have less to work on, and a processor having a different data transfer time (for example, because of a lower/higher capacity channel interface) should receive an appropriately adjusted share of the task. The molecular dynamics method for calculating the evolution of a system of N interacting atoms and molecules consists of two basic steps, force evaluation followed by numerical integration of the equations of motion. Most computational time is spent in the evaluation of potential energies and site forces (and site vinals). For the simulation system studied in this paper, the total interaction consists of the water-water interaction term, the ion-water and the watersoftwall terms, and the water-GA and the ion-GA interaction terms. Considering the pairwise terms used to compute these interactions, the number of pairs to be evaluated at each time step is N ( N - 1)/2, where N is the number of atoms or molecules involved in the calculation of the interaction pairs. For the simulated system of 512 GA atoms and 81 water molecules (the three-site water model is used) we need about 6 mega floating operations for a single time step to evaluate the pair interactions between GA atoms and water. If the model includes several phospholipids (the potential of which has been recently calculatedz5),thousands of water molecules are required due to the increased volume. Also, all the interaction pairs involving water, ion, GA, and phospholipids would enormously increase the computation time. The parallel code has four modules that execute in parallel modez5 (i) to obtain uniformity of subtask allocation (i.e., a task should be distributed equally among the processors) and (ii) to decrease the number of data transfers initiated during the execution of the multitasked codes, while at the same time (iii) minimizing the floating-point computations involving integer operations of extensive logic and indexing in the modules that evaluate the potential and force of interaction. The first module evaluates the water-water interaction term, the second module computes the ion-water and watersoftwall interaction, and the third parallel module calculates the water-GA and ion-GA interaction. After these interaction terms are computed, the parallel integration module is called to numerically integrate the equations of motion. We now arrive a t the following scheme for the multitasking of our “parallelable” MD algorithm. (1) At the beginning of each time step, the primary process (PP) sends to each secondary process (SP) the coordinates and velocities of the previous phase space point. (2) The primary and the secondary processes then send these data to the attached processors. Each attached processor (AP) evaluates its assigned portion of the site-site pair potential, site force, and site virial for the water-water interaction term. Table look-up evaluation of the potential, site force, and site virial are used. Upon completion, each AP transfers its portion of the potential, site force, and site virial for the water-water interaction term back to its corresponding process. The secondary processes then send to the primary process their portions of these calculated data, where elements calculated by each SP are put together. (3) The PP sends to each SP a dummy set of data. These data are sent to the AP’s. The AP’s then evaluate the energy, site force, and site virial for the ion-water and water-softwall interactions term. Upon completion, each AP transfers its portion of the potential, site force, and site virial for the ion-water and watersoftwall back to its secondary process. The secondary processes then send their received portions of the calculated data to the PP (24) Kogge, P. M. “The Architecture of Pipelined Computers”;McGrawHill: New York. (25) Swaminathan, P. K.; Vercauteren, D. P.; Kim, K. S.; Clement, E. J . B i d . Phys., submitted for publication.
where elements calculated by each SP are put together. (4) Similar procedures outlined in step 3 are repeated here to evaluate the potential, site force, and site virial for the water-GA and the ion-GA interaction term. The computed data are now sent from the AP’s to their corresponding processes. The SPs then send their portions of the computed data back to the PP, which puts these data together. (5) The PP sends to each SP a portion of the interaction site force. These data are sent to the AP’s. Each AP converts the site forces into CM force and torque, and then performs the integration of the equations of motion to find the coordinates and velocities of the phase space point. Various quantities, such as mean-square displacement, species potential energies, kinetic energies, step temperature, pressure, etc., can also be computed. (6) Periodic output may be generated at this point, and the cycle repeats from the first step for N numbers of time steps. When the cycle stops, final output for average values of computed quantities is done by the primary process. The splitting of the MD “parallelable” algorithm is the following. The pairwise potential energy matrix is the usual symmetric matrix; the site force matrix is a special skew-symmetric matrix in which its diagonal elements can be nonzero. Uniformity of subtask allocation and the double loop structure of the algorithm used in this routine suggest a sieving breakup of the following form (sequential)
(parallel)
where NAP is the number of attached processors used and the function Aid) denotes the pairwise potential energy or the corresponding force between molecules (atoms) i and j . The portion enclosed in parentheses on the right side of eq 3.1 is the part that runs on the mth process, with the summation over j stepping by NAP, rather than one. A slightly different modification scheme is taken for multitasking the numerical integration of the equation of motion. The single deloop structures appearing in this module are each broken up into equal portions to be taken as subtasks. Here we divide the terms among the processes according to the relation (sequential
(para I le1 1 (3.2)
i =1
m=l i=m.NAP
where the single sum over i on the left side of eq 3.2 is divided into a set of NAP approximately equal partial sums. On the right side of eq 3.2, the summation is performed over i stepping by NAP unless i > N . The summation over i on the right side of eq 3.2 is the part that runs on the mth process. When table look-up is done for evaluation of the site potential, site force, and site virial (in the second step), the following procedure is adopted to save transmission overhead. Once created, the tables are invariant throughout the simulation run. For this reason they are created at the beginning of the run, transfered to the AP’s, copied into local common blocks, and saved for subsequent use. Similarly, if a set of data is required by the AP’s for the computation, and, in particular, these data remain unchanged throughout the calculation, a local array is created in each AP and kept there for subsequent use. Attention should be given to the fact that the subtasks are executed in different processors. As such, their only relation to the main task is data communication. Once this link is established, the subtasks may arrange the data in any desired fashion; that is, as long as data communication is consistent, subtask optimization is completely independent of the main task. The parallel MD approach has a number of attractive features, including the fast floating point speed achieved by parallelism, which is essential, because MD simulations of systems with hundreds or thousands of atoms (molecules) are highly CPU bound. The parallel approach would make previously infeasible
Na+ and K+ Ion Transport
The Journal of Physical Chemistry, Vol. 89, No. 13, 1985 2813
TABLE I: Energetics for the Water-Cations (Na+ and K+)-GA System" W-W
w-s w-I
W-GA I-GA
PE temp, K
NaO
Na15
KO
K15
-25.01 f 0.46 0.09 f 0.05 -5.25 f 0.13 -50.33 f 0.30 -277.21 f 2.63 -83.91 f 0.06 303.6 f 18.3
-24.55 f 0.52 0.08 f 0.04 -9.37 f 0.22 -49.56 f 0.32 -160.88 f 3.20 -85.39 f 0.04 297.3 f 13.5
-25.14 f 0.37 0.21 f 0.06 -3.71 f 0.16 -50.25 f 0.26 -209.89 f 3.00 -81.48 f 0.05 303.1 f 14.7
-24.97 f 0.43 0.15 f 0.05 -5.86 f 0.27 -49.50 f 0.33 -183.79 f 3.41 -82.45 f 0.04 299.6 f 14.0
"Interaction energies between water and water (W-W), between water and softwall (W-S), between water and cations (W-I), between water and GA (W-GA), and between cation and GA (I-GA). Units are kJ/mol of water except for I-GA, the unit of which is kJ/mol of cation.
tx
I
Monomer-A&B
Monomer-B
Monomer-A
100.0
Figure 1. Features of a GA membrane channel. Projections of a GA channel (GA dimer) onto x-y plane (left) and x-z plane (right).
-15 1
I
H,O/GA
NaO
t
Figure 3. Temperature variation of NaO in GA.
5 -30 E
calculations (for example, studies of biomolecular structure-dynamics-function relationships in biological systems) a practical reality.
W
- 75 5 1-16 I
0
-8
I
-20 I
I
8 .16; ----+ZIAI ~
0 T
1
I
1
I
I 20 1
I
GA
6.3
X -6.3 Figure 2. Minimum interaction energies of water, K+,and Na+ with GA inside the channel as function of z (top) and the isoenergy maps for water, , 'K and Na+ in the x-z projection. (The carbonyl oxygens of GA are darkened.)
IV. Results As a visual aid in reporting our results, we begin by showing the GA structure and isoenergy maps for interaction with water and cations. Figure 1 shows the x-y and x-z projections of GA. This figure explicitly shows the small channel cavity. Figure 2 shows our previous results14for the minimum ab initio interaction energies of water, K+, or Na+'with GA as a function of z and for the isoenergy maps in the x-z plane, comparing with the backbone structure of GA. The minimum interaction energies were obtained by allowing the cation to move to the most energetically preferred position along the x and y directions for the given z value. For ion-water-GA systems experimental binding sites are located near Z = f 11 Although these binding sites are not shown in Figure 2, they may be explained theoretically by considering ion-water and ion-GA interaction energy contributions and also by including carbonyl librations and tail movement^.'^^^^ Next, we present and discuss the current MD results. For each cation, two different cases are simulated: one with the ion initially located near the center of the channel, and the other with the ion initially placed near the outside of the channel mouth. The former is around z = 0 A and the latter is around z = 15 A. To facilitate the following discussion, the former case for Na+ is denoted by NaO; the latter case for Na+, by Na15. For the same reasons, the former case for K+ is denoted by KO; and the latter case for K+, by K15. The temperature variation for the MD simulation NaO is shown in Figure 3. The fluctuation is reasonable for 81 water molecules a n d one cation; its standard deviation is only 18 K, while t h e average temperature is 304 K. Other simulated cases were as stable as shown in Table I; for Na15 the temperature variation is 297 f 14 K; for KO it is 303 f 15 K; and for K15 it is 300 f 14 K. Table I also shows the energetics for the water-cation-GA systems, reporting the interaction energies for water-water, (26) (a) Urry, D. W.; Walker, J. T.; Trapane, T. L. J . Membr. Biol. 1982, (b) Urry, D. W.; Venkatachalam, C. M.; Spisni, A.; Bradley, R. J.; Trapane, T. L.; Prasad, K. U. J. Membr. Biol. 1980, 52, 29-51. 69, 225-231.
2814
Kim et al.
The Journal of Physical Chemistry, Vol. 89, No. 13, 1985
Nal5
NaO
K15
KO
lorn
Figure 4. Ensemble-averaged positions for water molecules and cations. Cation positions are denoted by circles.
1
om - l ",
d
E0
.J
0 . 6 o, 7
z-15
00
0
KO
iJ ; I ob,-.------; R
a
0 0
0 0
u
0 -1
-1
01
03
T (P-c)
i , , 03
01
T (p*ec)
Figure 5. Trajectories of cations in GA.
water-softwall, water-cation, and water-GA, including total interaction potential energies. These values are very similar to the previous MC results.14 This is because the introduction of softwall instead of hardwall a t the cylindrical radial boundaries is of negligible difference to the energetics. The ensemble averaged positions of water molecules and the cations Na+ and K+ are shown in Figure 4; the cations are darkened. Inside the channel, about nine water molecules are found to be linearly positioned. Since our model includes many water molecules outside the channel, the water molecules near the channel mouth among those nine water molecules are not oriented toward the cations; this result differs from Wilson et al.? wherein the water molecules outside the channel were neglected. The water molecules outside the channel are more disordered, showing a three-dimensional network characteristic of the first and second solvation layers. Figure 5 shows trajectories of the cations Na+ and K+ near two selected points, near the center of the channel and outside of the channel. Within a 0.5-ps window the coordinates do not change much; the amplitudes are less than 0.6 A for either of the ions with the exception of the K15 case, where the amplitude is about 0.8 A. This latter phenomenon indicates more diffusivity for K+. The phase angle variation [tan-'b ( r ) / x ( t ) ) ]with time shows some oscillation, and radial distance (R) seems to change a little, but the radial distance per angular variation seems to be very characteristic. Inside the channel Na+ has R = 0.9 f 0.2 A, while K+ has R = 0.7 f 0.2 A. This result agrees with our previous MC simulation results;14 Na+ executes a more helical path motion along the channel than K+, since Na+ is more strongly bound to
NaO
I
0.1
0.3
Na15
I
0.1
T (pHc)
0.1
0.3
1
I
, T (wc)
Figure 7. Force variation of cations in GA. (Units for Nat and K+ are 3.30 X and 4.33 X N, respectively.)
GA and its van der Waals radius is smaller. Figure 6 plots the velocities with respect to time, while Figure 7 plots the forces with respect to time in Figure 7 to indicate the
Na+ and K+ Ion Transport
The Journal of Physical Chemistry, Vol. 89, No. 13, 1985 2815
Na' I Ionic Solnl
? - 0
\ A
O E
KO
\
[Ionic Solnl
n
t ?
n
-0 E
-051
Figure 8. Velocity autocorrelation functions of cations in GA compared with those of cations in ion solution.
characteristic time scales for changes in these key variables. For Na+ the forces are more irregular and often show strong peaks, while for K+ the forces are rather slowly varying and regular; this is related with the fact that K+ moves around near the channel axis more than Na+ does. The Na+ ion is physically closer to the GA frame and moves under the influence of a tighter cage, thus showing higher freqeuncy components (relative to the K+ ion) in its motion as seen from the trajectories, velocities, and forces. It must be remembered that the short trajectories shown cannot be representative of the rich possibilities that exist as time goes on; this is especially true in regions separated by potential barriers that have the universal effect of causing transmigrations to be relatively infrequent events, not showing up in short time spans. However, the time frame shown is representative of the unhindered regions that are local and accessible. Indeed, any flexibility in the x and y directions (toward the channel inside walls) must come from a much more elaborate analysis of the intramolecular modes of GA frame. However, the strongly caged state of the ion in the center region of the channel will be rendered more flexible along the channel axis when an external field is present as is in the real channel, for example. In the latter event, the ion is expected to travel down the channel driven by the field and be delayed only in regions that separate positions characterized by cages of very different strength such that the resultant barriers are not overcome by the field strength over that length. In Figure 8, the velocity autocorrelation function (VACF) for the Na+ and K+ cations in ion ~olution.~'The Na+ (NAO) and K+ (KO) cations inside the channel show fast oscillations of small damping with memory of the initial velocity lasting for a somewhat longer time. This reflects large local potential barrier that the cation motion inside the channel is strongly caged under the strong interaction force field from the GA and neighboring water molecules. The cations outside the channel (NA15 and K15) show some characteristic similar to the cations in bulk water (NaO and Na15). The behavior of Na15 and K15 VACF's seems partly to indicate that the inclusion of bulk water in this simple model is qualitatively reasonable in terms of number of water molecules outside the channel as well as the softwall description of its (27) (a) Impey, R. W.; Madden, P. A.; McDonald, I. R., in press.
t
,
LL
Y
-0 5
0 05
0 15
T (pa=)
0 05
0 15
T (PW
Figure 9. Force autocorrelation functions of cations in GA.
boundary region. When K+ is compared with Na+ from the six insets of Figure 8, K+ shows more diffusivity, losing its memory of previous velocity somewhat faster. Similarly, Figure 9 shows the force autocorrelation functions. These support the preceding conclusion pointing to the caged behavior. The time-dependent force autocorrelation function is related to the time-dependent friction in solution. The behavior of this function sheds light on the solvent effect on ion motions along different directions. For a harmonic oscillator that never loses its energy, this function is also a periodic nondecaying function repeating in time. Any sign of damping reflects dissipation of energy from the ion to the medium. A slight damping is noticed in Figure 9, but the oscillations reflect an overall slightly anharmonic coupled motion together with the surrounding cage in the GA channel. The contrast of the force autocorrelation functions between the inside and the outside of the channel in Figure 9 is once again striking as in Figure 8, reflecting that the
2876 The Journal of Physical Chemistry, Vol. 89, No. 13, 1985
t LL v
c \
r-
NaO(zz) A
A
t 5 -
=
0 5-
. L
9
\ -0 5 v
I
1
0.05
0.15
T (PW
Figure 10. Force autocorrelation functions of cations in GA for the x x , y y , and zz components.
cation inside the channel is caged. Figure 10 shows the components of force autocorrelation functions for the Na+ ion along and perpendicular to the channel axis. The zz component of NaO along the axis of the channel behaves very similarly to the x x and y y components for the NaO case shown; similar behavior is also observed for KO. This indicates that the motion of the water molecules does not couple strongly to the motion of the ion, possibly because the water molecules near the cation are strongly bound to the G A carbonyl oxygens. Alternately expressed, the interactions between water molecules by either dispersion forces of dipole-dipole forces fails to be a significant mechanism for energy dissipation from the ion by being carried into the bulk water. It should be mentioned, however, that such energy transfer is not impossible within present results, as an infrequent event. The origins of the force field causing the strong caging discussed above are clear from our model. From our simulations, we can conclude that the local barriers around the ion inside the channel are due to (1) the two strongly bound water molecules on either side of the cation (each held in position by carbonyl oxygens), (2) a deep potential well due to the strong interaction of the cation with GA carbonyl oxygens, and (3) the repulsive interaction of the cation with the GA channel wall.
V. Concluding Remarks The calculations of the dynamical quantities yield the insight necessary to determine the way in which the various control
Kim et al. variables considered in the sequence of models participate. The fact that, in the absence of an external field, strongly caged motions persist for long along the z axis is an important conclusion. This was our reason to study the system with no field applied at first. Without it one cannot be certain of the importance of speculated mechanisms acting under different conditions. As the models grow more complex, it becomes necessary to explore more aspects in more detail. It is expected that a sufficiently strong applied field may be used to induce motions along the channel even in this simple model. However, since unrealistic effects due to strong fields are of no interest, improved descriptions of the energy profiles must be achieved before examining field effects. It is reasonable at this stage to speculate that the GA backbone motions, head and tail motions, and carbonyl librations are important for the next stage in the building process presented here. The model should also include the phospholipids, whose interaction with the cations and bulk water play an important role. A reasonable model incorporating all these terms will constitute a more rewarding example for extensive study and must also investigate the infrequent events involved. Within our limited study, mean displacements or mobilities of ion and water inside the channel do not seem to be too different. At this stage, it does not allow us to conclude if an ion will push a column of water ahead of it in order to pass through the channel, because such a process requires a time order of magnitude larger than the time used in our simulation. Techniques of extending the range of time scales being described are crucial in this problem since there are limits on the number of reliable integration steps. For the GA channel system, the classical kinetics may not be applied in a simple manner. It should be noted that the experimental free energy profiles28obtained by applying Eyring's rate theoryz9to the channel transport are q u e s t i ~ n a b l e . ~The ~ presence of strongly caged oscillations observed in the present model forewarns the timedependent frictional forces are an important aspect of this problem. Acknowledgment. P.K.S. thanks the National Foundation for Cancer Research (NFCR) for financial support. We thank Dr. A. Laaksonen, Dr. M. Probst, and Prof. G. Lie for helpful discussions. Registry No. Na, 7440-23-5; K, 7440-09-7; gramacidin A, 1102961-1.
(28) (a) Hladky, S. B.; Haydon, D. A. Biochim. Biophys. Acra 1972,274, 294-312. (b) Eisenman, G.; Sandblom, J. P. Biophys. J . 1984, 45, 88-90. (c) Bambergnand, E.; Lauger, P. Biochim. Biophys. Acra 1974,367,127-133. (29) (a) Eyring, H.; Lumry, R.; Woodbury, J. W. Res. Chem. Prog. 1949, 10,1W114. (b) Woodbury, J. W. In 'Chemical Dynamics: Papers in Honor of Henry Eyring", Hirschfelder, J. 0.;Ed.; Wiley: New York, 1971; pp 601-661. (c) For the application to channel permeation, see for example Hille, B. In "Membranes: A Series of Advances",Vol3, Eisenman, G., Ed.;Marcel Dekker: New York, 1975; pp 255-323. (30) (a) Lauger, P.; Stephan, W.; Frehland, E. Biochim. Biophys. Actn, 1980, 602, 167-180; (b) Levitt, D. G . Biophys. J . 1982, 37, 575-587.
