INTRODUCTION

TOP ABSTRACT INTRODUCTION PNP THEORY COMPUTATIONAL IMPLEMENTATION PNP LATTICE MODEL OF... INFLUENCE OF MEMBRANE CHARGES... COMPARISONS WITH EXPERIMENTAL... CONCLUSIONS APPENDIX REFERENCES

Over the past decade methodological developments in NMR, x-ray crystallography, and electron spectroscopy have led to significant progress in determining structures of integral membrane proteins that form ion channels (Doyle et al., 1998; Song et al., 1996; Ketchem et al., 1997). This accumulation of high-resolution structural information has enabled better understanding of channel conductance, gating, and selectivity. Phenomenological descriptions of ion flow through channel proteins in terms of kinetic, electrodiffusion, and stochastic models have been given (Cooper et al., 1985). At a more atomistic level, equilibrium molecular dynamics simulations (Roux and Karplus, 1993) and nonequilibrium Brownian dynamics simulations (Chung et al., 1998, 1999; Corry et al., 1999) show long-term promise for elucidating detailed mechanisms of biological ion transport processes (Levitt, 1999).

Interest in biological ion channels is stimulated by the important role they play in regulating the electrical properties of cells, as well as their tendency to bind various antibiotic and toxin molecules (which can drastically alter their functional properties, particularly current flow through the channel). One of the most widely studied ion channels is the neutral pentadecapeptide gramicidin A (GA), which forms aqueous pores in lipid bilayers that selectively pass monovalent cations (Andersen, 1984; Koeppe and Andersen, 1996). Some of the functional properties of GA, such as single-file flux and ion valence selectivity, are found in physiologically important channels. This has made the GA channel a focal point of many theoretical and experimental studies.

Recent experiments on the gramicidin A channel have elucidated the influence of membrane electrostatics on its ion conductance. The dependence of GA channel conductance on lipid bilayer surface charge through the phosphatidylserine (PS) lipid was studied recently by Rostovtseva and co-workers (Rostovtseva et al., 1998). Two methods of varying surface charge were utilized, namely 1) titration of the surface charge by changing the pH of the bulk solution, and 2) reduction of lipid surface charge density by diluting the PS lipid with the neutral phosphatidylcholine (PC) lipid. The negatively charged PS membrane attracts cations to the surface of the membrane (and hence to the mouth of the channel), leading to increased cation conductance.

Busath and co-workers (Busath et al., 1998) studied the effect of electrostatic middle-range and long-range interactions on GA conductance, i.e., interactions with distant residues and noncontact pore waters, and with lipid molecules and bulk water, respectively. They measured single-channel current-voltage relationships for GA in two different planar bilayers, one being the dipolar PC membrane also utilized by Rostovtseva et al. (1998) and the other the glycerylmonoolein (GMO) bilayer, which is less dipolar than PC lipid. The experimental differences in ion conductance were explained in terms of a difference in interfacial dipole potentials for the two membranes.

In this paper we utilize Poisson-Nernst-Planck (PNP) theory to study the influence of membrane surface charge density and interfacial dipole potentials on the conductance of GA channel embedded in lipid bilayers. PNP theory combines the Nernst-Planck theory of electrodiffusion with the Poisson equation of electrostatics, including contributions from both fixed charges in the system and from the mobile charge density flowing through the channel. PNP theory has previously been applied to the study of ion transport in electrochemical liquid junction systems (Riveros et al., 1989) and electron transport in semiconductor devices (Markowich, 1986), as well as ion permeation through biological membrane channels (Eisenberg, 1998; Nonner et al., 1999).

There is some disagreement concerning the applicability of the PNP approach to microscopical channels (Miller, 1999; Corry et al., 1999) when the size of the permeant ion becomes comparable to that of the pore. The equilibrium counterpart of PNP, Poisson-Boltzmann (PB) theory, has been widely utilized, with considerable success, to calculate solvation and reorganization energies of molecules in solution (Sitkoff et al., 1994; Sharp, 1998). It is obvious that application of a continuum theory like PB or PNP to a microscopic system pushes the theory beyond the scope of its derivation, and thorough testing is needed to decide if the properties of the model resemble the properties of the real system. Until recently all applications of PNP theory to ion transport in biological systems either were one-dimensional (1D) or were performed in simplified geometries such that a detailed description of the protein structure and electrostatic properties was not possible. To compare the results of the calculations of a 1D model with the corresponding experimental values, the input parameters to the model must be regarded as fitting coefficients. Therefore, for a reduced 1D model it is difficult to ascertain which properties of the real system can be accounted for by the theory, given sufficient resolution of the structural details of the 3D system, and which ones result from parameter fitting (and thus have no physical meaning).

In the present work, we employ the lattice relaxation algorithm described in a recent paper (Kurnikova et al., 1999) to solve the PNP equations for a fully 3D model of the GA channel/charged membrane system. This procedure is based on a mapping of the protein and the embedding membrane onto a 3D cubic lattice, which defines dielectric boundaries, a distribution of fixed (partial) charges associated with the protein and lipid membrane, and a flow region for the mobile ions. The Poisson and Nernst-Planck equations are then solved self-consistently until convergence is achieved. This type of algorithm has been used before to solve the Poisson and the Poisson-Boltzmann equations for biophysical systems (Nicholls et al., 1990; Luty et al., 1992).

In Kurnikova et al. (1999) the accuracy of the 3D PNP algorithm was calibrated using both parallel-plate and cylinder models (for which numerically exact results could be obtained by 1D methodology). The 3D algorithm was then applied to a GA dimer channel embedded in a neutral membrane, focusing on the influence of the atomic partial charges of the protein upon the ion flux through the channel. Reasonable agreement with experimental current/voltage (I-V) results was obtained via PNP for this narrow channel, with the diffusion coefficients for the permeant ions chosen to have reasonable physical values.

In the implementation of 3D PNP theory presented by Kurnikova et al. (1999), some additional simplifying assumptions were made, namely that the concentrations of the mobile ions retain their bulk values right up to the channel mouths, and that the external potential associated with electrodes brought into proximity to the channel mouths can be fixed as a boundary condition rather close to (within 5 Å of) the mouth openings. These simplifications were invoked for ease of numerical implementation in those exploratory calculations. Further reflection suggests that both restrictions should be removed. Experimentally, electrodes used to probe applied voltages across the channel are situated at least a micron away from the channel openings. Because of the mobility (and hence the polarizability) of the electrolyte ions in water solvent, a uniform asymptotic voltage is obtained in the water outside of the immediate vicinity of the electrodes and the channel/membrane system. The electrolyte concentrations assume their bulk values in this regime.

We have recently studied how these asymptotic boundary conditions can be converted into boundary conditions appropriate for a finite-box real-space lattice PNP calculation (Graf et al., manuscript submitted for publication). For systems like the GA dimer/membrane complex at experimentally relevant ionic strengths and applied voltages, we found that the asymptotic state is obtained at a distance of ~15 Å from the membrane surface. Thus, fixing the mobile ions to their bulk concentrations and the electric potential to the applied electrode potential on a "simulation box" boundary at about this distance from the surface of the membrane and then solving the lattice PNP equations produces an I-V curve that does not change when the box boundaries are moved farther from the channel/membrane system.

Because the electrolytes are now mobile in the region between the box boundary and the channel membrane, it is important to include the possibility that different diffusion coefficients apply in the reservoir versus inside the channel pore. In the former region, we expect that the bulk diffusion constant values apply for each ion species. However, inside the channel diffusion is significantly constrained by the narrow pore walls. Intuition and related molecular dynamics experiments (Lynden-Bell and Rasaiah, 1996; Smith and Sansom, 1998) suggest that the diffusion coefficient inside the channel region can be significantly lower than its bulk analog. Consequences of these effects have also recently been elucidated by Graf and co-workers (Graf et al., manuscript submitted for publication).

The modifications just described, that is, including mobile ions in the reservoirs outside the channel and allowing for different diffusion constants in the reservoirs versus the channel interior, were incorporated into the calculations presented in this work. They are particularly important for the systems under study here, because the charges or dipoles embedded in the membrane affect the motion of ions near the membrane.

In the section titled PNP Theory we briefly review the PNP 3D equations. Their solution using a lattice relaxation algorithm is discussed in Computational Implementation, as is the basic model for the channel system. In the next section, PNP Lattice Model of the Gramicidin A Channel, the specific 3D lattice model used to represent the GA dimer is described. Then, in the following section, Influence of Membrane Charges and Interfacial Dipoles on GA Conductance, we study the influence of surface charges and interfacial dipoles on the GA channel conductance, using simple membrane models. In Comparisons with Experimental Observations we use these models to compare our results with recent experimental data. This is followed by concluding remarks.

	PNP THEORY

TOP ABSTRACT INTRODUCTION PNP THEORY COMPUTATIONAL IMPLEMENTATION PNP LATTICE MODEL OF... INFLUENCE OF MEMBRANE CHARGES... COMPARISONS WITH EXPERIMENTAL... CONCLUSIONS APPENDIX REFERENCES

The dynamical behavior of Brownian particles in the high-friction regime is governed by the Smoluchowski equation (Chandrasekhar, 1943),

(1)

where c(, t) is the concentration of these particles at position and the flux of particles is given by

(2)

Here D() is the spatially dependent diffusion coefficient (assumed to be isotropic), V() is the external potential energy acting on the particles, and ¹ = kT, where k is Boltzmann's constant and T is the absolute temperature. The first term of the right-hand side of Eq. 2 describes the motion of particles through a concentration gradient according to Fick's first law of diffusion. The second term accounts for the drift velocity, D()V()/kT, induced on the particles by the external force V().

For steady-state conditions, i.e., c(, t)/t = 0, the Smoluchowski equation (also called the Nernst-Planck equation) can be written as

(3)

If the concentration values at the boundary surfaces are known, this equation has a unique interior solution, which can then be input into Eq. 2 to determine the particle flux vector at any point in space and hence the ion current through any surface.

When the mobile particles are charged the potential energy that appears in the Nernst-Planck equation (Eq. 2) can have electrostatic and nonelectrostatic components, i.e.,

(4)

where U() is a nonelectrostatic potential energy assumed for simplicity to be the same for all ion species, z_ie is the charge of the ion species i (z_i is its valence and e is the magnitude of the electron's charge), and () is the electric potential. The potential U() represents repulsion by fixed objects (e.g., the pore walls) or the effect of short-range repulsions between mobile ions (Levitt, 1991a,b). [In fact, we prevent mobile ions from "going through the wall" via zero-flux boundary conditions (see below), and short-range ion-ion interactions are neglected here. So, in practice, the U() term is absent from the present version of our 3D PNP algorithm.]

The electric potential depends on the charge distribution of ions in the aqueous phase, any other fixed charged species in the system, the dielectric properties of the medium, and any external electric voltage applied across the system. The electrical potential profile is determined by solving the Poisson equation,

(5)

where () is the dielectric constant profile, _f() is the density of fixed charges in the system, and the second term accounts for the charge density of the N ionic species in the solution. Because, according to Eq. 5, the electrical potential depends on the concentration of ions in solution, we must solve it self-consistently with the steady-state Nernst-Planck equations (Eq. 3) for each ionic species in the liquid phase:

(6)

where D_i() is the spatially dependent diffusion coefficient appropriate to species i. Equations 4-6 comprise the essence of Poisson-Nernst-Planck (PNP) theory.

	COMPUTATIONAL IMPLEMENTATION

TOP ABSTRACT INTRODUCTION PNP THEORY COMPUTATIONAL IMPLEMENTATION PNP LATTICE MODEL OF... INFLUENCE OF MEMBRANE CHARGES... COMPARISONS WITH EXPERIMENTAL... CONCLUSIONS APPENDIX REFERENCES

To solve the 3D PNP equations for a protein channel/membrane system, we discretize the system onto a cubic lattice grid and then solve a finite-difference representation of the PNP equations on the grid. The self-consistent solution of these equations is obtained using a Successive Over-Relaxation (SOR) algorithm (Press et al., 1992; Coalson and Beck, 1998), including zero-flux boundary conditions for lattice points next to the channel and membrane walls. For the solution of the Poisson equation we used DelPhi (Nicholls et al., 1990), a finite-difference-based computer code widely utilized in applications to biophysical systems. The DelPhi program was modified to allow inclusion of a membrane slab around the protein channel, as well as different salt concentrations on either side of the channel and an electric potential difference across it. The algorithm for solving the NP equations has been described before (Kurnikova et al., 1999) for the case of a constant diffusion coefficient. In the Appendix of the present paper we show how this algorithm can be modified to allow spatial variation of the diffusion coefficient. This procedure was used in all calculations presented in this work.

The channel-membrane system is discretized on the grid: each grid point is characterized by the concentration of each ionic species, the electrical potential, the dielectric constant, and the diffusion coefficient. Ion concentrations are nonzero only in the aqueous region. The dielectric constant  takes the value _a = 80 in the aqueous region and _m = 2 in the protein/membrane region. Within the ion flow region, the diffusion constant for species i takes the value D_b,i in the bath region and D_c,i in the channel region (specific values for the systems under consideration in this work are noted below).

The boundaries separating the solvent molecules and mobile ions from the fixed membrane/protein region are taken to be the solvent molecule- and ion-accessible van der Waals surfaces. These are obtained using the method of Connolly (1983) as implemented in DelPhi.

We used uniform cubic lattices of 131-145 grid points per side to produce the results shown in this paper. Fixed values for the electric potential and ion concentrations were set on the upper and lower faces of the cubic box, and the protein channel axis was oriented perpendicular to those faces. As an initial "guess," the electric potential and ion concentration were taken to vary linearly along the channel axis between the fixed values at the faces. At the lateral walls of the box, a linear variation of the electric potential between the upper and lower faces and zero-concentration boundary conditions were imposed. It was verified that computed channel properties (e.g., I-V curves) were insensitive to lateral box face boundary conditions when the lateral box dimension significantly exceeded the width of the permeation channel.

The profiles were then updated according to the self-consistent SOR procedure outlined above. The concentrations were allowed to change in the entire solvent region (channel and reservoirs). Therefore, the condition of electrolyte charge neutrality was imposed only at the faces of the box. For the ion diffusion coefficients, a linear interpolation between bath and internal channel values was adopted (details are provided below).

Runs were performed on a DEC 21164a-clone and a set of IBM RS6000 workstations. Converged results for a point on the current-voltage curve took several hours, with the exact time depending on system details and initial conditions.

	PNP LATTICE MODEL OF THE GRAMICIDIN A CHANNEL

TOP ABSTRACT INTRODUCTION PNP THEORY COMPUTATIONAL IMPLEMENTATION PNP LATTICE MODEL OF... INFLUENCE OF MEMBRANE CHARGES... COMPARISONS WITH EXPERIMENTAL... CONCLUSIONS APPENDIX REFERENCES

Gramicidin A (GA) is a 15-amino acid peptide that dimerizes to form a monovalent cation-selective channel in the bacterial cell wall or in artificial lipid membranes (Venkatchalam and Urri, 1983; Wallace, 1990; Andersen and Koeppe, 1992; Busath, 1993). The secondary structure of the GA channel is a -helix head-to-head dimer, comprising two identical subunits, A and B (cf. Fig. 1 a), which forms a narrow water-filled pore. The 3D structure of the dimer is known from 2D NMR and NOE spectroscopy studies to a resolution of 0.86 Å (Arsen'ev et al., 1986). The GA chain consists of amino acids with alternating L and D stereochemistry, which permits nonpolar side groups to extend into the membrane while the pore is lined by the polar backbone peptide groups (see Fig. 1 a).

View larger version (27K): [in this window] [in a new window]   FIGURE 1   (a) Molecular representation of gramicidin A dimer, with the two peptides ("monomers") in dark gray (A) and light gray (B), respectively. (b) A 2D cut of the lattice representation of the gramicidin A channel in PC or PS lipid. The protein and membrane regions are in light gray and the solvent is shown in white. The entrance and exit regions are represented by a funnel-like pore with an external radius of 15 Å. (c) The same 2D cut for the lattice representation of GA channel in GMO lipid. The channel pore is taken to be cylindrical, with a radius of 2.5 Å.

Single-channel cation current for gramicidin A has been studied under a variety of conditions (Aidley and Stanfield, 1996; Hille, 1992; Busath, 1993). The abundance of experimental data available for the GA channel together with its small size have made it the preferred choice for testing theories of channel conductance (Barcilon et al., 1992; Roux and Karplus, 1993; Andersen and Feldberg, 1996). In previous work (Kurnikova et al., 1999), we used the GA channel to test the reliability and performance of our 3D PNP algorithm. Reasonable agreement between theoretical and experimental results for current-voltage curves was obtained. In the present work, we have refined our PNP calculation by including relaxation of ion concentrations in the bathing solutions and allowing for different diffusion constants in the channel interior and exterior. This enables us to study the influence on GA conductance of charged groups and interfacial dipoles that reside on the membrane surface and to compare our calculations with available experimental results (Rostovtseva et al., 1998; Busath et al., 1998).

Fig. 1 b shows how the gramicidin A channel, surrounding membrane (PC/PS bilayer model), and solvent baths are represented on the grid. This setup is somewhat different from the one used previously (Kurnikova et al., 1999). In the present model, the protein molecule is completely embedded in the membrane region, and the upper and lower faces of the computation box (abutting the two entrances to the pore region) are farther away from the channel mouths to allow ion density relaxation in the bath regions. For this 47-Å-wide lipid bilayer we assume that the membrane compresses because of hydrophobic matching between the lipid and the protein, as has been determined for another phospholipid bilayer (Harroun et al., 1999a,b). The exact details of PC/PS shrinkage are unknownwe used a funnel-like geometry for computational simplicity. To model the thinner (33 Å) GMO bilayer, which is characterized by a nondipolar surface, we used a cylindrically shaped membrane pore at the openings of the channel with a radius of 2.5 Å, as depicted in Fig. 1 c.

The coordinates of the heavy atoms of the protein (Arsen'ev et al., 1986) were taken from the Protein Data Bank (Bernstein et al., 1977). For the partial charges on GA atoms we used values from the AMBER86 united-atom force field (Pearlman et al., 1991), while values for the atomic radii were taken from DelPhi (Nicholls et al., 1990); the radii of the polar hydrogens were set at 1.0 Å. In the present calculations, the membrane and protein regions (gray area in Fig. 1, b and c) are described by the low dielectric constant _m = 2. The high dielectric constant _a = 80 characterizes the aqueous (channel and bath) region (white region of Fig. 1, b and c). For a more detailed discussion concerning the choice of these values see Kurnikova et al. (1999).

To obtain the results presented in Comparisons with Experimental Observations, where we compare our calculations with experimental results, we used as bulk diffusion constants for K⁺, Cs⁺, Cl, and H⁺ the known experimental values of 1.96 × 10⁵, 2.06 × 10⁵, 2.03 × 10⁵, and 9.31 × 10⁵ cm²/s, respectively (Hille, 1992). In the next section, the main purpose of which is to establish qualitative effects of surface charge and/or surface dipoles on channel properties (electrical potential and mobile ion distributions, I-V curves, etc.), we use a bulk diffusion constant of 1 × 10⁵ cm²/s for both the positive and negative ions. In general, values of diffusion coefficients inside the channel are difficult to obtain experimentally. A molecular dynamics simulation of ion transport inside a cylindrical channel yielded diffusion coefficients two or three times smaller than the bulk values (Lynden-Bell and Rasaiah, 1996). Another recent molecular dynamics study of permeation through several different ion channels found an ~10-fold reduction of the diffusion coefficient inside narrow channels (Smith and Sansom, 1998). To obtain agreement with experimental GA results (for the case of an uncharged membrane), channel diffusion coefficients 11 and 17 times smaller than the bulk values were required for Cs⁺ and K⁺ ions, respectively. This reduction is consistent with the decrease estimated by Smith and Sansom for Na⁺ ions (which are somewhat smaller than Cs⁺ or K⁺) permeating the narrow poly-Ala -helix bundle model (cf. Figure 3B of their paper, where reductions of ~15 times are indicated). Roux and Karplus found even larger reductions of the internal channel diffusion coefficient in their molecular dynamics simulation of permeation through gramicidin (Roux and Karplus, 1991), thus reinforcing the basic premise behind our modeling of the spatial diffusion constant profile in the present work.

	INFLUENCE OF MEMBRANE CHARGES AND INTERFACIAL DIPOLES ON GA CONDUCTANCE

TOP ABSTRACT INTRODUCTION PNP THEORY COMPUTATIONAL IMPLEMENTATION PNP LATTICE MODEL OF... INFLUENCE OF MEMBRANE CHARGES... COMPARISONS WITH EXPERIMENTAL... CONCLUSIONS APPENDIX REFERENCES

The ion transport characteristics of protein channels can be modified by differences in electrical potential at the bulk water-membrane surface induced by dipoles lining the membrane surface. These dipoles are created by a nonrandom orientation of lipid headgroups, fatty acid carbonyl groups, and water (Gawrisch et al., 1992). Another factor that can strongly affect the channel conductance is the presence of charged polar groups in some lipid bilayer surfaces forming ion-channel systems (Green and Andersen, 1991). Certain lipids like phosphatidylserine contain both negatively charged carboxylic oxygens and interfacial dipoles on their surface. The ion current of GA embedded in this membrane (PS) is larger than the corresponding current through the uncharged phosphatidylcholine (PC) (Rostovtseva et al., 1998).

Here we use 3D PNP theory to study ion conduction through the GA channel under several membrane electrostatics conditions. Specifically, we compare the conductance properties of a charged membrane, a dipolar membrane, a charged/dipolar membrane, and a neutral membrane.

We modeled the charged membrane by including negatively charged "dummy" atoms on the surfaces of the bilayer (Fig. 2 a, black spheres). The positions of these atoms were attached to the coordinate file for the GA atoms, and their number was chosen to correspond to a surface charge density of 0.021e/Ų (a value estimated for the PS membrane by Rostovtseva et al. (1998)). The specific geometrical arrangement chosen, namely three concentric squares (cf. Fig. 2 b, black spheres), is somewhat arbitrary because of the lack of knowledge about the structure of lipids in bilayer membranes.

View larger version (49K): [in this window] [in a new window]   FIGURE 2   Molecular representation of gramicidin A dimer with negative charges (black) and dipoles (light gray) embedded in the membrane. (a) Lateral view. (b) Top view. The negative charges (84 in total) and the positive charges of the dipoles (360 in total) are placed inside the membrane (+ sign in a). The negative charges of the dipoles are placed on the aqueous side of the membrane-liquid interface ( sign in a).

To model the interfacial dipoles at the membrane surface we added dipoles at each side of the bilayer membrane (Fig. 2 a, light gray spheres). It has been determined for a wide range of membrane structures that the potential is positive inside the membrane with respect to the aqueous phase (Haydon and Myers, 1973; Jordan, 1984). Accordingly, we placed the positive ends of the dipoles inside the membrane and separated this layer of charge by 5 Å from a negatively charged layer, comprising the negative ends of the dipoles, which abut the bath. The absolute values of the membrane potential can be estimated from a simple capacitor model. For such a dipolar interface, the dipole potential V_d can be approximated as (Flewelling and Hubbell, 1986)

(7)

where _µ is the density of dipoles at the surface (in Ų), µ is the individual dipole moments (in Debyes), and _eff is an effective dielectric constant at the interface, yielding V_d in volts. We used a charge of ±0.031e on the ± components of each dipole and a dipolar density of 1 dipole per 9.4 Å² (cf. Fig. 2 b, light gray spheres). Assuming _eff = 20 (a value intermediate between the bulk value for water and a rigid structure; Flewelling and Hubbell, 1986), Eq. 7 gives a dipole potential of ~160 mV for the parameters given above. This value is similar to the experimental value of (+)120-200 mV for phosphatidylcholine (PC) lipid (Jordan, 1984; Flewelling and Hubbell, 1986; Busath et al., 1998) (compared with glycerylmonoolein (GMO) lipid) (Haydon and Myers, 1973; Busath et al., 1998).

We performed the computations of this section using the membrane model of PC shown in Fig. 1 b. Calculations were made for a membrane without charges or interfacial dipoles (neutral membrane), with dipoles (no charges), with charges (no dipoles), and with both dipoles and charges (cf. Fig. 2). As noted at the end of the previous section, we used a bulk diffusion coefficient of D₊ = D = 10⁵ cm²/s. To account for the expected reduction of permeation mobility inside the pore a channel diffusion coefficient of D₊ = D = 10⁶ cm²/s was employed. A linear variation of the diffusion constant between these two values in the funnel-shaped openings connecting the bulk solution to the "inner" (protein) channel was imposed. We checked that this choice of interpolation scheme for connecting internal and external flow region, which is necessarily somewhat arbitrary, does not qualitatively alter the computed channel properties.

In Fig. 3, current-voltage curves for GA embedded in the neutral, dipolar, charged, and dipolar/charged membrane are shown for a symmetrical 0.2 M bathing solution. The presence of negative charges on the membrane surface increases the channel ion conduction compared with the neutral membrane, in agreement with experimental observations (Apell et al., 1979; Rostovtseva et al., 1998). The presence of a dipolar interface decreases the ion conductance relative to that of the corresponding neutral membrane, also in agreement with experiment (Busath et al., 1998), but this effect is moderate compared with the increase in current due to the presence of charges in the membrane.

View larger version (16K): [in this window] [in a new window]   FIGURE 3   Current-voltage curves for GA embedded in a neutral membrane (solid line), dipolar membrane (dashed line), charged membrane (long dashed line), and charged and dipolar membrane (dot-dashed line). The PNP equations were solved for a bulk salt concentration of 0.2 M on both sides of the channel and diffusion coefficients D+ = D = 105 cm2/s (bulk) and D+ = D = 106 cm2/s (channel).

The electrostatic potential and positive and negative ion concentrations along the channel axis are shown in Figs. 4 and 5, respectively, for the different surface dipole/charge configurations under consideration at an applied voltage of 100 mV. The presence of negative charges on the membrane surface deepens the potential at the entrance and exit from the channel by ~29 mV; it also lowers the potential well inside the channel by ~18 mV (Fig. 4). Consequently, more positive ions are attracted into the pore (Fig. 5), thus increasing the ion current. The presence of dipoles at the membrane-liquid interface produces a reduction in the depth of the potential well (by ~5 mV) with a concomitant decrease of positive ion density inside the pore (Fig. 5). This reduction of the electric potential has previously been pointed out by Jordan (1984) to explain the effect of membrane dipoles on channel conductance. The dipoles also destabilize the entrance of positive charges into the pore (manifested as a small bump in the electrostatic potential curve between the membrane surface and the actual protein pore entrance). The positive side of the interfacial dipoles located in the membrane interior leads to this increase in the electric potential, which in turn reduces the ion current through the channel. In all cases, the anion distribution density inside the channel is negligible compared with the cation density, in agreement with the observed cationic selectivity for GA channels. Thus the issue of shielding of cations inside the channel by negative counterions (Corry et al., 1999) does not arise here.

View larger version (21K): [in this window] [in a new window]   FIGURE 4   Electrostatic potential along the channel axis. The symbols are the same as in Fig. 3. The vertical lines indicate the position of the membrane surfaces (external lines) and the protein mouths (internal lines). The applied external potential drop is 100 mV, and the bulk concentration is 0.2 M.

View larger version (26K): [in this window] [in a new window]   FIGURE 5   Positive and negative ion concentrations along the channel axis (the anion concentrations are essentially zero inside the channel for all four cases). The symbols are the same as in Figs. 3 and 4. The external potential drop is 100 mV, and the bulk concentration is 0.2 M.

In the following section we test the validity of these membrane electrostatic models and of 3D PNP theory by comparison with recent experimental studies on GA conductance in different membrane environments.

	COMPARISONS WITH EXPERIMENTAL OBSERVATIONS

TOP ABSTRACT INTRODUCTION PNP THEORY COMPUTATIONAL IMPLEMENTATION PNP LATTICE MODEL OF... INFLUENCE OF MEMBRANE CHARGES... COMPARISONS WITH EXPERIMENTAL... CONCLUSIONS APPENDIX REFERENCES

Titration of membrane charge in PS

It is well known that charged polar groups on the surface of lipids can interact with ions in solution and thus modify the physical properties of the bilayer system (Eisenberg et al., 1979). Recently, Rostovtseva and co-workers (Rostovtseva et al., 1998) used this fact to manipulate the surface charge density of the membrane. In one set of experiments, they titrated the lipid charge by changing the pH of the bulk solution. In another set, they modified the membrane charge density by diluting the charged lipid PS with uncharged PC. They then utilized a model based on Gouy-Chapman theory to study the surface titration process. Taking into account chemical binding of counterions to the negatively charged lipid groups and the corresponding decrease in effective membrane charge, they estimated a value for the intrinsic dissociation constant of the PS polar headgroups (see below). They used three different models to determine the cation concentration at the channel mouths: an analytical solution for the planar Poisson-Boltzmann equation, a linearized Poisson-Boltzmann equation solution, and construction of a Gibbs dividing surface between the solvent bath and surface membrane. The conductance versus pH predictions of this model were ultimately based on the measured conductance of GA in neutral PC lipid (cf. equation 10 of Rostovtseva et al. (1998)).

Here we use a 3D PNP algorithm to study the effect of titration on the ion conductance of the GA channel. Any ambiguity in the determination of the electric potential at the surface of the charged membrane and at the center of the channel mouths is eliminated when the PNP equations are self-consistently solved in the 3D system. Another advantage of the 3D PNP approach is that its predictions of the pH dependence of ion current through a GA channel in the PS membrane are independent of the ion conductance (measured or calculated) in neutral PC lipid.

We modeled the uncharged PC and charged PS lipids as indicated in Fig. 1 b. We assumed that the two membrane bilayers have the same width (47 Å) and the same dipolar interface (Rostovtseva and co-workers made this assumption implicitly because they did not include effects of the dipolar interface in their conductance study of PC and PS). These two assumptions can be justified by the similarities between the zwitterionic headgroup of phosphatidylcholine and phosphatidylserine and the approximately equal sizes of the hydrocarbon chains in the two phospholipids. Dipoles were included on both sides of the membrane bilayer model for PC and PS. To model the charged PS we embedded surface charges in the membrane (Fig. 2). The same values of surface dipolar and charge density utilized in the previous section were adopted here. We did not include interfacial dipoles or surface charges along the funnel-like entrance/exit to the channel because the structural modifications of the lipid molecules due to the membrane shrinkage at the pore are not known.

We determined the current-voltage relation for GA in both charged (PS) and uncharged (PC) membranes. Results based on a channel diffusion constant of 1.79 × 10⁶ cm²/s for both cations and anions are shown in Fig. 6. This diffusion constant value was chosen to fit the experimental data for GA/neutral membrane conductance in a 1 M CsCl solution (Rostovtseva et al., 1998) (inset of Fig. 6, filled triangles). The increase in current for GA embedded in the charged membrane compared to that obtained for the same ion channel embedded in an uncharged membrane under otherwise identical conditions is in good agreement with the experimental data. The negative charges in the lipid bilayer surface decrease the electrostatic potential along the channel axis (cf. Fig. 4), which attracts cation density to the pore (Fig. 5), thus increasing the ion current compared with the case of an uncharged (but still dipolar) membrane.

View larger version (24K): [in this window] [in a new window]   FIGURE 6   Calculated current-voltage relationship for GA embedded in an uncharged PC membrane (solid symbols) and a charged PS membrane (open symbols) at neutral pH. The electrolyte concentrations are 0.1 M (diamonds) and 1.0 M (triangles). The inset shows experimental results (Rostovtseva et al., 1998), using the same symbol convention as in the main panel (it also shows experimental results at pH 1 (solid circles and squares)).

To study the titration of the negative carboxyl groups of PS lipid with bulk solution protons we follow the analysis of Rostovtseva et al. relating the membrane charge density to the local counterion concentrations and dissociation constants. Cs⁺ and H⁺ can bind to the carboxyl groups of the bilayer PS lipids:

with dissociation constants

(8)

Parentheses here denote two-dimensional surface concentrations (with subscript s), and 3D concentrations of mobile ions at the surface interface are labeled with subscript int. The lipid surface charge density  is

(9)

and the maximum charge density ^max is proportional to the concentration of carboxyl groups in PS:

(10)

Combining Eqs. 8-10, the following expression for the charge density at the lipid surface is obtained (Rostovtseva et al., 1998):

(11)

where [...]_b denotes bulk concentration and it has been assumed that [H⁺]_int/[Cs⁺]_int = [H⁺]_b/[Cs⁺]_b. Equation 11 relates the charge density of the membrane to the bulk concentrations of the metal ion and proton, the local concentration of Cs⁺ at the interface of the surface membrane and in solution, and the dissociation constants K_a and K_Cs. The bulk concentrations are known experimentally, but the electrolyte concentrations near the surface and dissociation constants are not. Rostovtseva et al. used a simple 1D solution of the PB equation for the local concentration, using the Gouy-Chapman expression to estimate the electrical potential at the surface/solution interface. Here we solve the 3D PNP equations self-consistently with Eq. 11 for the surface charge density (which depends implicitly on the local concentration of mobile positive ions). We used a value of 6.0 × 10⁵ cm²/s for the proton diffusion coefficient in the GA channel. This value produced agreement with the channel conductance measured in uncharged PC bilayers at pH 1.0 in the absence of Cs⁺ (table 2 of Rostovtseva et al.). For K_Cs we used the value 20 M (Eisenberg et al., 1979), and we varied the value of K_a to fit the experimental data as discussed below.

In Fig. 7 a, the pH dependence of the GA channel conductance in the charged membrane calculated with the 3D PNP algorithm is compared with the experimental result. (We included the proton current in addition to the metal ion current for our calculation of the pH dependence of channel conductance; for the rest of the calculations presented in this section, which were carried out in pH neutral solutions, only metal ion current was included.) Our results reproduce the experimental curve rather well in the pH range between 1.5 and 6, with more deviation at lower pH. There is a slight reduction of the conductance as the pH is reduced from 8.0 to ~2.75, which can be traced to a decrease in the effective membrane charge density caused by protons binding to the carboxyl groups (Fig. 7 b). Below pH 2.75 the conductance increases sharply, even though the net surface charge density continues to decrease (Fig. 7 b), because of the rise in proton concentration and hence in the proton current (the proton has a larger diffusion constant than the Cs⁺ ion). We obtained best fits to the experimental data when K_a = 1.8, which is somewhat smaller than the value obtained by Rostovtseva et al. (K_a = 2.5).

View larger version (17K): [in this window] [in a new window]   FIGURE 7   (a) Conductance and (b) effective negative surface charge density, calculated using 3D PNP theory versus pH for a gramicidin A channel embedded in a PS membrane. The bulk electrolyte concentration is 0.1 M, and the conductance values were calculated for a potential difference of 100 mV. Shown in the inset of a (open and filled circles) are the experimental values (Rostovtseva et al., 1998).

We also studied the effect of titrating the charged PS lipid by mixing in neutral PC lipid. The conductance of GA varies because of the change in membrane charge density induced by "diluting" the membrane with uncharged lipid. Surface charge densities at several PS-PC compositions were estimated by Rostovtseva et al. with the nonactin method (McLaughlin et al., 1970). As in the case of pure PS, we modeled these mixtures by embedding charged particles on both surfaces of the (already dipolar) bilayer membrane to get the appropriate charge densities at different PS-PC lipid compositions and then calculated the conductance. In Fig. 8, we show the results obtained with our 3D PNP algorithm. The agreement with the experimental curve (shown in the inset of the figure) is reasonable. The calculated conductance is in general smaller than the corresponding experimental value, especially for the uncharged, PC-rich mixtures. Nevertheless, the qualitative similarities between the results of our model and experiment suggest that the variations in ion current are induced mainly by the changes in membrane charge and less so by any additional lipid-dependent structural factors.

View larger version (21K): [in this window] [in a new window]   FIGURE 8   Ion conductance for gramicidin A in mixed bilayers as a function of PS/(PC + PS) ratio in a pH neutral solution with 0.1 M CsCl. These conductances were calculated for an applied external potential of 100 mV. The corresponding experimental results (Rostovtseva et al., 1998) are shown in the inset.

I-V curves for GA in neutral and dipolar membranes

Busath et al. (1998) measured K⁺ and Na⁺ conductance through a GA channel embedded in planar diphytanoylphosphatidylcholine (DPhPC) and GMO bilayers. The interfacial dipole potential is significantly larger for DPhPC than for GMO lipid (cf. previous section). We performed 3D PNP calculations, using our dipolar membrane model of DPhPC (Figs. 1 b and 2), and for a model of GMO, which comprises a neutral membrane model 33 Å wide containing a narrow cylindrical pore of radius 2.5 Å (cf. Fig. 1 c). (The membrane shrinkage is probably smaller for this bilayer, which has a width only slightly greater than the length of GA (de Planque et al., 1998).) As in the case of the DPhPC membrane, a linear interpolation between bulk and inner (GA) channel diffusion constants was imposed, although for GMO the spatial extent of the interpolation region is much smaller (cf. Fig. 1, b and c). The calculated current-voltage curves for different salt concentrations are plotted in Figs. 9 and 10 for the dipolar (DPhPC) and neutral (GMO) membrane, respectively. A channel diffusion coefficient of 1.12 × 10⁶ cm²/s was used for K⁺. (This value fits the calculated I-V curve with the experimental value in neutral GMO lipid for a cation concentration of 1 M.) The corresponding experimental I-V curves (Busath et al., 1998) are shown as insets in the figures. The reduction in channel current observed with the dipolar DPhPC lipid relative to neutral GMO lipid is well reproduced by our model. The interfacial dipoles increase the electrostatic potential along the channel and at the pore openings (cf. Fig. 4), which results in decreased flux of positive ions through the channel. The difference in membrane width also affects the I-V relationship: compare the rather modest effect of the dipoles in Fig. 3 with the more significant dipole-induced changes observed in Figs. 9 and 10. These latter can be traced to a >10% increase in current through the narrow GMO pore relative to the current obtained for the "hypothetical" nondipolar DPhPC considered in Fig. 3. Our results also show a very slight superlinearity and sublinearity of the I-V curves at high and low concentrations, respectively, compared with the experimental results. Saturation and superlinear behavior of the experimental I-V curve in DPhPC (and superlinearity in GMO) at concentrations 2 M and higher, probably caused by multiion interference, is not reproduced by primitive PNP theory, which assumes a continuum description of the (infinitesimal) ions in the system.

View larger version (24K): [in this window] [in a new window]   FIGURE 9   Current-voltage relationship in DPhPC (dipolar) membrane. The electrolyte concentrations are 0.1 (open triangle), 0.2 (closed circle), 0.5 (open square), 1.0 (open circle), and 2.0 M (closed square). The inset shows the experimental results (Busath et al., 1998), with the same symbol convention, except that these authors used dot-filled squares for 0.1 M.

View larger version (27K): [in this window] [in a new window]   FIGURE 10   Current-voltage relationship in GMO (nondipolar) membrane. The symbol legend is the same as in the previous figure. The experimental results (Busath et al., 1998) are shown in the inset.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION PNP THEORY COMPUTATIONAL IMPLEMENTATION PNP LATTICE MODEL OF... INFLUENCE OF MEMBRANE CHARGES... COMPARISONS WITH EXPERIMENTAL... CONCLUSIONS APPENDIX REFERENCES

In this paper we have utilized a 3D PNP algorithm to study the effect of membrane surface charges and dipoles on the conductance properties of the gramicidin A channel. Good agreement with a range of experimental results was obtained. Three-dimensional PNP calculations can provide valuable insights concerning the influence of membrane electrostatics on ion flux through the bulk solution and channel.

We have found, for example, that the observed current decreases when dipoles are included on the membrane surfaces, but these dipoles have only a small effect on ion permeation through the GA channel (all other factors being equal). By switching the dipoles off in the DPhPC membrane (a useful exercise, which is more easily done in a numerical simulation than a laboratory experiment!), we find that the observed difference in the currents through GA/GMO and GA/DPhPC systems is in large part due to the difference in the thickness of the GMO and DPhPC membranes, indicating that it is important to take into account the specific details of the membrane to understand ion flow through these systems.

Charges on the membrane surface have a larger effect on the ionic conductance through the GA channel. The observed increase in current when there are negative charges on the membrane surface was traced to a drop in the electrostatic potential inside the channel and the channel mouths. The rather good agreement between experimental and PNP results in this case suggests that the difference in conductance between GA/PC and GA/PS (modeled with the same membrane thickness) is due mainly to the presence of charges on the surface of the latter.

Nevertheless, limitations of PNP theory should be kept in mind. These include a continuum description of the solvent and permeant ions, treating the latter as infinitesimal point particles, and an intrinsic reliance on a mean field approximation that disregards ion-ion correlations (Corry et al., 1999). These limitations render models based on primitive PNP theory incapable of predicting properties that depend critically on the correlated motion of finite-size particles, such as multiion interference effects observed at high salt concentrations (see, for example, Fig. 9). Furthermore, primitive PNP does not account for the energetic barrier encountered by a mobile ion when it moves through a narrow channel, which is a consequence of the finite radius of the ion. When there is a large dielectric discontinuity between the aqueous pore and the surrounding protein channel/membrane structure, this barrier can be large (Dieckmann et al., 1998, Chung et al., 1998, 1999), particularly for narrow channels. The success of 3D PNP modeling in accounting for a range of experimentally observed permeation data on the GA channel thus suggests that long-range structural electrostatic forces (induced by the shape of the embedding membrane and the distribution of charges attached to it) significantly influence the dynamics of the permeant ions, and that such considerations, rather than delicate atomistic-level details of the migration of particles through the narrow channel, dominate the behavior of observable I-V curves in GA systems. Another general source of difficulties with PNP theory (and others, too; Chung et al., 1998, 1999; Corry et al., 1999) is the lack of experimental knowledge concerning the values of some parameters, e.g., the dielectric constants of water in the channel and the diffusion coefficients for ions inside the channel, as well as the lack of high-resolution information concerning the lipid bilayer structures. Despite these limitations, our results suggest that solutions of the 3D PNP equations, combined with intuition (based on experimental evidence) concerning structural details of the protein channel/membrane system, may be useful for understanding a range of biological ion transport processes.