Using Theory and Simulation to Understand Permeation and Selectivity in Ion Channels

METHODS: A Companion to Methods in Enzymology 14, 342–351 (1998) Article No. ME980589

Using Theory and Simulation to Understand Permeation and Selectivity in Ion Channels Eric Jakobsson Department of Molecular and Integrative Physiology, Center for Biophysics and Computational Biology, National Center for Supercomputing Applications, Beckman Center for Advanced Science and Technology, University of Illinois, Urbana, Illinois 61801

It is clear that the function of ion channels must flow from their structure. With recent advances in computational power and methodology, it appears feasible to correlate structure to ion channel permeation at an atomistically detailed level of description. The overall strategy is to structure the calculations in a hierarchy, ranging from coarse-grained thermodynamic and kinetic descriptions to fine-grained molecular dynamics descriptions with atomic detail. Each level of description is connected to the others by appropriate statistical mechanical theory. The coarse-grained descriptions can be correlated directly with electrophysiological experiment. The fine-grained descriptions are used to parameterize the coarse-grained descriptions and to describe the permeation process at the most detailed level. This strategy has so far had varying degrees of success. It has successfully described water permeation through lipid bilayers and gramicidin channels. It has revealed the essential events of ion permeation through gramicidin channels at an atomistically detailed level. The role of channel protein motions in permeation has been elucidated. However, it appears that force fields used to describe molecular dynamics must be refined further to achieve completely accurate predictions of the permeation of such small ions as sodium. Channels with more complex structure and more multiion occupancy than gramicidin pose major computational challenges with respect to sampling protein conformations and ion distributions involved in the permeation process. Possible approaches to meeting these challenges are discussed. q 1998 Academic Press

Molecular biophysicists strongly believe that function is a consequence of molecular structure. In principle it should be possible to predict the electrophysiology of an ion channel precisely if one knows the structure, plus the molecular details of the channel’s environment, i.e., its relationship to the surrounding membrane and electrolyte. However, even such a small biological system as a single ion channel is a formidably complicated structure with many degrees of motional freedom. The smallest known ion channel, gramicidin,

contains just 30 amino acids, and has a precisely known structure. With the fastest machines at this writing, to do a direct molecular dynamics simulation of a gramicidin channel in a lipid environment requires approximately 1 CPU h of computer time for 1 ps on a supercomputer. To do a ‘‘brute force’’ simulation of ion permeation; i.e., to see enough ions go through the channel to get a direct measure of the current for one point on an I–V curve, one would need to simulate about one-tenth of a microsecond, or about 100,000 h of CPU time. This is obviously an impractically large amount of time with present technology and with technology available in the foreseeable future. To overcome this obstacle to direct simulation, it is necessary to use statistical mechanics combined with simulation to structure the problem into hierarchies of descriptions at different levels of detail. The more finegrained the description, the shorter the time over which a simulation can be carried out to achieve a description of the system dynamics. More coarse-grained descriptions make it feasible to do simulations over a longer period. A promising strategy is to structure the problem into hierarchies. Then the results from the very fine-grained description can be used to provide input parameters to the more coarse-grained description, which can then be used to simulate or calculate directly a biologically important attribute, such as a conductance. This hierarchical strategy is diagrammed in Figure 1. This article provides an outline of how to plan and execute this strategy, to the end of understanding in detail the relationship between structure and function in ion channels. It should be kept in mind that this work, while well grounded in computational technique and statistical mechanical theory, is workin-progress. Even for the above-mentioned gramicidin channel, no completely successful atomistically detailed calculation has been done that accurately predicts ion conductance. Quantitatively successful calculations of water permeability have been done, however, both for permeation through a lipid bilayer (1) and for water permeation through a gramicidin channel (2).

342

AID

1046-2023/98 $25.00 Copyright q 1998 by Academic Press All rights of reproduction in any form reserved.

Methods 0589

/

6208$$$$41

03-23-98 16:31:15

metha

AP: Methods

PERMEATION AND SELECTIVITY IN ION CHANNELS

I. APPLICABLE METHODS a. Electrodiffusion Theory In this section we quickly outline particular aspects of the electrodiffusion theory that are critical in the application of this theory to understanding ion permeation in channels. The classic connection between electrophysiologically measured currents and physical theory is the electrodiffusion, or Nernst–Planck, equation,

343

the diffusion coefficient for ionic species i, Ci is the concentration of ionic species i, zi is the valence of ionic species i, F is Faraday’s constant, R is the gas constant, T is absolute temperature, and V is voltage. The integral form of the Nernst–Planck equation for ion permeation, which gives the ionic current as a function of the applied voltage, is Ii Å ziF

Cout,i 0 Cin,iexp(ziFDV/(RT))

*

d

,

[2]

dx [exp(ziFV/(RT))/Di]

0

Fluxi Å 0Di(0dCi/dx / (ziCiF/RT)dV/dx),

[1]

where Fluxi is the flux density of ionic species i, Di is

where Ii is the current density carried by species i, Cout,i is the extracellular concentration of species i, Cin,i is the intracellular concentration of species i, DV is the transmembrane potential difference, d is the membrane thickness or ion channel length, and the coordinate system is chosen such as to place the extracellular membrane face at x Å 0. Setting the current to zero in the Nernst–Planck equation yields the Nernst equation for the electromotive driving force for a particular ionic concentration gradient: Vi Å

S D

Cout,i RT ln . ziF Cin,i

[3]

The precise measurement of specific ionic currents and voltages was pioneered by Hodgkin and Huxley (3). However, the currents measured were not through individual ion channels but through populations of channels in excitable membranes. Early on, it was noticed that for many types of channels, the current-versus-voltage curves for open channels were linear, giving rise to the use of the conductance equation to correlate electrophysiological data: Ii Å gi(Vi 0 DV).

[4]

Here, gi is the conductance of the membrane for species i. The conductance equation, which is equivalent to Ohm’s law, is a solution to the Nernst–Planck equation for particular boundary conditions. Specifically, if one postulates that Di is constant within the membrane channels, and imposes the boundary condition on Eq. [2] that the current Ii must be a linear function of the transmembrane voltage DV, then the integral form of the Nernst–Planck equation becomes (4) Ii Å FIG. 1. Hierarchy of calculational methods and theories that can be used to correlate measured ion fluxes through channels with channel structure.

AID

Methods 0589

/

6208$$$$42

z2i F Cout,i 0 Cin,i rPir (Vi 0 DV), RT ln(Cout,i/Cin,i)

[5]

where Pi, the permeability of the membrane to species

03-23-98 16:31:15

metha

AP: Methods

ERIC JAKOBSSON

344

i, is defined as the ratio of the diffusion coefficient to the membrane thickness, i.e., Pi Å Di/d. Comparing Eqs. [4] and [5], it is seen that the relationship between the conductance and the permeability for a channel that gives linear I–V relationships is gi/Pi Å

z2i F Cout,i 0 Cin,i r . RT ln(Cout,i/Cin,i)

[6]

When the I–V relationship through an ion channel is curved, it is often fit with the constant-field equation, which is derived by solving Eq. [2] with the assumption that the voltage profile within the channel is linear: Ii Å z2i FPiDV

Cout,i 0 Cin,iexp(ziFDV/(RT)) . exp(ziFDV/(RT)) 0 1

permeation is the electrostatic interaction between the ion and the channel, the permeability ratios as determined experimentally and analyzed by Eq. [8] should equal those as determined by Eq. [9]. In fact they often do not (6). This discrepancy is termed the ‘‘anomalous mole-fraction effect.’’ The anomalous mole-fraction effect has been been taken to imply that ions interact with each other inside channels (7). Beyond that, the mole-fraction effect suggests the necessity of generalizing the Nernst–Planck equation to account for interactions beyond electrostatics. The second term on the right-hand side of Eq. [1], the term with the voltage gradient, should be changed to a free energy gradient. Therefore the Nernst–Planck equation should be rewritten as

[7]

S

Fluxi Å 0Di

Levitt (5) showed how to extend the electrodiffusion theory (originally considered by channel biophysicists to be a bulk continuum theory) to the case of a single ion channel occupied by only one ion at a time. Consider the situation of an ion channel that is partially selective between two ions of the same valence, for example, univalent cations such as sodium and potassium. Independently of any particular assumption about the potential profile within the channel (as in deriving Eq. [7]) or about constants of integration (as in deriving Eq. [5]), the Nernst–Planck equation leads to two different ways of calculating the selective permeability. One way is to measure the reversal potential, i.e., the potential at which the total current through the channel is zero. For ions of the same valence, the denominator on the right-hand side of Eq. [2] will be the same. Summing the fluxes of two univalent cations to zero as described by Eq. [2], using the identity of the denominator of the right-hand side, yields

D

dCi Ci d(DG) / r , dx RT dx

[10]

where the generalized force, d(DG)/dx, is the change in free energy of the system as the ion is moved through the channel. The integrated form of Eq. [10], analogous to Eq. [2], is Ii Å ziF

Cout,i 0 Cin,iexp(ziFDV/(RT))

*

d

.

[11]

dx[exp(DG/(RT))/Di]

0

[9]

Note that the numerator on the right-hand side of Eq. [11] is the same as on the right-hand side of Eq. [2]. This is because the relative free energy of ions in two electrolyte solutions of similar strength is dominated by the electrostatic component. Thus, Eq. [11] preserves in its original form the Nernst equation for the reversal potential when there is only a single permeant ion (Eq. [3]). But the denominator on the right-hand side of Eq. [11] now includes not only electrostatic interactions between the channel and the ion but also ion–ion interactions within the channel, ion–channel– water interactions within the channel and as the ion enters and leaves the channel, and averaging over the various conformational changes the channel undergoes as the ion translocates the channel. The entire problem of understanding the relationship between ion channel structure and ion permeation can be summarized as the problem of evaluating the integral in the denominator of the right-hand side of Eq. [11] for a given channel structure. An appropriate name for the denominator on the right-hand side of Eq. [11] is the diffusive resistance.

where it is understood that the currents for 1 and 2 are measured at the same concentrations and voltages. If the determining factor for the energy profile of ion

b. Stochastic Dynamics Based on the homology between solutions to the diffusion equations and the distribution of particles undergoing a random walk, Einstein showed that the mi-

P1 Cin,2exp(FDVr/(RT)) 0 Cout,2 Å , P2 Cin,1exp(FDVr/(RT)) 0 Cout,1

[8]

where DVr is the measured reversal potential at specified concentrations of ions 1 and 2. According to the Nernst–Planck equation, there is a second way to determine the relative permeability, and that is to measure currents carried by each of the ions in the absence of the other. In this case, the ratios of the permeabilities should equal the ratios of the currents, P1 I1 Å , P 2 I2

AID

Methods 0589

/

6208$$$$42

03-23-98 16:31:15

metha

AP: Methods

PERMEATION AND SELECTIVITY IN ION CHANNELS

croscopic basis for diffusion is random walk (8) or Brownian motion. This leads to a rewriting of the differential form of the Nernst–Planck equation into a Brownian dynamics algorithm for moving ions through the channel,

Dx Å 0

Di d(DG) r rDt / grnr(2DiDt)1/2, RT dx

[12]

where Dx is the distance the ion moves in a time increment Dt, and grn is a random number chosen from a Gaussian distribution such that the mean is 0 and the mean of the squares is 1. Depending on how rough the potential surface is, it may be necessary to use a finergrained version of Brownian dynamics, in which the particle velocity is continually updated by the random number generator,

S

kbT 1 d(DG) rv 0 r m m dx

Dv Å Dt 0

/ grnr

S D 2Dt Di

1/2

kbT , m

r

D [12*]

where Dv is the change in the ion velocity in time Dt, kb is Boltzmann’s constant, and m is the mass of the ion. For a detailed description of the use of Eq. [12*] in the context ion motion in channels, see Chiu et al. (9). In either Eq. [12] or [12*], the motion of the ions is not completely deterministic. Rather, the concentration gradient term in the original Nernst–Planck equation is replaced by a random number generator, whose output is scaled according to the diffusion coefficient, or mobility. The current is calculated by running the program for a large number of time steps and counting how many ions cross the channel each way. This requires an algorithm to have ions enter the ends of the channel in a fashion that corresponds to the electrolyte concentration. The details of this algorithm are described by Bek and Jakobsson (10). For the special case where only one ion can enter a channel at a time, it has been shown that the numerical solutions of Eq. [12] (Brownian dynamics) and Eq. [10] (continuum electrodiffusion) agree with each other (11). In this instance, it is obviously more efficient to use the continuum equation to achieve the number relating the flux to the free energy profile. An experimental example of this situation is sodium permeation through the gramicidin channel (12). It has also been shown that if the ion–ion interactions within the channel are reduced to zero (possible in the simulation but obviously not in a real multiion channel), the Brownian dynamics simulation with a linear potential profile within the membrane reproduces the constant field solution to the Nernst–Planck equation (13). A problem

AID

Methods 0589

/

6208$$$$43

345

arises when ion–ion interactions within the channel are introduced. In this case the logical way to introduce these interactions into the Nernst–Planck equation is to solve it simultaneously with Poisson’s equation, Ç2u Å 0r/e (14). This was first done by deLevie et al. (15). When the continuum calculations were compared with the corresponding Brownian dynamics simulations, it was found that the currents were markedly different, as were the statistics of the ion distributions within the channel. The continuum calculations had systematically calculated ion–ion repulsions that were too large and therefore had ion concentrations in the channel that were too small (13). The problem is with the manner of doing the accounting for the ion–ion interactions in the multiply occupied channels. Consider the case of a channel that may have any number of ions up to n. When there is one ion in the channel, it sees only the potential due to the channel. When there are two, they see the potential due to the channel and due to each other. So there is a free-energy profile associated with double occupancy that is different from that for single occupancy. This issue is iterated up to n-ion occupancy. The net result is that there is not just one free energy profile, but n 0 1 free energy profiles. In this situation, an exact calculation of the diffusive resistance for a continuum formalism requires the evaluation of not just a single integral but n-multiple integrals. Levitt (16) has considered this issue and suggested procudures for a reasonably efficient approximate numerical solution of this problem. While computationally more intense than Levitt’s approximate solution of the multiion electrodiffusion problem, the Brownian dynamics computations are straightforward and unambiguously correct. In comparing the Brownian dynamics method with continuum methods for computing fluxes in multiply occupied channels, a critical point is that the computing time for the Brownian dynamics goes up roughly linearly with the number of ions in the channel. The computing time for evaluating nested integrals such as are necessary for an exact solution of the continuum equations would go up roughly exponentially with the number of ions in the channel. Thus for multiply occupied channels, we would reasonably expect it to be both simpler and also more computationally efficient to use Brownian dynamics to compute the fluxes, rather than to attempt an exact continuum solution for the flux equations. With respect to the domain of applicability of different methods, it should be noted that the single-occupancy situation can be made to pertain at low permeant ion concentrations. For example, while it is clear that at least two cesium ions can get into a gramicidin channel, at low cesium concentrations one sees characteristic single-occupancy behavior. Tracer-labeled flux ratio studies do not show evidence for multiple occupancy in voltage-gated sodium channels under physiological

03-23-98 16:31:15

metha

AP: Methods

ERIC JAKOBSSON

346

conditions (17). On the other hand, sodium channels do exhibit a dependence of PNa/PK on ionic concentrations, which is commonly taken as evidence of multiple occupancy [see review of Begenisch (18)]. Most potassium channels are clearly multiply occupied, by all experimental criteria, with the exception of sarcoplasmic reticulum K channels [see review of Yellen (19)]. c. Molecular Dynamics The finest-grained level of description feasible for describing biomolecular dynamics is as a set of charged balls (the atoms), connected by springs (the bonds). The charges on the atoms, the polarizability of the atoms, the sizes of the atoms, and the energetics of deforming the bonds are collectively known as the force field of the system. Force fields for biomolecular simulations are derived partly from spectroscopic studies on small molecules that are essentially parts of macromolecules, partly from quantum-mechanical calculations, and partly are adjusted to fit known constitutive behavior, such as self-diffusion coefficient and dielectric constant. Computing molecular dynamics trajectories is conceptually simple, since atoms simply move according to Newton’s laws, but is computationally quite complex, because of the large number of atoms in biomolecular systems and near environment and the various types of forces involved. Several programs are conveniently organized to simulate biomolecular systems. Three of the most commonly used are AMBER, GROMOS, and CHARMM. Given a structural molecular model of a biomolecule, the essence of a molecular dynamics simulation is to assemble the model molecule (i.e., assign spatial coordinates to each atom in the molecule) together with as much of the environment as is needed to provide realistic behavior, initiate thermal motion by assigning random velocities to the atoms with a Maxwellian velocity distribution characteristic of the temperature at which the simulation is to be run, and then let the system move according to Newton’s laws. By collecting statistics on the motions, conformations, and the energetics of the system as it moves, one attempts to gain insights into the relationship between the molecular structure and its function. A problem with molecular dynamics is that it is only feasible for simulations of short durations. Even for a very reduced model of a small channel (e.g., a gramicidin channel without explicit surrounding membrane), a simulation of a nanosecond duration will take months. For much larger channels, or for a small channel with a more completely modeled environment, simulations would take correspondingly longer. As indicated in the Introduction, we must use molecular dynamics in some more clever way than simply simulating the system long enough to let lots of ions pass through it.

AID

Methods 0589

/

6208$$$$43

It is possible to use sampling techniques and statistical mechanics to infer mobilities and free energy profiles from molecular dynamics. These parameters can then be inserted into more coarse-grained theories, such as stochastic dynamics and the Nernst–Planck equation. A key element to add to the molecular dynamics is importance sampling. Importance sampling gets around a key problem in molecular dynamics. We expect in general that a free energy profile along the ion permeation pathway consists of a series of potential wells and barriers. For example, it would be expected that negatively charged amino acid residues would provide potential wells for cations in the permeation pathway. In the gramicidin channel, periodically disposed carbonyl oxygens from the protein backbone provide the potential wells. The rate of permeation is defined by how frequently the ions cross the free energy barriers, but the natural tendency of the ‘‘raw’’ molecular dynamics trajectories is to spend practically all their time in the potential wells. Thus important rate-limiting regions of the free-energy profile are not sampled. In importance sampling, the ion is artificially restrained in an added harmonic potential arranged so as to force the ion to spend equal time near each position along the permeation pathway. Then in each small region in which the ion is restrained, the shape of the free energy curve is derived from applying Boltzmann statistics to the distribution of ionic positions in that region, according to the equation

S D

DGx/d px Å RT ln , px/d DG x

[13]

where the subscripts x and x / d represent two positions very close to each other along the ion permeation pathway and px and px/d represent probabilities of finding the ion at those two positions as computed by the molecular dynamics, following a correction for the artificial restraints that hold the ion in the vicinity of x. [This procedure and analysis are described in detail in the context of ion permeation by Roux and Karplus (20).] The entire free energy curve is pieced together from the shapes of the curve in each region. The mobility of the ion in each region of the channel can be ascertained by analyzing the magnitude of the short-term fluctuaions in its trajectory. From these fluctuations one computes the mean-square-deviation (MSD) function:

K

L

MSD(t) Å (x(t / t) 0 x(t))2 .

[14]

To evaluate the MSD, one compiles, from the fluctuating positions of the ion in the molecular dynamics calculations, the mean distance that the ion moves along

03-23-98 16:31:15

metha

AP: Methods

PERMEATION AND SELECTIVITY IN ION CHANNELS

the permeation pathway in a time interval t. The diffusion coefficient in the region is just twice the slope of the computed MSD-versus-time lag: d(MSD) . dt

Di Å 2r

[15]

Methods 0589

/

6208$$$$44

2. APPLICATION TO KNOWN OR POSTULATED CHANNEL STRUCTURES a. Gramicidin: A Known Structure i. Water Permeability

For a detailed description of the application of Eqs. [14] and [15] to molecular dynamics simulation data from an ion channel, see Chiu et al. (9). For the special case in which each ion moves independently of each other ion (single occupancy) the diffusive resistance (Eq. [11]) may be calculated from the free energy curve (calculated with Eq. [14] from MD simulations) and the diffusion coefficient along the permeation pathway (calculated with Eq. [15] from molecular dynamics simulations). In principle this solves the problem of computing permeation from structure, with a couple of significant caveats. These are (1) that the force fields and other aspects of the methodologies that go into the molecular dynamics simulations are accurate enough to give good free energy and mobility profiles to match the data, and (2) that the importance sampling is efficient enough to give us good sampling of the conformational states of the system in a reasonable computation time. We return to these issues later in the article when we review the results of actual computations that have used these methods. In the event that the channel can be multiply occupied, the theoretical and methodological issues become more complex. Now, in place of the analytical solution to Eq. [11], one must turn to stochastic dynamics simulation as a solution. Beyond that, one has the task of calculating the ion–ion interactions in the channel. Because of the dielectric heterogeneity of the channel– water–lipid environment, the ion–ion interactions will not be well described by any simple electrostatics (21). One could postulate a dielectric constant intermediate between the aqueous and lipid values, as has been done in calculations in which the goal is to determine the qualitative phenomena that arise from ion–ion interactions in the channel (10). But that is not likely to be accurate enough for determining structure–function relationships. The preferred procedure would be to solve the Poisson–Boltzmann theory for the induced potential at all positions in the channel due to an ion at any position in the channel, and then put the results into the stochastic dynamics program as a lookup table for ion–ion forces at any instant, depending on the ionic positions. Recently, fast algorithms have been developed for solving the full nonlinear Poisson–Boltzmann equation (22), and these algorithms have been incorporated into the simulation package UHBD. Thus, although these calculations are numerically complicated, they have been well automated so that they are readily accessible.

AID

347

The structure of the gramicidin channel is known to a high resolution from solid-state NMR studies (23, 24). Therefore it is a good system to test the methodologies of structure–function prediction. Of all known ion channels, it is the one that has the highest amount of obligatory water movement associated with the ion motion, due to its long narrow pore shape (25). A variety of simulation studies have been done on water motion in the channel, based on molecular dynamics computations with explicit lipid removed from the system, to save computer time [see review by Roux and Karplus (26)]. All studies agree that the water motion is obligatorily single file and that the water mobility in the channel is lower than in bulk, but higher than the mobility of ions in the channel. (Another way of stating the relative mobility of ions and water in the channel is to state that the mobility of a chain of water molecules in the gramicidin pore is dramatically reduced by the presence of one or more ions.) In the studies in our laboratory, we exploited the obligatory single filing of the waters to analyze the fluctuations of the center-ofmass of the water chain in the pore, as if this chain were a single particle. From the MSD analysis (Eq. [15]), we showed that the local friction governing diffusive motion of the waters in the channel was about twice the diffusive friction in bulk water (9). Furthermore, the chain of waters in the channel was seen to be moving in a potential well of height about 3.6 kT and ˚ , which caused the effective diffusion width about 2.9 A coefficient for water to translocate the channel to be several times smaller than the diffusion coefficient de˚ width of the potenfined by the local friction. The 2.9-A tial barrier is about size of the water–water spacing in the channel. This correspondence leads to the hypothesis that the source of the barrier is the energetic cost of a water at the end of the chain entering/leaving one end of the channel and a water at the other end leaving/entering the end of the channel. Another significant finding is the extent to which the water motion depended on the protein motion. In the simulations one can do an experiment that is impossible in ‘‘real life’’: one can hold the protein motionless, effectively at 0 K, and keep the water at 300 K. In this computer experiment, the effective diffusion coefficient for water translocation drops to a very low value (27). Effectively, the barrier for water to enter and leave the channel became very high. In another computer experiment, artificial restraints of different strength were placed on the side-chain motions of the channel,

03-23-98 16:31:15

metha

AP: Methods

ERIC JAKOBSSON

348

while leaving the backbone (that forms the channel lining in gramicidin) unrestrained (28). It was seen that the side-chain restraints had a major effect on the magnitude of the effective diffusion coefficient for water translocation through the channel. Since the side chains are interacting with the lipid environment, it became important to do simulations with the gramicidin embedded in a lipid bilayer, to understand the water permeation in this system. We have done simulations of water transport through the gramicidin channel in an environment of excess lipid (dimyristoyl phosphatidylcholine) and excess water (2). Our results are based on a 1.2-ns simulation. Although this is not a sufficiently long time for individual water molecules to translocate the entire channel, it is sufficient time for five water molecules to make complete transitions between the channel environment and the bulk aqueous environment. By considering the motion of the waters in the channel to be essentially a ‘‘shaking stack’’ (29) and inferring an effective diffusion coefficient for the stack by the frequency with which waters come on and off the ends of the stack, we calculated that the five events in 1.2 ns corresponded to an effective diffusion coefficient for water in the channel of 3.5 1 1007 cm2/s. The measured effective diffusion coefficient for water-permeating gramicidin channels in a phosphatidylcholine (PC) membrane is 2.8 1 1007 cm2/s (30). Thus the calculation is quantitatively successful, suggesting that the force fields and other aspects of the molecular dynamics methodology are accurate enough for computing water permeabilities. By computing the MSD correlation for the waters within the channel we were able to compute the fraction of the system diffusive resistance that came from the inside of the channel, as opposed to the region just outside the channel mouth, where the water had to move through a space impinged on by phospholipid head groups. We found that only about one-sixth of the total diffusive resistance came from the channel itself, with the rest coming from the interactions of the water just outside the channel with the phospholipid head groups. This result from the simulations seems to explain the discrepancy between the above-cited experimental water permeability for gramicidin in PC membrane and a measurement six times larger in glycerol monolein membrane (31). Presumably the difference lies in the relative lack of head groups that restrict access to the mouth of the channel. ii. Ion Permeability Protein motions are as critical for ion permeation as for water permeation. Molecular dynamics simulations with an ion near the mouth of the channel reveal a deformation of the b-helix that apparently provides a low-energy intermediate for moving the ion from the bulk-solvated environment to the channel-solvated en-

AID

Methods 0589

/

6208$$$$44

vironment without the need for crossing a significant energy barrier (32). Roux and Karplus (20) have used the importance sampling technique combined with molecular dynamics as described above to construct a potential of mean force profile for sodium permeation in gramicidin. For this singly occupied situation, it has been possible to use electrophysiology theory to assess the extent to which the calculated potential of mean force corresponds to the actual situation (8, 33, 34). The result of this analysis is that the shape of the potential curve is about right but the depths of the potential wells are too large, so the computed permeability would be far too low. It seems likely that a major part of the problem is the inadequacy of the standard force fields for proteins and water to deal with the large electric fields in the vicinity of such a small ion as sodium. There needs to be a way to include polarizability of both water and protein groups near the ion. Another indication of this need is the comparison between computational and NMR experimental studies of sodium binding in gramicidin (35, 36). In this instance there is close agreement between the simulation and experimental studies as to the location of the binding site but some disagreement as to the extent to which the channel protein is distorted by the ion binding. b. Model-Built Channels i. Water Permeability Computational studies have been done on solvation and water mobility in channel-like featureless cavities (37) and in a model of the pore region of the voltagegated sodium channel [model constructed by Guy and Durell (38), solvation studies by Singh et al. (28)]. For most of their length the modeled sodium channel is wider than the gramicidin channel, and is lined with amino acid side chains, rather than the protein backbone as is the case with gramicidin. Notable features of the sodium channel simulations include the following observations: The water molecules tend to be oriented within the channels according to the intrinsic polarization of the channel. Even for channels wide enough to accommodate a few water molecules across, the mobility of the water within the channel is substantially reduced from bulk value. In the model sodium channel simulations (28), it was also shown that the water could be effectively completely immobilized by taking the protein temperature close to absolute zero, even if the water remained at 300 K. The water continued to make very small thermal fluctuations, but fluctuations of a magnitude necessary for translational diffusional motion are prevented by immobilizing the protein. ii. Ionic Permeability Studies have begun on the ionic permeation of the sodium channel model described above. The early results are as follows.

03-23-98 16:31:15

metha

AP: Methods

PERMEATION AND SELECTIVITY IN ION CHANNELS

Prior to doing molecular dynamics, a study was done in which an ion was moved through the permeation pathway and the system was energy minimized with the ion at different locations (39). From the comparative system energies with different ionic species at different locations, it was clear that the model channel would select for cations over anions, but there was no clear basis as revealed by this method for sodium over potassium selectivity. An interesting result was seen for one particular residue, a lysine situated at the narrowest region of the channel according to the model. Although it may seem odd for a positive residue to be situated in the narrow part of a cation channel, experimental evidence suggests that this residue is critical in the determination of the channel selectivity. When this residue is mutated to a negatively charged glutamate, the sodium channel becomes calcium selective (40). When the lysine is mutated to a positively charged arginine, the channel selects for univalent cations, but does not select well between sodium and potassium (41). In the energy minimization calculations, this lysine is seen to undergo a distinctive side-chain transition as a sodium ion is moved through the narrow part of the channel, effectively changing the channel diameter by approximately twofold as a function of the ion position. We presume that specific interactions between the ion, water molecules, and the protein in this region, especially the lysine, are critical for determining the channel selectivity. However, the energy minimization calculations do not sample conformation space well enough to give a good understanding of the interaction. The next step was to do importance sampling combined with molecular dynamics, in a fashion similar to the work previously done on the gramicidin channel by Roux and Karplus (20). These calculations are currently being analyzed. Because of side-chain conformation transitions as the ion permeates, the free energy profile is much more complex than for permeation of gramicidin. It is not yet clear whether the same degree of sampling that is adequate for gramicidin will suffice to give an accurate free energy profile in the Na channel.

III. SUMMARY a. Accomplishments to Date There has been developed a well-defined body of theory and calculational methods that give promise of being able to compute permeation rates and selectivity from ion channel structure. These methods include a combination of molecular dynamics, statistical mechanical theory, electrostatics, stochastic dynamics, and the Nernst–Planck equation with the electrical

AID

Methods 0589

/

6208$$$$45

349

potential profile generalized to a free energy profile. These methods will be linked in a hierarchy in which the description that explicitly includes experimentally determined quantities (Nernst–Planck) will be parameterized by stochastic dynamics results, which in turn will be parameterized by molecular dynamics results. So far the complete connection between the molecular dynamics and an experimentally measured permeability is quantitatively shown only for water permeation of lipid and gramicidin. Sodium permeation in gramicidin is well understood but the numbers computed for conductance directly from molecular dynamics simulations give values too low for conductance, probably because of the inaccuracies that arise in molecular dynamics when proteins and water are subjected to the large fields such as are generated very near such a small ion as sodium. For more complicated situations it reasonable to believe that our theoretical and computational methods are powerful enough to relate channel structure to conductance and selectivity, but the details of conformational sampling in the channel are not yet proven to work. An important general finding from the work to date is that the motions of the channel protein are critical to channel function. The classical view that channels, in contrast to carriers, are fixed structures, is not correct. The motions of channel pores and selectivity filters, while undoubtedly of smaller amplitude than the conformational changes of carriers, are no less vital to channel functioning. b. The Next Challenges The first immediate challenge is to improve the force fields for the interactions among ions, water, and proteins. The likely avenue of improvement is to include polarizability in the protein and water force fields, so that the energies of the three-way ion–water–protein interaction can be calculated accurately enough to predict experimental conductance and selectivity directly from molecular dynamics computations. One good model system for this effort will continue to be sodium permeation of gramicidin. This system has a known structure, has the relative simplicities of single occupancy and backbone lining the pore lumen, and tests the force fields severely because of the large fields around the small sodium ion. The second immediate challenge is to develop and demonstrate the adequacy of sufficiently powerful sampling methods so that protein conformations can be adequately sampled to provide an accurate free energy profile for ion permeation. In addition to the importance sampling combined with molecular dynamics, described above, other methods of conformational sampling for the protein should be explored. Recent improvements in Monte-Carlo sampling techniques should be explored for this purpose (42). Model-built structures for the sodium channel are reasonable sub-

03-23-98 16:31:15

metha

AP: Methods

ERIC JAKOBSSON

350

jects for this type of study (37, 43). Although the structure is not precisely known, there are enough experimental constraints on possible structures that reasonable models can be made, and the channel has the simplifying virtue of less extreme multiion interactions than K channels, so that treatment as a single-ion channel may be a reasonable first approximation. For multiply occupied channels, the sampling problem is compounded by the need to sample over distributions of ions within the channel. Gramicidin permeated by alkali metals such as potassium, cesium, and rubidium is a good system for such studies, because the electrophysiology clearly indicates multiple occupancy for these ions (44). However, the gramicidin structure is not typical of other ion channels, in that the pore lumen is lined by backbone rather than side chains. Another potentially good class of systems for testing sampling methods is porins. There are known high-resolution structures that can provide the basis for realistic molecular dynamics studies (45). This class of channels has the complexity of side chains lining the pore and, based on the geometry, would be expected to be multiply occupied. c. Possibilities for the Future It seems likely that it will be feasible to correlate channel permeability and selectivity from channel structure, using refinements and further extensions of the methods developed so far. As the simulation methods and associated theory described in this article improve in parallel with improvements in membrane protein structure determination and structure prediction, it is hoped that the combination of these methods will ultimately make it feasible to predict channel properties from the channel protein sequences.

ACKNOWLEDGMENTS I am grateful to many colleagues over the years for helpful discussions and collaborative interactions. Professor Shankar Subramaniam made helpful comments on an earlier draft of this paper. Work from our laboratory cited in this paper was supported at different times by the National Institutes of Health, the National Science Foundation, the American Heart Association, and the Pfizer Pharmaceutical Company. Computer time was provided by a grant from the National Center for Supercomputing Applications.

REFERENCES 1. Marrink, S. J., and Berendsen, H. J. C. (1994) J. Phys. Chem. 98, 4155–4168. 2. Jakobsson, E. G., Chiu, S. W., and Subramaniam, S. (1998) Biophys. J. 74, A232. 3. Hodgkin, A. L., and Huxley, A. F. (1952) J. Physiol. 116, 449– 472.

AID

Methods 0589

/

6208$$$$46

4. Jakobsson, E., and Jaslove, S. (1981) Biophys. J. 29, 129A. 5. Levitt, D. G. (1986) Annu. Rev. Biophys. Bioeng. 15, 29–57. 6. Heginbotham, L., and MacKinnon, R. (1993) Biophys. J. 65, 2089–2096. 7. Hille, B. (1992) Ionic Channels of Excitable Membranes, 2nd ed., Sinauer, Sunderland, MA. 8. Einstein, A. (1956) Investigations on the Theory of the Brownian Movement, Dover, New York (republication in translation of original paper of 1908). 9. Chiu, S.-W., Novotny, J. A., and Jakobsson, E. (1993) Biophys. J. 64, 98–108. 10. Bek, S., and Jakobsson, E. (1994) Biophys. J. 66, 1028–1038. 11. Jakobsson, E., and Chiu, S.-W. (1987) Biophys. J. 52, 33–45. 12. Procopio, A., and Andersen, O. (1979) Biophys. J. 25, 8a. 13. Cooper, K., Jakobsson, E., and Wolynes, P. (1985) Prog. Biophys. Mol. Biol. 46, 51–96. 14. Chen, D., and Eisenberg, R. (1993) Biophys. J. 64, 1405–1421. 15. De Levie, R., Seidah, N. G., and Moreira, H. (1971) J. Membr. Biol. 10, 171. 16. Levitt, D. G. (1991) Biophys. J. 59, 271–277. 17. Hodgkin, A. L., and Keynes (1955) J. Physiol. 128, 61–88. 18. Begenisich, T. (1987) Annu. Rev. Biophys. Biopphys. Chem. 16, 247–263. 19. Yellen, G. (1987) Annu. Rev. Biophys. Biophys. Chem. 16, 227– 246. 20. Roux, B., and Karplus, M. (1991) Biophys. J. 59, 961–981. 21. Jordan, P. (1982) Biophys. J. 39, 157–164. 22. Holst, M. J., Kozack, R. E., Saied, F., and Subramaniam, S. (1994) Proteins Struct. Funct. Genet. 18, 231–245. 23. Ketchem, R. R., Hu, W., and Cross, T. A. (1993) Science 261, 1457–1460. 24. Koeppe, R. E., II, Killian, J. A., and Greathouse, D. V. (1994) Biophys. J. 66, 14–24. 25. Finkelstein, A. (1987) Water Movement through Lipid Bilayers, Pores, and Plasma Membranes: Theory and Reality, Wiley–Interscience, New York. 26. Roux, B., and Karplus, M. (1994) Annu. Rev. Biophys. Biomol. Struct. 23, 731–761. 27. Chiu, S.-W., Jakobsson, E., Subramaniam, S., and McCammon, J. A. (1991) Biophys. J. 60, 273–285. 28. Singh, C., Sankararamakrishnan, R., Subramaniam, S., and Jakobsson, E. (1996) Biophys. J. 71, 2276–2288. 29. Schumaker, M. F. (1992) Biophys. J. 63, 1032–1044. 30. Rosenberg, P. A., and Finkelstein, A. (1978) J. Gen. Physiol. 72, 341–350. 31. Dani, J. A., and Levitt, D. G. (1981) Biophys. J. 35, 501 – 508. 32. Bobak, M., Chiu, S. W., and E. Jakobsson (1991) J. Mol. Graphics 9, 44–45. 33. Chiu, S-W., and Jakobsson, E. (1989) Biophys. J. 55, 147 – 157. 34. McGill, P., and Schumaker, M. F. (1996) Biophys. J. 71, 1723– 1742. 35. Woolf, T. B., and Roux, B. (1997) Biophys. J. 72, 1930–1945. 36. Tian, F. K., Lee, C., Hu, W., and Cross, T. A. (1996) Biochemistry 35, 11959–11966. 37. Sansom, M. S. P., Kerr, I. D., Breed, J., and Sankararakrishnan, R. (1996) Biophys. J. 70, 693–702. 38. Guy, H. R., and Durell, S. R. (1995) in Ion Channels and Genetic

03-23-98 16:31:15

metha

AP: Methods

PERMEATION AND SELECTIVITY IN ION CHANNELS Diseases (Dawson, D. C., and Frizzell, R. A., Eds.), Rockefeller Univ. Press, New York. 39. Sankararakrishnan, R., Subramaniam, S., and Jakobsson, E. (1997) Biophys. J. 72, A360. 40. Heinemann, S. H., Terlau, H., Stu¨hmer, Imoto, K., and Numa, S. (1992) Nature 356, 441–443. 41. Favre, I., Moczydlowski, E., and Schild, L. (1996) Biophys. J. 71, 3110–3125.

AID

Methods 0589

/

6208$$$$46

351

42. Scott, L. R. (1996) in Biological Membranes: A Molecular Perspective from Computation and Experiment (Merz, K. M., and Roux, B., Eds.), pp. 83–104, Birkha¨user, Boston. 43. Lipkind, G. M., and Fozzard, H. A. (1994) Biophys. J. 66, 1–13. 44. Finkelstein, A., and Andersen, O. S. (1981) J. Membr. Biol. 59, 155–171. 45. Watanabe, M., Rosenbusch, J., Schirmer, T., and Karplus, M. (1997) Biophys. J. 72, 2094–2102.

03-23-98 16:31:15

metha

AP: Methods