INTRODUCTION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The goal of elucidating the structure-function relationships in biological ion channels has gained a new impetus with the determination of the KcsA potassium channel structure (Doyle et al., 1998). Most of the theoretical efforts in modeling the KcsA channel have so far focused on molecular dynamics (MD) simulations of potassium ions in the channel (Guidoni et al., 1999, 2000; Allen et al., 1999, 2000; Shrivastava and Sansom, 2000; Åqvist and Luzhkov, 2000; Luzhkov and Åqvist, 2000; Berneche and Roux, 2000; Roux et al., 2000). These studies provide valuable information on the selectivity mechanism and the energetics of ion permeation in the channel, but do not make predictions about the quantity that can be directly measured by experiment, namely the conductance. In a recent 100-ns MD simulation, Crozier et al. (2001) calculated the conductance of a simplified channel in somewhat extreme conditions (1 M solution with a 1.1-V applied potential). This gives hope that it may be possible to determine conductance of biological channels from MD studies under physiological conditions in the not too distant future. Currently, however, typical MD simulations of biological channels can be run for ~10 ns, which is too short to estimate the channel conductance, or even to explore the dynamics of a single-conduction event. Of course, this is not a new problem, and permeation models of lower resolution such as Brownian dynamics (BD) (Cooper et al., 1985; Kuyucak et al., 2001) and Poisson-Nernst-Planck equations (Levitt, 1986; Eisenberg, 1999) have long been considered in the literature. The latter approach has recently been shown to be invalid in a narrow pore environment because it neglects the self-energy of ions (Corry et al., 2000). There has also been some initial development of a "solvent primitive model" in infinite cylinders with no dielectric boundaries, in which uncharged solvent molecules are explicitly simulated (Tang et al., 2001). This technique has not been applied to biological channels so far, and, as such, its accuracy is yet to be determined. This leaves BD simulations as a computationally tractable tool for calculating a channel's conductance from its structure.

BD simulations were first proposed as a way to study ion channels by Cooper et al. (1985). The early simulations involved one-dimensional (1D) studies of schematic channels (Jakobsson and Chiu, 1987; Bek and Jakobsson, 1994), but the extension to three-dimensional (3D), necessary for realistic modeling, has not been achieved until recently. The difficulty lies in the calculation of the forces on ions at each time step, typically found from the solution of Poisson's equation, which is computationally too expensive if done numerically (Hoyles et al., 1998a). Thus, the first 3D BD simulations were performed by Li et al. (1998) using a torus-shaped channel, for which analytical solutions of Poisson's equation are available (Kuyucak et al., 1998). For a channel with an arbitrary shape, this problem was finally resolved by storing the potential and electric field values in a set of lookup tables, and interpolating the required values during simulations from the table entries (Hoyles et al., 1998b). Since then, 3D BD simulations have been used in model studies of acetylcholine receptor (Chung et al., 1998), KcsA potassium (Chung et al., 1999, 2002; Allen et al., 2001), L-type calcium (Corry et al., 2001), and Porin channels (Schirmer and Phale, 1999; Im et al., 2000; Phale et al., 2001).

Recently, questions have arisen about the methods of implementing the boundaries in BD simulations of ion channels. Another problem, distinct from the issue of the boundaries, is that of accurate representation of the forces on ions in the channel. In our simulations, we have concentrated on representing the forces acting on ions accurately because, ultimately, they are responsible for driving the ions across the channel. Also, the calculated conductance values could be very sensitive to errors in electric fields and potentials, e.g., conductance has an exponential dependence on energy barriers in channels. In contrast, we regard the implementation of reservoir boundaries, required to maintain ion concentrations and create driving potentials, as a secondary issue. The purpose of the reservoirs is to move the necessarily unphysical system boundaries away from the critical part of the simulation. Provided the reservoirs are large enough, a simple implementation of the boundaries should then suffice. We implement the boundaries by applying a uniform electric field across the channel and keeping a fixed number of ions in the reservoirs. The chosen concentration values in the reservoirs are maintained by recycling ions from one side to the other whenever there is an imbalance due to a conduction event: this process mimics the current flow through a closed circuit. There has, however, been a great deal of debate in the field about the appropriateness of such a simple stochastic boundary. Most recently, Im et al. (2000) have proposed a more elaborate treatment of boundaries using a grand canonical Monte Carlo (GCMC) method. In this paper, they also question the validity of the simple method, but unfortunately do not support this criticism with any hard evidence, such as a comparison of the two methods. A separate issue is that their treatment of forces on the ions is less sophisticated than in our simulations, in that the self-energy of the ions (or the reaction field) is ignored. Although the authors recognize that this is an important deficiency and plan to improve their technique, the work of Im et al. (2000) has been deemed to be the most rigorous so far for BD simulations of ion channels (e.g., Phale et al., 2001; Tobias, 2001; Tieleman, 2001).

The source of this preoccupation with boundaries in BD appears to arise from association with MD simulations where the correct treatment is known to be crucial (Sagui and Darden, 1999). However, the nature of the electrostatic forces in BD is very different from those in MD: first, an ion's electric field (or potential) in water is reduced by a factor of 80 due to the dielectric shielding, and second, shielding due to the counterions completely annuls the remaining field strength beyond 4 Debye lengths. For physiological concentrations (150 mM), this length scale is about 30 Å. With this provision on the reservoir dimensions, we believe that a simple boundary method should be adequate for the purpose of calculating the conductance of a channel from BD.

In view of the debate outlined above, however, it seems prudent to perform additional tests on the validity of our simple stochastic boundaries. The work of Im et al. (2000) provides us with an opportunity to do so. We have modified our computer programs to deal with the more sophisticated boundaries, and have carried out BD simulations of model channels using the two different methods. We have also experimented with different methods of representing the transmembrane potential. The purpose of these tests is to determine whether the GCMC boundaries (or the alternative representations of the membrane potential) make any difference to the conductance of the model channel, or to the concentrations of ions near the mouths of the channel. If not, we feel safe in concluding that the reservoirs are adequately insulating the channel from the boundary conditions, and that our simulations accurately reflect the physical processes taking place in ion channels.

	METHODS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Brownian dynamics

Because the application of BD simulations to realistic 3D channel geometries has been discussed in detail in our earlier papers (e.g., Li et al., 1998; Hoyles et al., 1998b; Chung et al., 1998, 1999; Corry et al., 2001), we give only a brief description of the method here and focus on the boundary methods that are to be compared. In BD, the motion of individual ions are simulated using the Langevin equation,

(1)

where m_i, q_i, and v_i are the mass, charge, and velocity of the ith ion. In Eq. 1, the effect of the surrounding water molecules is represented by an average frictional force with a friction coefficient m_ii, and a stochastic force F_R arising from random collisions. The last two terms in Eq. 1 are, respectively, the electric and short range forces acting on the ion.

The total electric field at the position of the ion is determined from solution of Poisson's equation, and includes all possible sources due to other ions, fixed and induced surface charges at the channel boundary, and the applied membrane potential. For the proposed channel boundaries used in this study, Poisson's equation can only be solved numerically. This is achieved using the boundary charge method (Levitt, 1978), which is improved by including the effect of curvature of sectors in the solutions (Hoyles et al., 1998a). Because solving Poisson's equation at each time step is computationally prohibitive, we store precalculated values of the electric field and potential due to one- and two-ion configurations in a system of lookup tables, and interpolate values from these during simulations (Hoyles et al., 1998b). The short-range forces are used to keep the ions in the system and also to mimic other interactions between two ions that are not included in the simple Coulomb interaction. These include short range repulsion and hydration effects as described previously (Corry et al., 2001).

The Langevin equation (Eq. 1) is solved at discrete time steps following the algorithm devised by van Gunsteren and Berendsen (1982). To simulate the short-range forces more accurately, we use a multiple time-step algorithm in our BD code. A shorter time step of 2 fs is used across the channel where short range ion-ion and ion-protein forces have the most impact on ion trajectories, whereas, elsewhere, a longer time step of 100 fs is used.

Simulations under various conditions, each lasting for one million time steps (0.1 µs), are repeated numerous times. Initially, a fixed number of ions are assigned random positions in the reservoirs with velocities also assigned randomly according to the Boltzmann distribution. The cylindrical reservoirs have a fixed radius of 30 Å, and their height is adjusted to obtain the desired concentration. The concentrations in each of the reservoirs are maintained using one of two stochastic boundary techniques described below. The current is determined from the number of ions crossing the channel during the simulation period.

The BD program is written in FORTRAN, vectorized, and executed on a supercomputer (Compaq AlphaServer SC). The time to complete the simulations depends on the number of ions, how often ions enter the short time-step regions, and whether the simple or GCMC boundaries are used. A temperature of 298 K is assumed throughout and a list of the other parameters used in the BD simulations is given in Table 1. (Note that the diffusion coefficient is related to the friction coefficient m in Eq. 1 by the Einstein relation, D = m/kT.)

Reservoir boundaries

To maintain the specified concentrations in the reservoirs, we apply stochastic boundaries. Here, we compare the use of a simple boundary that maintains a fixed number of ions in the system with a more sophisticated GCMC boundary that allows fluctuations in the number of ions.

Simple stochastic boundary

Grand canonical Monte Carlo method

(2)

(3)

(4)

View larger version (16K): [in this window] [in a new window]   FIGURE 1   Diagram of the channel and reservoir system used in BD/GCMC simulations. The protein and membrane forms a 3D channel when rotated about the central axis by 180°. In this case the channel is cylindrical, although it may have any shape. Attached to each end of the channel are cylindrical reservoirs. During BD simulations, ions move within this channel and reservoir system. When the GCMC procedure is used to maintain ion concentrations in the reservoirs, ions are created and destroyed in the buffer zones around the outside edge of the reservoir indicated by the dashed lines. The dimensions shown are those for the cylindrical channel discussed in the text.

(5)

Transmembrane potentials

A second issue to do with reservoir boundaries is how to apply a potential difference that drives ions through the channel. There are at least three possibilities that have been considered in the past, and there has been much debate as to which is most appropriate.

In all our recent BD simulations, we have created a transmembrane potential by simply applying a uniform electric field through the system. This applied field is included in the solution of Poisson's equation with the dielectric boundaries so that it induces surface charges of it own. The resulting potential is far from being linear across the systemit drops much more rapidly through the channel than in the reservoirs.

Another approach is to fix the potential at the desired values along the far ends of the reservoirs. This creates an equipotential surface at each end, which can be set independently to create a potential drop across the channel. This is similar to placing electrodes at the far end of each reservoir, though, in an actual experiment, the electrodes would be much farther from the channel than in a typical simulation. To use such a scheme, we solve Poisson's equation with the specified boundary potentials using a finite difference method (Moy et al., 2000). The results due to the transmembrane potential, fixed charges in the protein wall, and charges induced by these are stored in a 3D lookup table.

A final method, which has been introduced by Im et al. (2000), is to make the potential more realistic by moving the fixed potential surfaces far away from the simulation system. The electrolyte solution in between the reservoir and the fixed potential surface is treated as a continuum by solving the Poisson-Boltzmann (PB) equation in this region. The potential in the system is thus calculated using a modified PB equation, which reduces to Poisson's equation in the BD simulation system where ions are treated explicitly. We implement this procedure by using a finite difference method to solve the modified PB equation when constructing the 3D lookup table.

For the purpose of comparing the stochastic boundary techniques used for maintaining concentrations in the reservoirs, it is important that we use the same applied potentials. Thus, in all the simulations discussed here, except for the final section on transmembrane potentials, we utilize the first "uniform field" approach.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

Equilibrium ion distribution

We first demonstrate that the BD simulations with either the simple or GCMC procedure maintains the desired equilibrium conditions by examining the relative distribution of ions in bulk solution. For this purpose, we set all dielectric constants in the system equal to 80 so that there are no dielectric boundaries in the system, and ignore ions that are within 8 Å of the reservoir boundaries to avoid edge effects in sampling. In Fig. 2, we show the radial distribution functions for K-Cl (A) and K-K (B) ion pairs, obtained from BD simulations of a 500-mM KCl solution for 10⁶ time steps (0.1 µs), in one case with the GCMC routine in place (triangles) and in another without (filled circles). When testing the GCMC procedure, the buffer regions are enlarged to occupy the entire reservoirs (so that ions can be created or destroyed anywhere) because the test is for a bulk solution. The curves agree closely and depict the peaks due to the contact and solvent-separated minima in the potential of mean force. They also closely follow the results obtained from the hypernetted chain equations indicated by the solid line. This agreement indicates that the equilibrium structure of the electrolyte is accurately reproduced in BD simulations with or without the GCMC routine.

View larger version (12K): [in this window] [in a new window]   FIGURE 2   Radial distribution function for (A) K+-Cl and (B) K+-K+ ion pairs moving in the reservoirs as found from BD with the GCMC routine in place throughout the reservoirs (triangles), BD simulations without the GCMC routine (filled circles), and the HNC equations (solid line). The BD plots were calculated by sampling from simulations using a cylindrical channel with no dielectric boundaries and ignoring ions within 8 Å of the reservoir boundaries.

Cylindrical channels

Because we are interested in comparing different treatments of the boundary in BD simulations, it matters little which channel model we use. Thus, for simplicity, we first make our comparisons in a simple cylindrical channel, before demonstrating the robustness of these results in the more complex potassium channel model we have studied previously.

The cylindrical channel model is formed by rotating the curve shown in Fig. 1 about the central axis. The channel radius is set to 3 Å and its entire length to 35 Å. First, we set the dielectric constant of the protein to 80, equal to that of the electrolyte in the channel and reservoirs. In this case, because there is no dielectric boundary, no reaction field can be induced. Although not realistic, this provides the simplest case in which to test the stochastic boundary methods, and it also helps in showing the importance of the reaction field when these results are compared to those with a realistic choice of dielectric constant. The overall concentration in the simple boundary simulations is set to be the same as the average concentration during the GCMC simulations. For compatibility with earlier simulations, all studies in the cylindrical channels are carried out using NaCl solution.

Ion distributions and fluctuations

(6)

View larger version (17K): [in this window] [in a new window]   FIGURE 3   Fluctuations in cation numbers in a region around the channel mouth comprising 50% of the volume of one reservoir when (A) the simple stochastic boundary is used and (B) when the GCMC procedure is used. A dielectric constant of 80 is used everywhere. The number of ions in the region is sampled every 100 BD steps and the relative frequency is calculated during a 0.2-µs simulation period. The average concentration in the simulation is 280 mM, corresponding to 15 ions of each type in each reservoir. The average number of ions in the region and the standard deviations are indicated by and , respectively. The dashed line shows the corresponding binomial distribution from Eq. 6 with N = 15 and x = 0.5.

View larger version (40K): [in this window] [in a new window]   FIGURE 4   Fluctuations in cation numbers as in Fig. 3 except in a smaller region near the channel mouth containing 25% of the entire reservoir volume using (A) the simple boundary and (B) the GCMC boundary.

Channel currents

View larger version (13K): [in this window] [in a new window]   FIGURE 5   Comparison of the current-voltage curves in a cylindrical channel using the simple stochastic boundary (filled circles) or the GCMC boundary (triangles). A 3-Å-radius channel is used with  = 80 everywhere and 265 mM NaCl solution. Ions are driven across the channel by an electric field of 2 × 107 V/m, corresponding to a potential drop of approximately 200 mV across the system. Sodium and chloride currents are plotted separately and each set of results is fitted by the solid line.

View larger version (11K): [in this window] [in a new window]   FIGURE 6   Comparison of the current-voltage curves in a cylindrical channel with fixed charges using the simple stochastic boundary (filled circles) or the GCMC boundary (triangles). A 3-Å-radius channel is used with eight monopoles placed at each end as described in the text. A dielectric constant of 2 is assigned to the protein and 80 everywhere else.

More complex channels

We have seen that the choice of stochastic boundary used to maintain concentrations in the simulation reservoirs has no effect on the currents flowing through simplified cylindrical channels. As a final test, we check to make sure that this conclusion is valid in a more complex and realistic multi-ion channel that we have modeled previously, and checking it at a range of concentrations. To do this, we use the KcsA potassium-channel model that has been described in an earlier paper (Allen and Chung 2001). An open-state KcsA-channel shape has been constructed using MD from the known closed-state crystal structure (Allen et al., 2000). A dielectric interface is then constructed by tracing out a boundary using a water molecule and assigning the dielectric constant a value of 2 in the protein and 60 in the channel. The final shape, and the pore-forming peptide helices are shown at the top of Fig. 7. Charges are assigned at positions corresponding to the protein atoms using the extended CHARMM-19 parameter set. More details can be found in the above references.

View larger version (44K): [in this window] [in a new window]   FIGURE 7   Conductance concentration curve in a model KcsA potassium channel found from BD simulations using the simple stochastic boundary (filled circles) and GCMC boundary (triangles). An electric field of 2 × 107 V/m is used to drive the ions across the channel. The results represent the averages of 1.5-2.0-µs simulations. The inset shows the shape of the channel and indicates the major peptide helices.

In Fig. 7, we plot the current-concentration curve in the KcsA potassium channel surrounded by KCl solution under a driving potential of 200 mV. The results of our simulations show a saturation of channel current with increasing conductance, and are fitted by a Michaelis-Menton curve to indicate this. The data from simulations carried out using the simple stochastic boundary, indicated by the filled circles, are those published elsewhere (Allen and Chung, 2001). When the GCMC boundary is used, the new data shown by the triangles reproduces this curve well at all concentrations studied. Thus, once again, the choice of stochastic boundary has little effect on channel currents, even over a large range of concentrations.

Transmembrane potentials

So far, we have examined techniques for maintaining ion concentrations in the reservoirs. A separate but related issue arises when we consider applying a transmembrane potential to drive ions through the channel, because this may rely on setting boundary values for the potential at the edges of the reservoir. In reality, the driving field or potential arises from ion clouds on each side of the membrane or a distant electrode. Thus, setting the potential at the back edges of the reservoirs is not quite correct because it fails to allow for variations at these positions caused by the movement of mobile ions. Im et al. (2000) claim that their use of the modified PB equation avoids these difficulties by taking into account the effects of mobile ions when determining the potentials at each end of the BD simulation. The uniform field approach does not include the cause of the transmembrane potential, rather just takes it as given, a bias that could be created by ionic clouds or polarization external to the BD system.

But, rather than entering a debate as to which is the most realistic way to create the transmembrane potential, we simply demonstrate here that it again makes no difference which method is used, by directly comparing the three. In Fig. 8, we plot the average potential along the central axis of the channel during a BD simulation using the uniform field, fixed potential, and modified PB approaches. In all cases, the majority of the potential drop occurs in the channel, with the potential remaining relatively flat in the reservoirs. When the fixed potential or modified PB methods are used, a slightly greater charge separation occurs in the reservoirs due to ions being attracted to the potential generating electrodes. This is especially true near the outside edge of the reservoir, and leads to the drop in the magnitude of the potential there. The boundary effects created in these techniques are, however, quickly screened out by the mobile ions in the system. It is worth noting that, in the modified PB method, we solve the modified PB equation only once before doing the BD simulation. Thus, the electrolyte outside of the reservoirs does not react to the presence of the explicit ions in the BD simulation. If the modified PB equation was solved at each BD time step, it is possible that the electrolyte would act, on average, to cancel some of the charge separation seen in the BD simulation. This, of course, would only act to bring the potential closer to the uniform field case, and would not alter our conclusion. Inside the channel, the potentials are almost identical in all cases, and good agreement is maintained until ~20 Å from the channel mouth. Thus, no matter which technique is used to create the transmembrane potential, the average potential seen by ions in or near the channel is the same. Because it is the potentials in and around the channel that drive ions through it, the choice of technique for creating a transmembrane potential is, therefore, irrelevant when it comes to calculating the current passing through the channel in a BD simulation.

View larger version (12K): [in this window] [in a new window]   FIGURE 8   Average potential profiles along the channel axis during a BD simulation when the transmembrane potential is set via the uniform field (solid line), fixed potential (dashed line) and modified PB (dotted line) methods. The cylindrical channel model is used, with a dielectric boundary but no fixed charges. The potential is plotted between the two ends of the reservoirs. In all cases, the potential at each point is averaged over a 0.1-µs BD simulation.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

We have presented a number of results that demonstrate that it does not matter whether the simple stochastic or the GCMC stochastic boundary is used to maintain ion concentrations in the reservoirs during BD simulations. In both cases, the edge of the reservoirs or the GCMC buffer zones must be at a reasonable distance from the channel, ideally 3-4 Debye lengths, such that any unphysical edge effects or particle creation or destruction are screened from the channel. When these precautions are followed, the number of ions near the channel fluctuates according to the binomial distribution, and the current passing through the channel is the same with either method. Similarly, it does not matter how the transmembrane potential is set. Provided the reservoirs are large enough, mobile ions redistribute themselves, causing the potential drop across the channel to be the same in all cases.

The simple boundary method is conceptually simpler, involves less calculations, and is considerably faster. For a typical simulation presented here with a 300-mM solution, a 1-µs simulation takes ~45 h of CPU time to complete using the simple boundary, and 115 h with the GCMC boundary method. The greater time is due to the potential energy of the system having to be recalculated for each GCMC creation or destruction step. The simple boundary also allows one to specify beforehand the exact concentration that will be used in a given simulation. Thus, for BD simulations of solutions at the usual physiological concentrations, the added complexity of the GCMC method provides no perceptible advantages compared to the simple boundary method.

One situation where the GCMC boundary does have an advantage over the simple boundary method is the simulation of solutions at very low concentrations. For example, to simulate an ion species in the micromolar range using the simple boundaries would require reservoirs thousands of times larger than those described here to contain at least one ion of this type. This is not only cumbersome, but also makes including a second ion species at a higher concentration (say in the millimolar range) problematicthe reservoir would have to contain thousands of ions of the second type, making it too slow to simulate. The GCMC boundary can reach such low concentrations using a small reservoir because there need not always be an ion of each type in the simulation. The low concentration species is simply created in the buffer regions at a proportionally lower rate.

Of course, it is possible to treat the boundaries in other ways not discussed here. For example, a periodic boundary could be used to maintain ion concentrations and potentials, such as that typically used in nonequilibrium molecular dynamics simulations (Crozier et al. 2001; Tang et al. 2001). These techniques have their own advantages, such as avoiding any explicit potential boundary, and disadvantages such as only being able to model symmetric solutions at each end of the channel. However, from what we have seen here, it should be apparent that the choice of boundary does not matter, provided some common-sense precautions are observed.

Although the GCMC boundary opens up some new avenues for simulation at low concentrations, it is not, despite the claims of Im et al. (2000), a more accurate method than the simple stochastic boundary. The results presented here support the obvious expectation that the physics taking place near the channel is the main determinant of channel currents, not the way the boundary conditions are implemented. Thus, as long as the reservoirs are large enough so that the edge effects are completely screened out near the channel, one need not worry about the exact implementation of the system boundaries. Instead, it is more important to describe the ion dynamics in and near the channel accurately, including the effects of image forces.