INTRODUCTION

TOP ABSTRACT INTRODUCTION SIMULATION ALGORITHM COMPUTATIONAL ILLUSTRATIONS CONCLUDING DISCUSSION APPENDIX A APPENDIX B REFERENCES

In recent years, molecular dynamics (MD) simulations with explicit solvent and membranes have provided an unprecedented amount of information about a number of ion channels such as gramicidin A (Roux and Karplus, 1994; Woolf and Roux, 1994, 1996, 1997; Chiu et al., 1999), LS3 (Zhong et al., 1998a,c), alamethicin (Tieleman et al., 1999), the transmembrane domain of the influenza A virus M2 coat protein (Sansom et al., 1997; Zhong et al., 1998b), OmpF Escherichia coli porin (Suenaga et al., 1997, 1998; Tieleman and Berendsen, 1998), and KcsA (Guidoni et al., 1999; Shrivastava and Sansom, 1999; Allen et al., 1999; Bernèche and Roux, 2000). Such computer simulations arguably represent the most detailed method to study complex biomolecular systems. Nevertheless, despite the progress in computational methodologies, investigations of ion channels are still confronted with difficult fundamental problems (Roux, 1998). One of the main problem is due to the time scale of the permeation process: the translocation of a single ion across a channel takes on the order of a microsecond (Hille, 1992), which is extremely long compared to the typical length of calculated MD trajectories. To simulate the permeation process it is necessary to rely on approaches that are simpler and computationally less expensive than fully detailed MD.

Brownian dynamics (BD) provides an attractive computational approach for simulating the permeation process over long time scales without having to treat all the solvent molecules explicitly (Cooper et al., 1985). The approach, which consists of integrating stochastic equation of motions with some effective potential function (Ermak, 1975), is physically reasonable because the movement of ions in aqueous solution is essentially chaotic and diffusive (see "Elementary properties of ions in solution" in Hille (1992) and references therein). There is vast experience with the application of BD simulations to macromolecular polyions (Simon and Zimm, 1969), ionic solutions (Ermak, 1975; Wood and Friedman, 1987; Jardat et al., 1999), enzyme-substrate encounters (Wade et al., 1994; Madura et al., 1995), and ion channels (Läuger and Apell, 1982; Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989; Bek and Jakobsson, 1994; Li et al., 1998; Chung et al., 1998, 1999; Schirmer and Phale, 1999). Nonetheless, important difficulties with BD simulations of ion channels are still unresolved.

A first difficulty concerns the implementation of the boundary conditions corresponding to the local ion concentration in the bulk region at the channel entrance. The channel is an open system with a fluctuating number of ions and the simulations require a modelization of the entrance and exit of ions at the channel mouth. The "driving force" responsible for a steady state of ion flow with non-equilibrium boundary conditions of concentrations cannot be represented directly in terms of a mechanical force, because such a force arises effectively from the statistical difference in the number of particles per unit volume. A net ion flux is established because it is more likely for ions to enter the channel from the side with the higher concentration. To simulate concentration effects, Jakobsson and co-workers designed the "entrance tube" algorithm by means of which one ion is inserted randomly into a tube of variable length at the entrance of the channel (Cooper et al., 1985; Jakobsson and Chiu, 1987). The entrance tube algorithm was used to simulate one-ion (Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989) and multi-ion pores (Bek and Jakobsson, 1994). In those simulations, the ion movements were restricted to one dimension (1D) along the pore axis. A related theory, based on a random walk in 1D, was developed by McGill and Schumaker (1996). BD simulations of ion permeation were extended to three dimension (3D) by Chung and co-workers (Li et al., 1998; Chung et al., 1998, 1999) for channel models of simplified shapes. In those simulations, the total number of ions was fixed. To constrain the ion concentration in the bulk region under conditions of steady flow, ions are removed randomly from the bulk region on one side of the membrane and taken to the bulk region on the other side. Recently, Schirmer and Phale (1999) simulated the trajectories of individual Na⁺ and Cl ions going through atomic models of Escherichia coli porins using the program UHBD (University of Houston Brownian Dynamics) developed by McCammon and co-workers (Madura et al., 1995); one ion at a time was considered in those simulations and finite-concentration effects were not modeled explicitly.

A second difficulty with BD of ion channels concerns the implementation of the transmembrane potential. At the microscopic level, the transmembrane potential arises from a very small imbalance of net charges of the mobile ions in the bulk solutions on each side of the membrane (interfacial polarization). A modified Poisson-Boltzmann (PB) equation based on statistical thermodynamics was developed to calculate the transmembrane potential (Roux, 1997, 1999). Detailed numerical calculations using the modified PB equation showed that the transmembrane potential is roughly linear in the case of a narrow pore such as the gramicidin channel (Roux, 1999). This calculation provides a validation of the constant field approximation that had been used in several BD simulations (Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989; Bek and Jakobsson, 1994; McGill and Schumaker, 1996). The linear field is, however, probably inaccurate in the case of wide aqueous pores and it is then necessary to rely on the numerical solution of the modified PB equation. The transmembrane potential and its relation to continuum electrostatic theory is often misunderstood. For example, in the BD simulations of Chung and co-workers the transmembrane potential is calculated as the solution to Poisson's equation with boundary values (Li et al., 1998). However, such treatment is inconsistent with the microscopic origin of the transmembrane potential as an interfacial polarization phenomenon, because it does not take into account the influence of mobile counterions in the solution.

Our goal is to formulate and develop a non-equilibrium simulation algorithm enabling one to simulate ion permeation through ion channels. A valid non-equilibrium simulation algorithm must result in a satisfactory description of the equilibrium state of the ion channel in the absence of net fluxes. However, none of the previously proposed BD algorithms for simulating ion channels (Jakobsson and Chiu, 1987; Chiu and Jakobsson, 1989; Bek and Jakobsson, 1994; McGill and Schumaker, 1996; Li et al., 1998; Chung et al., 1998, 1999; Schirmer and Phale, 1999) satisfies this requirement. The recently formulated statistical mechanical equilibrium theory of ion channels (Roux, 1999) provides a sound basis for establishing a correct algorithm for simulating non-equilibrium ion transport across membrane channels. In particular, it was shown that the number of ions in the channel follows the statistics of an effective Grand Canonical Ensemble. To simulate ion channels with appropriate boundary conditions allowing fluctuations in the number of ions, we combine the Brownian dynamics (BD) simulation method with the Grand Canonical Monte Carlo (GCMC) algorithm (Norman and Filinov, 1969; Allen and Tildesley, 1989). We call the combined method the GCMC/BD algorithm. The present GCMC/BD algorithm is closely related to the dual-volume-control molecular dynamics method (DCV/MD) (Heffelfinger and Ford, 1998; Thompson et al., 1998; Pohl and Heffelfinger, 1999; Thompson and Heffelfinger, 1999), which has been used to study diffusion problems in chemical physics.

In the next section, the algorithm is developed and is illustrated with some simple examples and with a simulation of a detailed model of the OmpF porin of E. coli in 1 M KCl. The paper is then concluded with a brief discussion.

	SIMULATION ALGORITHM

TOP ABSTRACT INTRODUCTION SIMULATION ALGORITHM COMPUTATIONAL ILLUSTRATIONS CONCLUDING DISCUSSION APPENDIX A APPENDIX B REFERENCES

Ion dynamics

We consider a system of volume V corresponding to the neighborhood of an ion channel embedded in a lipid membrane and in contact with surrounding aqueous salt solutions. In the bulk solution, the density of ions of type  ( = 1, 2, ...) is , and the excess chemical potential is . We are seeking a simple algorithm to simulate the dynamical movements of the ions in the channel system. Because the system is in contact with a bulk solution, the number of ions in the volume V can vary with time.

Let us assume for the sake of simplicity that the number of ions in the system is fixed. The Cartesian coordinates of the n ions of type  in the volume V are represented by R  (r⁽¹⁾, r⁽²⁾, ... , r^(m)). (For example, the position of the second ion of type 1 is given by r₁⁽²⁾.) Assuming that the ions are initially in the configuration {R₁(t), R₂(t), ...}, we wish to predict the position of the ions a short time later. If all the degrees of freedom were explicitly included in the theory, the deterministic trajectory of the system would simply obey Newton's law of classical mechanics. However, although all the ions in the system are present explicitly in order to account accurately for ion-ion interactions, it is desirable to incorporate the influence of the solvent, channel, and membrane implicitly to allow long time-scale simulations. If those degrees of freedom relax rapidly compared to the movements of the ions, a physically reasonable description of the dynamics is provided by the Langevin equation (Chandrasekar, 1943; McQuarrie, 1976)

(1)

where m is the mass, is the friction coefficient, F⁽ⁱ⁾ is the effective microscopic net force, and f⁽ⁱ⁾(t) is a random Gaussian force obeying the fluctuation-dissipation theorem f⁽ⁱ⁾(t) · f⁽ⁱ⁾(0) = 6k_BT(t). When the friction is large and the motions are overdamped, the inertial term (mass times acceleration) may be neglected. This approximation yields the BD equation (Chandrasekar, 1943)

(2)

where D = k_BT/ is the diffusion constant of ions of type , and ⁽ⁱ⁾(t) is a random Gaussian noise with ⁽ⁱ⁾(t) · ⁽ⁱ⁾(0) = 6D(t). Equation 2 can be integrated numerically with a finite time step using the algorithm proposed by Ermak (1975) to generate a BD trajectory.

The effective microscopic force acting on the ith ion, F⁽ⁱ⁾, is directly related to the gradient of the multi-ion potential of mean force (PMF)

(3)

The PMF is a statistical mechanical concept that plays an important role in the description of complex molecular systems at equilibrium (Kirkwood, 1935), and at non-equilibrium (Chandler, 1978). In the context of ion channels, the multi-ion PMF, (n₁, n₂, ... ; R₁, R₂, ...), corresponds to the reversible thermodynamic work function (free energy) to insert the (n₁, n₂, ...) ions at their position (R₁, R₂, ...) into the system while averaging over the configuration of the solvent molecules, channel, and membrane (Roux, 1999). It should be stressed that the mean force F⁽ⁱ⁾ is expressed explicitly in terms of all ions present (i.e., there is no averaging over the position of the ions in the channel system). However, the ions located outside the volume V are "averaged out" and their influence is taken into account implicitly by the multi-ion PMF. This is a direct consequence of the Heaviside function H_p(X_r) appearing in Eq. 26 of Roux (1999), which prevents the ions other than (R₁, R₂, ...) from entering into the pore region.

BD trajectories generated with Eq. 1 or Eq. 2 with a fixed number of ions cannot describe the permeation process in a satisfactory manner. In reality, the total number of ions must be able to fluctuate under the influence of specific non-equilibrium boundary conditions because ions can enter and leave the system.

Grand canonical Monte Carlo

It has been shown previously that under equilibrium conditions the variable number of ions in the volume V follows the effective Grand Canonical Ensemble (GCE) (Roux, 1999)

(4)

An important prerequisite in constructing a valid non-equilibrium simulation algorithm is that the system relaxes to the correct statistical properties when the channel is submitted to equilibrium boundary conditions. The GCMC algorithm is a procedure allowing the simulation of a fluctuating number of particles in accord with Eq. 4 (Norman and Filinov, 1969; Allen and Tildesley, 1989). Briefly, it consists of constructing a random walk (discrete time Markov chain) of the configuration of the system during which particle creation and destruction can occur. To attempt the creation of an ion of type  in a system that contains n ions of that type, one inserts a new ion at a random position and calculates the probability

(5)

where  = V is the expectancy for the number of ions of type  in the region of volume V, and  [(... , n + 1, ...)  (... , n, ...)] is the change in free energy (PMF) of the system due to the new ion. The creation is accepted if a uniform random number between 0 and 1 is less than or equal to P_creat, otherwise it is rejected and the ion is removed from the system. To attempt destruction of an ion of type , one randomly chooses one ion out of the n particles, removes it from the system, and calculates the probability

(6)

where  [(... , n, ...)  (... , n  1, ...)] is the change in free energy (PMF) of the system due to particle removal. The destruction is accepted if a random number between 0 and 1 is less than or equal to P_destr, otherwise it is rejected and the ion is put back in its original position. To ensure microscopic reversibility, one must attempt particle creation and destruction with equal probability. Equations 5 and 6 are constructed from the constraint that the probability of creation/destruction obeys detailed balance (the formulas are derived in Appendix A). The GCMC procedure yields a statistical ensemble rigorously consistent with Eq. 4.

The GCMC/BD algorithm

To simulate the time course of the system realistically, one must carefully handle the unphysical particle creation and destruction inherent to the GCMC procedure by restricting it to "buffer regions" sufficiently remote from the main region of interest. Then, it may be expected that the GCMC/BD simulation algorithm can provide a useful approximation to study dynamical properties involving microscopic events far away from the boundary of the system.

The implementation of the GCMC/BD algorithm for simulating transmembrane ion channels is schematically illustrated in Fig. 1. A spherical geometry is shown for convenience, though other choices are also possible. The channel system is divided into five specific spatial regions: the inner region, the buffer regions on sides I and I, and the outer regions on sides I and II. The trajectory of all the ions in the inner and buffer regions is calculated with Eq. 2 according to the BD algorithm. The force acting on the ions is calculated from the derivative of the multi-ion PMF of the system taking into account the interactions between all the ions present in the inner and buffer regions and the influence of the channel, as well as the transmembrane potential. To ensure that constant conditions are maintained on the boundaries of the inner region, the ions in the buffer regions I and II on both sides of the membrane are kept in equilibrium with the bulk solution with which they are in contact. This is enforced via the GCMC procedure with particle creations and destructions. The creation and destruction of particles is implemented using the probabilities (Eqs. 5 and 6) where the chemical potential of ions of type  is specified by

(7)

It should be noted that the excess solvation energies µ(I) and µ(II) are influenced by ion-ion interactions in the bulk solution and are, thus, concentration-dependent (see below). In Eq. 7 it was assumed without loss of generality that the membrane potential is zero on the extreme left of the membrane (side I) and V_mp on the extreme right (side II).

View larger version (50K): [in this window] [in a new window]   FIGURE 1   Schematic representation of the ion channel-membrane system with the various regions for GCMC/BD simulations. The membrane potential is zero on the extreme left of the membrane (side I) and Vmp on the extreme right (side II). The channel system is divided into five specific spatial regions: the inner region, the buffer regions on sides I and II, and the outer regions on sides I and II. In the inner and buffer regions, the ions are treated explicitly and their trajectories are calculated according to BD with Eq. 2. In the buffer regions on sides I and II, the concentrations of the ions are maintained in equilibrium with their corresponding bulk solution using the GCMC procedure. In the outer region the ions are treated implicitly. The electrostatic potential is calculated using the modified Poisson-Boltzmann equation expressed in Eq. 10.

The GCMC procedure has the effect of enforcing boundary conditions corresponding to constant electrochemical potential in the two buffer regions. In contrast, no particle creation or destruction is taking place in the inner region, and the time course of the inner system is governed by BD. When the system is at equilibrium, the electrochemical potential of any ion is the same in all the regions of the system, and there is no net flow. However, when non-equilibrium conditions are imposed at the boundaries, a stationary state is simulated as particles flow from the regions with a high value of electrochemical potential to the regions with lower values. Because the buffer regions cannot run out of particles or be filled by particles, they essentially act as infinite thermodynamic reservoirs and sinks for the particles with respect to the central inner region. In practice, one performs several steps of GCMC in the buffer regions for each step of BD to prevent any ion accumulation or depletion in the buffers.

Microscopic model

The combined GCMC/BD algorithm described above is a general approach that allows the implementation of non-equilibrium boundary conditions of concentration and transmembrane potential across an ion channel. To have a practical simulation algorithm, we need to choose a particular microscopic model describing the effective ion-ion and ion-channel interactions. The multi-ion PMF must be computationally efficient because it is needed repeatedly during the GCMC/BD simulation for generating the BD trajectory according to Eqs. 2 and 3, and for calculating the creation and destruction probabilities according to Eqs. 5 and 6. Furthermore, the multi-ion PMF must also be carefully constructed because it is effectively the microscopic foundation of a given model.

For the sake of simplicity, we choose to represent the solvent as a structureless dielectric medium and incorporate its influence on the multi-ion PMF implicitly. We assume that the interaction between two ions separated by a distance r in the bulk solution is

(8)

where and are the parameters of the Lennard-Jones 6-12 (LJ) potential, q and q are the charges of the ions, and _bulk is the dielectric constant of bulk water. Equation 8 corresponds to the well-known restricted primitive model with a soft core that has been extensively used in statistical mechanical studies of ionic solutions (Ramanathan and Friedman, 1971; Ermak, 1975; Wood and Friedman, 1987; Jardat et al., 1999). A rigorous rationale for this type of approximation is provided by the MacMillan-Mayer theory of ionic solutions (McMillan and Mayer, 1945). Following this approximation for our microscopic model, the multi-ion PMF for a channel system is

(9)

where (r) is the electrostatic potential arising from the charge distribution of the channel and the transmembrane potential, and U_core is a repulsive potential preventing core-core overlap of the ions with the channel and membrane. The electrostatic potential (r) in Eq. 9 is calculated with the modified Poisson-Boltzmann (PB) equation (Roux, 1997, 1999)

Inner & buffer regions, sides I & II
Outer region, side I

(10)

Outer region, side II
where (r) and ²(r) are the space-dependent dielectric constant and Debye-Hückel screening factor, _p(r) is the charge density of the protein, and V_mp is the equilibrium electrode potential far away on side II. It may be noted that the space-dependent dielectric constant and Debye-Hückel screening factor are uniform within each region (e.g., inner, buffer, or outer regions) and material (e.g., protein, membrane, and bulk solution). The standard PB equation is recovered when V_mp = 0. Although the present approximation ignores the reaction field of the ions caused by the dielectric boundaries, it includes the dominant electrostatic effects. An efficient method for treating the reaction field in the GCMC/BD algorithm is currently in preparation. It should be noted that, in accord with our previous statistical mechanical analysis (see Eq. 26 of Roux (1999)), there is no screening factor ²(r) in Eq. 10 for the PB equation corresponding to the regions where the ions are represented explicitly (inner and buffer regions). This is in contrast with the treatment of Schirmer and Phale (1999), who included ionic screening at a mean-field level in the region of the aqueous pore using the PB equation. The modified PB Eq. 10 can be solved numerically using standard finite-difference methods (Warwicker and Watson, 1982; Klapper et al., 1986) and the result stored in memory for efficient computer simulations. The forces (first derivative of the PMF) can be calculated analytically from the grid-based potentials (Im et al., 1998). The core repulsive potential is also calculated on a discrete grid and stored for computational efficiency. Noting that its form is isomorphic to that of an ion of charge unity interacting with a positive electrostatic potential, the energy and forces are calculated with the same method. Further computational details about grid-based force calculations are given in Appendix B.

The excess chemical potential used in the GCMC creation and destruction probabilities is needed to implement correct GCE boundary conditions in the buffer regions. The value of µ(I) and µ(II) for ions of type  in Eq. 7 must be consistent with the microscopic model that was chosen (here, the primitive model). For computational convenience, the excess chemical potential of the cations and anions was calculated using the hypernetted chain (HNC) integral equation (Hansen and McDonald, 1976)

(11)

complemented by the Ornstein-Zernike equation

(12)

where the * represents a 3D spatial convolution and h(r) and c(r) are the ion-ion pair correlation function and direct correlation function, respectively. The concentration-dependent excess chemical potential based on the HNC equation is (Morita and Hiroike, 1960; Singer and Chandler, 1985)

(13)

The excess chemical potential of K⁺ and Cl for a primitive model of an aqueous KCl salt solution was calculated for concentrations varying from 10 mM to 1 M using standard computational techniques for solving the HNC equation (Roux et al., 1990). The LJ parameters of the K⁺ and Cl ions, given in Table 1, were taken from previous work (Roux, 1996). Standard combination rules were used, i.e., _ij =  and _ij = (_ii + _jj)/2. The calculated excess chemical potentials are given in Table 2.

View this table: [in this window] [in a new window]   TABLE 2   Excess chemical potential µ for the primitive model of KCl salt solution (kcal/mol)

	COMPUTATIONAL ILLUSTRATIONS

TOP ABSTRACT INTRODUCTION SIMULATION ALGORITHM COMPUTATIONAL ILLUSTRATIONS CONCLUDING DISCUSSION APPENDIX A APPENDIX B REFERENCES

Equilibrium structure of ionic solution

As a first illustration of the GCMC procedure, we consider the equilibrium structure of an aqueous 150 mM KCl salt solution. A spherical system of 50 Å radius with one Cl ion kept fixed at the origin was simulated. Using the values for 150 mM KCl given in Table 2, the excess chemical potentials are 0.18703 and 0.19008 kcal/mol for K⁺ and Cl, respectively. Because our aim is to test the ability of the GCMC algorithm to generate the correct equilibrium structure, the entire simulation system was treated as a buffer region (i.e., we allowed particles to be created or destroyed everywhere in the system). No BD steps were generated for this equilibrium system. After some equilibration, a GCMC simulation of 10⁴ cycles was generated; one cycle corresponds to 10 creation/destruction steps of GCMC and 100 steps of Metropolis Monte Carlo with random displacements of 1 Å in all three Cartesian directions.

In Fig. 2 (top) we compare the radial distribution function of K⁺ and Cl calculated from the GCMC/BD simulation relative to the fixed Cl ion with the results obtained from the HNC integral equation. There is an excellent agreement between the distribution functions, indicating that both the absolute number and relative positions of the ions are simulated correctly. The entire system is electrically neutral. The average number of K⁺ is 47.2 and the average number of Cl is 46.2 (excluding the fixed Cl ion). The fluctuation in the number of K⁺ and Cl is on the order of seven to eight particles, i.e., ~15% of the total number of each ion type. The magnitude of the fluctuations indicates that a fixed number of ions would provide a poor approximation to an open system obeying the statistics of the GCE.

View larger version (18K): [in this window] [in a new window]   FIGURE 2   Radial distribution function between a fixed Cl ion and the surrounding K+ and Cl ions of an aqueous salt solution. In the top plot, the radial distribution function calculated from the HNC equation is compared with the results of the GCMC simulation. In the bottom plot, the distribution function is shown for the entire spherical simulation system with 50 Å radius.

As shown in Fig. 2 (bottom), there is a slight decrease of the ion density near the boundary of the simulation system. The edge effect is caused by the neglect of the interactions with the ions localized in the outer region. In the present microscopic model of the PMF given by Eq. 9, we ignore the interactions of an ion approaching the edge of the simulation region with the other ions that are located in the outer region. In the present simulation the ions "see" a vacuum beyond the radius of 50 Å while there should be a bulk ionic solution extending to infinity. In a more elaborate microscopic model of the PMF, the influence of the salt solution in the outer region could be incorporated implicitly, as an electrostatic reaction field. Nevertheless, even without a reaction field the ion density is accurately simulated except within 3 Å from the edge of the system.

Transmembrane potential

As a second illustration of the GCMC procedure, we simulated an impermeable membrane surrounded by symmetric 150 mM KCl salt solution on both sides in the presence of a transmembrane potential of 1.0 V. A spherical system of 50 Å radius with an impermeable membrane of 20 Å thick was simulated. As in the previous test, the entire simulation region outside the membrane was treated as a buffer with the GCMC procedure and no BD steps were generated. After some equilibration, a GCMC simulation of 10⁶ cycles was generated: one cycle corresponding to 10 creation/destruction steps of GCMC and 100 steps of Metropolis Monte Carlo with random displacements of 1 Å in all three Cartesian directions. A stepwise repulsive potential of 200 kcal/mol was used to represent the repulsive membrane region.

The resulting density profile of the ions along the Z-axis perpendicular to the membrane plane is shown in Fig. 3. The interfacial polarization gives rise to a decrease of K⁺ ions on side I and an increase on side II. As expected, the opposite trend is observed in the case of Cl. The properties of the electric double layer are reproduced with a relatively small number of ions: there are ~32 K⁺ and Cl ions in the system on average, with fluctuations on the order of 7 to 8. There are 16.7 Cl on the left and 15.2 on the right, and there are 16.2 K⁺ on the right and 15.7 on the left.

View larger version (13K): [in this window] [in a new window]   FIGURE 3   Ion density in the electric double layer for a transmembrane potential of 1.0 V (the density is given in ion/Å3).

As in the case of a simple ionic solution, there is a slight decrease of the ion density near the edge of the simulation system because interactions with the ions localized in the outer region are not included. Nonetheless, the effect is small and localized far away from the membrane surface. The density profile was calculated by counting the average number of ions within slices of Z at a given position Z, and then normalizing by the cross-sectional area of the sphere at that position (the average density profile near Z = ±50 Å is determined with more statistical uncertainty because the circular cross-sectional area tends to zero at the edge of the system).

The GCMC procedure is implemented such that both sides I and II are in equilibrium with the salt solution with which they are in contact. The excess chemical potential of K⁺ and Cl, corresponding to the values for a 150 mM KCl solution given in Table 2, is the same on both sides of the membrane. In addition, the creation or destruction of an ion of type  in the buffer on side II is based on the chemical potential (II), which includes an offset of +qV_mp according to Eq. 7. Although there are no static charges present in this simple system (membrane or channel), the solution of Eq. 10 with _p(r) = 0 yields a non-zero electrostatic potential (r) because of the membrane potential V_mp. The potential , arising from the interfacial polarization of the ions located in the outer region, acts on all the ions in the inner and buffer region (side I and side II). On side II, the electrostatic potential arising from the ions in the outer region is nearly 1.0 V. On side I, the potential is around 0.0 V. The transmembrane potential is nearly linear across the membrane. The calculated potential is in fact very similar to an analytical solution for a planar membrane (Läuger et al., 1967; Roux, 1997) and this simple form was used in the present simulation for the sake of simplicity.

Particle diffusion with concentration gradient

So far we have only simulated systems in equilibrium with no ion flow. We now use the GCMC/BD algorithm to simulate non-equilibrium particle diffusion under a concentration gradient. For the sake of clarity and simplicity, a 50 Å × 50 Å × 50 Å cubic system without any particle-particle interactions is considered. A first buffer at a concentration of 500 mM is set in the region 25 < Z  20 Å on side I, and a second buffer at a concentration of 100 mM is set in the region 20 < Z < 25 Å on side II. Because it is more likely for a particle to enter the system on the side with the higher concentration, a stationary state with a net flux of particles is established in the +Z direction. A simulation of 10⁶ steps with a time step of 50 fs was performed; 10 steps of GCMC were performed for each step of BD to maintain the equilibrium in the buffer regions. The total length of the simulation is 50 ns.

The concentration profile calculated from the simulation is shown in Fig. 4. As expected, the simulated concentration profile is linear and corresponds exactly to that of a simple 1D diffusion problem along the Z-axis. The concentration gradient is established over a distance d of 40 Å, between 20 and +20 Å (the thickness of each buffer region is 5 Å). The linear concentration profile is established despite the relatively small number of particles involved in the simulation system. The concentrations of 500 mM and 100 mM correspond to 3.0 × 10⁴ part/ų and 6.0 × 10⁵ part/ų, respectively. This yields an average of 3.76 and 0.75 particles in the two buffers. The number of particles in each buffer varies according to the GCE as shown in Fig. 5. Because there are no particle-particle interactions, the number of particles in each buffer follows exactly a Poisson distribution with P_n = ⁿ/n! exp[n/].

View larger version (12K): [in this window] [in a new window]   FIGURE 4   Concentration profile for gradient diffusion (in part/Å3). On side I a buffer is set at a concentration of 500 mM (3.0 × 104 part/Å3). On side II a buffer is set at a concentration of 100 mM (6.0 × 105 part/Å3). The gradient is established over a distance of 40 Å.

View larger version (14K): [in this window] [in a new window]   FIGURE 5   Probability distribution of the number of particles in the buffer regions I and II.

The net flux of particles calculated from the GCMC/BD simulation is 0.0029 part/ps. The mean-square deviation is estimated to be on the order of 0.0003 part/ps based on nine independent runs. According to Fick's law (Fick, 1855; Hille, 1992), the net flux of particles is given by

(14)

where A is the XY cross-sectional of the system (2500 Å²). The particle flux calculated from Eq. 14 with a diffusion constant of 0.196 Å²/ps is 0.003 part/ps, in excellent accord with the result from the GCMC/BD simulation. The fluctuation in the number of particles is essential to obtain the correct (non-integer) average number of particles in the buffer regions (3.76 and 0.75), and thus the correct boundary conditions of concentration for the gradient diffusion.

Ion flux through the OmpF porin

As a last illustration of the GCMC/BD algorithm, we now consider the non-equilibrium flow of K⁺ and Cl ions through the porin OmpF of E. coli. The OmpF channel is well-characterized, both structurally and functionally, which makes it an excellent system for testing computational methods (Schirmer, 1998). The 3D structure of OmpF has been determined by x-ray crystallography (Cowan et al., 1992); each OmpF monomer of the trimeric structure possesses a wide ion-conducting aqueous pore. The conductance of the monomers has been measured for a range of salt concentrations (Prilipov et al., 1998).

An atomic system was constructed to perform GCMC/BD simulation with an orthorhombic region (41.6 Å × 41.6 Å × 86.0 Å) corresponding roughly to one OmpF monomer. Buffers of 1 M KCl salt solutions were implemented from 43 Å to 40 Å and from +40 Å to +43 Å using the GCMC procedure; the excess chemical potentials were taken from Table 2; they are 0.309 and 0.334 kcal/mol for K⁺ and Cl, respectively. The electrostatic potential arising from the protein and the transmembrane voltage was calculated with the modified PB Eq. 10 using a fully detailed atomic model of OmpF for transmembrane potentials of +100 and 100 mV. To set up the atomic system, the entire trimer with its symmetry axis oriented along the Z-axis was embedded in an explicit lipid bilayer using the protocol of Woolf and Roux (1996). The default state was used for all ionizable residues except for Glu-296 and Asp-312, which were protonated as in previous work (Karshikoff et al., 1994; Schirmer and Phale, 1999). The atomic charges were taken from the all-atom PARAM22 force field (MacKerell et al., 1998) of the CHARMM program (Brooks et al., 1983). The net charge of the OmpF trimer is 30e. A dielectric constant of 2 was used for the interior of the protein and membrane regions. A dielectric constant of 80 was assumed for the bulk solvent region (including the aqueous pore region of OmpF). A membrane thickness of 24 Å was used. The optimized atomic radii for proteins were used to set up the dielectric boundaries (Nina et al., 1997); CHARMM radii were used for the lipid molecules (MacKerell et al., 1998). The radii of the hydrocarbon atoms at the center of the bilayer were increased to fill any remaining empty cavities in the hydrocarbon core of the membrane. The electrostatic potential was first calculated with a coarse grid (131³ points with a grid spacing of 1.0 Å) centered on the entire OmpF trimer with periodic boundary conditions imposed in the X and Y directions. The result of the coarse calculation was then used to set the potential on the edge of the box to perform a second calculation using a finer grid (107 × 107 × 217 with a grid spacing of 0.4 Å) centered on a single OmpF monomer. A Debye-Hückel screening factor ²(r) corresponding to a 1 M KCl salt concentration is assumed in the outer region. An ion exclusion Stern layer of 1.5 Å in addition to the protein atomic radii was used to set the spatial dependence of the ionic screening factor. However, in contrast to the BD simulations of Schirmer and Phale (1999), ²(r) is set to zero in the inner and buffer regions because the explicit ions are included explicitly (see Eq. 10). All continuum electrostatic calculations were performed using the PBEQ (Poisson-Boltzmann Equation) module (Roux, 1997; Im et al., 1998) incorporated in the biomolecular simulation program CHARMM (Brooks et al., 1983).

The GCMC/BD simulations were generated with a time step of 25 fs; 10 steps of GCMC were performed for each step of BD. The electrostatic energy and forces were calculated using a 2nd-order B-spline interpolation scheme according to the method described in Appendix B. The core repulsion potential U_core was calculated and stored on a grid similar to the electrostatic potential; a 3rd-order B-spline interpolation scheme was used to calculate the repulsive energy and forces (see Appendix B). A single uniform repulsive potential of +20 kcal/mol based on the LJ radius of K⁺ was used for both ions. The molecular surface (including contact and reentrant surface) was used to ensure that no region in the protein interior is accessible to the ions in the grid-based repulsive core potential. To obtain statistical convergence, 10 independent GCMC/BD simulations were generated with different initial configuration and seed numbers for the random number generator. The total length of each simulation is 0.5 µs. One GCMC/BD simulation took ~40 h on a single Pentium II 400 mHz processor. The results of the GCMC/BD simulations of OmpF averaged over the 10 runs are summarized in Table 3.

Fig. 6 shows a typical snapshot of the OmpF system for GCMC/BD simulations. Several K⁺ and Cl are observed in the channel and vestibule regions. More specifically, one K⁺ and one Cl are in the pore near the constriction zone formed by the L3 loop (Schirmer, 1998; Schirmer and Phale, 1999). The Cl ion is close to the so-called "arginine cluster" (Arg-42, Arg-82, and Arg-132) while the K⁺ ion is close to the acidic residue Glu-117 (Cowan et al., 1992; Schirmer, 1998). On average, the entire GCMC/BD system contains ~86 ions (45 K⁺ and 41 Cl) for V_mp = 100 mV and 88 ions (46 K⁺ and 42 Cl) for V_mp = +100 mV. The interior of the pore was defined as the part of the aqueous channel corresponding to the -barrel of OmpF. According to this definition, the OmpF channel nearly spans 32 Å along the Z-axis. The wide vestibule formed by the long extracellular loops was not included in the definition of the channel because it is extensively exposed to the bulk solution. The average number of K⁺ ions inside the pore (6.7 and 7.1) is significantly larger than the number of Cl ions (1.8 and 1.6) reflecting the specificity of OmpF for cations (Cowan et al., 1992). The small difference in the number of K⁺ and Cl for +100 mV and 100 mV is caused by the asymmetric charge distribution and shape of the OmpF pore. It may be noted that the difference in the averages is much larger than the fluctuations (see Table 3).

View larger version (68K): [in this window] [in a new window]   FIGURE 6   Molecular graphics view of the GCMC/BD simulations of the OmpF pore. The orthorhombic simulation box is highlighted in green. The 3 Å thick regions corresponding to the 1 M KCl buffers are shown. K+ and Cl ions are observed in the pore and in the vestibule of the monomer. A pair of K+ and Cl ions is in the pore near the constriction zone formed by L3. The Cl ion is close to Arg-42, Arg-82, and Arg-132, and the K+ ion is close to the acidic residue Glu-117 (ball-and-stick model representation).

Fig. 7 shows the number of K⁺ and Cl ions crossing the plane Z = 0 in the +Z direction (i.e., from the bottom buffer to the top one) for V_mp = ±100 mV. Each figure shows parts of two different GCMC/BD trajectories. As expected, the total number of ions flowing in the direction of the transmembrane electric field is mostly increasing with time, although there are some fluctuations. The observed flow of ions is not symmetric with +100 and 100 mV. The calculated conductance is 1.88 nS for +100 mV and 2.40 nS for 100 mV. This can also be observed in Fig. 7 where the total number of K⁺ at t = 150 ns is 145 in +100 mV and 180 with 100 mV. The asymmetry of the ion flux in response to an applied potential reflects the charge distribution and shape of the OmpF pore. Interestingly, the difference between the partial average occupancy numbers of K⁺ and Cl for +100 mV and 100 mV is smaller than the corresponding partial fluxes. This may reflect the fact that probability of occupancy is determined by the entire volume of the aqueous pore, whereas the rate of transport is governed by the energy barrier in the constriction zone.

View larger version (17K): [in this window] [in a new window]   FIGURE 7   Number of ions crossing the plane at Z = 0 in the center of the OmpF monomer as a function of time.

The calculated conductance, ~2.0 nS, is larger than the experimental value of 0.84 nS determined at the same salt concentration (Prilipov et al., 1998). The disagreement with the experimental value, which is not surprising given the simplicity of the present model, may be due to a number of factors. In particular, bulk values were used for the diffusion constant of K⁺ and Cl. This is certainly an overestimate. Several computational studies suggest that the diffusion constant of ions and water molecules is reduced inside a molecular pore (Sansom et al., 1996; Roux and Karplus, 1991a; McGill and Schumaker, 1996; Chiu et al., 1993; Tieleman and Berendsen, 1998). In addition, the present treatment of GCMC/BD based on Eq. 9 ignores the effect of the reaction field of the ions caused by the dielectric boundaries and the ions present in the outer region. The reaction field may have complex effects. However, it should on average make it more unfavorable for the ions to leave the high dielectric bulk region to enter the narrow pore. Therefore, it seems likely that the influence of the reaction field should tend to reduce the calculated conductance. Similarly, we expect that the average number of ions inside the pore would be reduced if the influence of the reaction field were included.

	CONCLUDING DISCUSSION

TOP ABSTRACT INTRODUCTION SIMULATION ALGORITHM COMPUTATIONAL ILLUSTRATIONS CONCLUDING DISCUSSION APPENDIX A APPENDIX B REFERENCES

A novel algorithm based on GCMC and BD for simulating ion channels has been described and implemented. The combined GCMC/BD algorithm described above is a general approach that allows the implementation of arbitrary non-equilibrium boundary conditions of concentration and potential for studying ion fluxes across a membrane channel system. Although the algorithm was developed and illustrated with BD, dynamical inertial and memory effects could be incorporated by using the Langevin equation (LE) or generalized Langevin equation (GLE), respectively (McQuarrie, 1976).

We have illustrated the GCMC/BD algorithm with a few simple cases such as an equilibrium salt solution, an equilibrium ionic double layer under a transmembrane potential, a non-equilibrium concentration gradient, and lastly, non-equilibrium ion flux through a detailed atomic model of porin OmpF of E. coli. The approach will allow routine simulations of detailed models of molecular pores with various non-equilibrium conditions of salt concentration and transmembrane potential for comparison with experimental data.

In the present simulations we used a continuum electrostatic approximation related to the restricted primitive model of ionic solutions (Ramanathan and Friedman, 1971). While this is certainly an approximation, it incorporates the main features of ion-ion and ion-channel interactions that govern permeation process. However, the multi-ion PMF appears as an input in the GCMC/BD simulation framework and, therefore, more sophisticated choices are possible. For example, the PMF could, in principle, be extracted from detailed MD of the channel with explicit solvent and membrane (Roux and Karplus, 1993; Woolf and Roux, 1997). Although this is arguably the best route for studying ion permeation at the microscopic level, there are a number of difficulties with this approach. In particular, the evaluation of multi-ion PMF using biased sampling and free energy methods requires extensive calculations. Furthermore, the final results may be disappointing. Although the PMF calculated by MD simulations is sufficiently accurate to provide the salient features about the position of barriers and wells, it is often not sufficiently accurate to yield a quantitative description of ion fluxes and selectivity because ion permeation is extremely sensitive to small variations in free energies. For these reasons, other approximations for the multi-ion PMF have, and will continue, to play an important role, e.g., mixed continuum-discrete models (Dorman et al., 1996) or simplified atomic models (Roux and Karplus, 1991b; Chung et al., 1999). In particular, adjustments and optimization of phenomenological free energy profiles (based on MD results) to fit experimental data such as performed by McGill and Schumaker (1996) will be an important approach to examine the sensitivity of the transport properties upon the detail of a model. Future challenges will consist of combining information from MD, statistical mechanical theories of ionic solutions, continuum electrostatics, and experimental data to construct effective and accurate approximations to the multi-ion PMF. It is our hope that the GCMC/BD algorithm will serve as a rigorous framework for testing these approximations and progressing in our understanding of ion channels.