PNP Equations with Steric Effects: A Model of Ion Flow through Channels

Department of Applied Mathematics, Feng Chia University, 100 Wen-Hwa Road, Taichung, Taiwan 40724
Department of Mathematics, Taida Institute for Mathematical Sciences (TIMS), No. 1, Sec. 4, National Taiwan University, Roosevelt Road, Taipei 106, Taiwan
Department of Mathematics, Pennsylvania State University University Park, Pennsylvania 16802, United States
Department of Molecular Biophysics and Physiology, Rush University, 1653 West Congress Parkway, Chicago, Illinois 60612, United States
J. Phys. Chem. B, 2012, 116 (37), pp 11422–11441
DOI: 10.1021/jp305273n
Publication Date (Web): August 17, 2012
Copyright © 2012 American Chemical Society

Abstract

Abstract Image

The flow of current through an ionic channel is studied using the energetic variational approach of Liu applied to the primitive (implicit solvent) model of ionic solutions. This approach allows the derivation of self-consistent (Euler–Lagrange) equations to describe the flow of spheres through channels. The partial differential equations derived involve the global interactions of the spheres and are replaced here with a local approximation that we call steric PNP (Poisson–Nernst–Planck) (Lin, T. C.; Eisenberg, B. To be submitted for publication, 2012). Kong combining rules are used and a range of values of steric interaction parameters are studied. These parameters change the energetics of steric interaction but have no effect on diffusion coefficients in models and simulations. Calculations are made for the calcium (EEEE, EEEA) and sodium channels (DEKA) previously studied in Monte Carlo simulations with comparable results. The biological function is quite sensitive to the steric interaction parameters, and we speculate that a wide range of the function of channels and transporters, even enzymes, might depend on such terms. We point out that classical theories of channels, transporters, and enzymes depend on ideal representations of ionic solutions in which nothing interacts with nothing, even in the enormous concentrations found near and in these proteins or near electrodes in electrochemical cells for that matter. We suggest that a theory designed to handle interactions might be more appropriate. We show that one such theory is feasible and computable: steric PNP allows a direct comparison with experiments measuring flows as well as equilibrium properties. Steric PNP combines atomic and macroscales in a computable formulation that allows the calculation of the macroscopic effects of changes in atomic scale structures (size 10–10 meters) studied very extensively in channology and molecular biology.

Introduction


Ion channels are protein molecules that conduct ions (such as Na+, K+, Ca2+, and Cl that might be named bioions because of their universal importance in biology) through a narrow pore of fixed charge formed by the amino acids of the channel protein.(-2) Membranes are otherwise quite impermeable to natural substances, so channels are gatekeepers for cells. Channels are natural nanovalves that control a wide range of biological function.(3) Channels open and close stochastically, allowing ionic current to flow and forming a path for solute movement when they are open.(4-6) Only electrodiffusion moves ions through channels, so this biological system is like a hole in a wall that we should be able to understand physically.(7-9)
Ion channels are responsible for signaling in the nervous system. Ion channels are responsible for the coordination of muscle contraction and the transport of dissolved substances and water in all tissues. Each of these functions has been so important for so long that evolution has probably produced a nearly optimal adaptation within physical constraints and conserved it using the same design principle again and again.
Investigation of the physical mechanisms of current flow has just begun, although there is no shortage of descriptive metaphors in the literature of structural and molecular biology and biophysics.(2) The fundamental problem in a physical analysis is one of scale.(10) Mutations in single amino acids, which sometimes change only a handful of atoms involving perhaps just one permanent charge (radius of 0.1 nm), have dramatic biological effects. Such sensitivity comes as no surprise to the biologically oriented chemist or physicist.
Theories and simulations must account for the sensitivity of macroscopic function to atomic detail. Ion channels are nanovalves designed so that a few atoms, coded by the genetic blueprint of the protein, can control macroscopic function: that is what nanovalves (and channels) are all about. Theories and simulations must deal with 0.1 nm structural changes in charged groups that produce changes on the macroscopic scale of function. Structures as small as 0.1 nm move, and cannot be stopped from moving, via thermal (nearly Brownian) motion in 10–16 s. Their trajectories reverse direction infinitely often, whereas biology moves on time scales of 10–3 s or so in much simpler trajectories. The central physical issue is how to preserve this sensitivity to tiny structures while averaging over trajectories with such complex behavior over a 1013 range of time. Other scales also pose problems. Physics and biology and simulations of physics and biology must cope with a wide range of time scales and concentrations(10) as well as the immense range and strength of the electric field.(11-15)
The extremes of length, time, and concentration scales are all involved in the natural function of ion channels (or any nanovalve), so theory and simulations must deal with all of these extremes together. It is not likely that atomic-scale simulations by themselves will be able to deal with these, all together in finite time. Rather, reduced models of the type used widely in the physical sciences are more likely to be helpful in the foreseeable future.
A useful reduced model will include atomic-scale structural variables that determine macroscopic function. Sensitivity functions, determined by the theory of inverse problems, can help evaluate and construct reduced models. Biological function will be sensitive to important parameters and insensitive to others. The utility of these models can be evaluated by solving the relevant inverse problem for channels(16-18) using general methods.(19, 20)
So far, the most studied reduced model for ion flow in the bulk and in channels is the Poisson–Nernst–Planck (PNP) equation.(11, 12, 21-26) Although this model has had some success in dealing with experimental data,(22, 27-51) it does not include correlations introduced by the finite diameter of ions,(52) and these are of great importance in determining the selectivity of channels(9, 53, 54) and the properties of ionic solutions in general.(55-69) Crudely speaking, PNP is to nonequilibrium systems (such as channels) what Poisson–Boltzmann is to static systems: both are first approximations, useful in showing the crucial role of the electric field,(11, 12, 14, 70) the ionic atmosphere, and screening.(71) Neither are adequate models for ionic solutions such as seawater or the related solutions inside and outside biological cells.(52, 72-74)
Recently, work by Eisenberg, et al.,(75-78) built on the energetic variational theory of complex fluids,(77, 79-86) has developed a new set of PNP equations to implement and generalize an approach to selectivity started by Nonner and Eisenberg.(87-92) Nonner and Eisenberg considered a simplified model with ions (and side chains of the channel protein) represented as spheres of different finite sizes. They have shown in a long series of papers that important (static) selectivity properties of some significant types of ion channels can be explained with this model (reviewed in refs 9, 53, and 93). They have in fact constructed a single model with two adjustable parameters (diameter of channel, dielectric coefficient of protein, both set only once to unchanging values) using a single set of (crystal) radii of ions that fits the detailed and complex selectivity properties of two quite different types of channels, the heart CaV calcium channel(90, 94-98) and the nerve NaV sodium channel.(9, 91, 99) The theory accounts for the properties observed in solutions of different composition, with concentrations ranging from 10–7 to 0.5 M. When the side chains in the model are amino acid residues Asp-Glu-Lys-Ala, the channel has net charge of −1 although it is very salty (the magnitude of net charge is 3). The channel then is a DEKA sodium channel. When the side chains are Asp-Glu-Glu-Ala, then the channel is even saltier (the channel has net charge of −3). It is a DEEA calcium channel with quite different properties, although no parameters in the model are changed whatsoever, except the side chains that determine selectivity.
Most importantly, channels have been built according to the prescription of this theory, and they behave as predicted.(100-102) In studying another channel, the ryanodine receptor (of enormous biological importance as the final regulator of Ca2+ concentration in muscle and thus of contractions), Gillespie successfully predicted subtle and complex properties of selectivity and permeation before experiments were done and also the properties of drastic mutants in a wide range of solutions before the experiments were performed,(103-112) extending work on an earlier unsuccessful reduced model of the receptor that did not deal with the finite diameter of ions.(30, 32, 113, 114) One of the Meissner–Gillespie mutants reduces the permanent charge from 13 M to zero, yet Gillespie’s theory fits current voltage relations in several solutions with nearly the same approximately eight parameters as for wild type.
The calculations reported here extend the pioneering calculations of Hyon(75, 76, 78) by applying the energy variational approach to ion channels.(115) The bath and boundary condition treatments are somewhat different. A full 3D treatment is needed before the appropriate 1D approximation (particularly boundary conditions) can be determined without ambiguity.(25, 43, 88, 116-121)
Here, we simulate the properties of the family of calcium channels CaV (reviewed in refs 9 and 98) and sodium channels NaV(91, 99) using parameters already shown to fit a wide range of stationary (time-independent) experimental data in a variety of “symmetrical” solutions, which are designed so that current does not flow. Our results agree with previous equilibrium binding results and extend them to the world of current–voltage relations using a model and numerical methods that can be easily implemented on inexpensive computers. Current–voltage relations can be computed in a few hours on a notebook system.

Mathematical Model


Poisson–Nernst–Planck (PNP) Equations with Size Effects
The energy functional and the procedures for handling it with variational calculus are central to the energetic variational approach (EnVarA) formulated by Liu, more than anyone else. Liu’s approach is described in refs 75, 7982, 84, 115, and 122126. The “energy” of EnVarA is shown to correspond to the Helmholtz free energy of classical thermodynamics (in applicable equilibrium systems) in a recent article.(127) The application of EnVarA to membranes,(128) biological cells and tissues,(77) and ions, in channels and bulk solution, is described in detail in refs 1 and 75−78 in a way that is accessible to physicists and chemists without extensive experience with variational methods.
The energy functional for the ion channel is defined by(1)where ci and zi are concentration and valence for the ith ion (i = 1,···, N – 1); cN = cO–1/2 is the concentration for side chain O–1/2 (as in the glutamate side chain) with valence zN = zO–1/2 = −1/2 located in the filter only; is the electrostatic potential; kB is the Boltzmann constant; T is the absolute temperature; N is the number of ions; e is the unit charge; ρ0 is the permanent charge density; cO–1/2 is the concentration of the spherical “side chain” with valence zO–1/2 = −1/2 located in the filter only; V is the restraining potential that keeps the side chain inside the filter at all times; ai and aj are the radii of ions i and j; and εij is the energy coupling constant between ion i (including side chain O–1/2) and j (including side chain O–1/2). The last term is the repulsive part of the hard sphere potential that keeps ions apart.
The hard sphere repulsion characterizes the finite-size effect of ions and side chains inside the filter. These repulsive terms obviously depend on the chemical species and are called combining rules when they describe the interactions of different species. We discuss the combining rules later.
The basic reasoning is that the ion filter is so narrow that extra energy is needed to crowd ions into its tiny volume.(-9, 53, 54, 87, 89, 95-97, 99, 101, 102, 104, 129-131) Without finite size effects, the total energy will yield the traditional PNP equations.
The Euler–Lagrange equation (eq 1) will introduce PNP equations with size effects that depend on the global properties of the problem in the following way(2)(3)where flux Ji is(4)There is extra flux from the restraining potential that keeps side chains within the selectivity filter:(5)These equations are very similar to the drift-diffusion equations of semiconductors,(11, 12, 23, 24, 116, 121, 132-149) with the first term in the flux being the diffusion term and the second one being the drift term driven by the electrostatic potential of the field.
The third term involves a mutual repulsive force and interparticle hard sphere potential that is not typically found in semiconductor equations (although the semiconductor literature is so large that the volume exclusion of finite-sized holes and electrons is probably found somewhere that we do not know). This term includes forces usually called Lennard-Jones and depends globally on the properties of the solution everywhere because of the range of the integral on the right-hand side of eqs 4 and 5.
We call attention to the important role that the coefficients of these steric terms will have in determining biological function. The role of these steric terms will be somewhat different in our calculations from those in classical equilibrium analysis of ionic solutions using Monte Carlo simulations, for example. The cross terms in our expression appear as part of partial differential equations. These terms will then have effects on all terms in the solution of those partial differential equations. The integration process spreads out the effects of the cross terms. They propagate into everything as the partial differential equations are solved. The usual classical equilibrium treatment of Monte Carlo simulations is likely to produce different radial distribution functions from those produced by our differential equations, but a detailed comparison of MC and EnVarA calculations of absolutely identical models is necessary to determine the significance (or even existence) of this effect.
The steric cross terms, often called combining rules, turn out to have significant effects on the properties of ion channels. As mentioned, we use the Kong combining rules(150) to describe the repulsion between ionic spheres because they seem to be more accurate and justified(151, 152) than the more common Lorentz–Berthelot rules. It will turn out that the biological properties of ion channels are quite sensitive to these terms, but the choice of parameters seems to require detailed fitting to the properties of specific biological channels and transporters. We do not know what the effects of attractive terms (known to be present in bulk solutions) will be when we include them. We reiterate that these terms do not change the diffusion coefficients in our model.
One might think that simulations in full atomic detail (of molecular dynamics) would give good estimates of the combining rules, but sadly that is not the case. These simulations of molecular dynamics use combining rules (similar or identical to what we use in our reduced models) in the force fields of their own calculations.(150-152) We cannot use molecular dynamics to justify our combining rules because molecular dynamics itself uses those rules. It is possible that no one knows what cross terms to use in bulk solution. It seems likely that no one knows what cross terms should be used inside a channel or between side chains and ions. Indeed, it is difficult to conceive of experiments that might measure these inside channels with reasonable reliability.
Returning to the mathematical issues, we note that the singular convolution integral term can be regularized by a cutoff in the integration domain or simply by letting the integrand be 0 when |xy| ≤ ai + aj. However, the regularized term still produces numerical difficulty and is very time-consuming to compute, particularly at high dimensions, even if fast Fourier transform methods are used. We have tried. In addition, computing the convolution term generates an artificial boundary layer with a length of several grid spacings that needs to be filtered out and may have troublesome qualitative effects that are not so easy to remove by any local filtering because it has some of the properties of aliasing. Aliasing has devastating effects if not handled properly in both temporal and spatial systems that are treated as periodic when they are not.
We turn now to a simplified steric model that is much easier to compute because it uses only a local representation of interatomic forces. As we shall see, this steric model allows the computation of a large range of interesting phenomena, despite this simplification.
Local PNP Equations with Steric Effects (PNP–Steric Equations)
(6)and (7)where δ is a small number for the cutoff length, cδ is a dimensionless integration factor associated with δ, and d is the dimension. Here, the symmetry εij = εji has been assumed for notational convenience. To obtain this model, we have two important considerations: (1) the localization of the nonlocal size effects and (2) the finite truncations, which make the term local. Compared to the standard PNP equations, the PNP–steric equations have extra nonlinear differential terms (in the spatial variables) called steric effects. These represent the effective averaging/coarse graining of microscopic size effects for the macroscopic/continuum scales.
It should be clearly understood that coarsening terms of this sort are used throughout the chemistry literature, including within the simulations of molecular dynamics. The potentials of molecular dynamics simulations of proteins are not transferrable from quantum mechanical simulations of interatomic forces. The force fields that are used in every time step of an atomic-scale simulation include terms such as our εij justified only in the way that we have. Thus, molecular dynamics simulations depend on effective parameters, as do ours. Molecular dynamics simulations are no more derivable from quantum mechanics than are our models.
The main difficulty with eqs 25 comes from the convolution integral of the energy functional E with the following form(8)Usually, one would approximate the above integral by truncating the kernel 1/|xy|12 with the cutoff length δ, which causes kernel 1/|xy|12 to have a flat top when |xy| ≤ δ. To approach kernel 1/|xy|12, the length δ must be set as a small number tending to zero. One may expect that the smaller the cutoff length δ, the better the approximation. However, because of the effect of high-frequency Fourier modes, the approximation may lose the accuracy of numerical computations and makes numerical simulations difficult and inefficient.(153)
To deal with the effect of high-frequency Fourier modes, band-limited functions are used to cut off high-frequency Fourier modes. The functions act as optical filters by selectively transmitting light in a particular range of wavelengths. Band-limited functions play important roles in the design of signal transmission systems with many applications in engineering, physics, and statistics.(154) See also any textbook on digital signal or image processing. In ref 1, a class of band-limited functions depending on the length δ is found to approximate the kernel 1/|xy|12 and allow the derivation of the PNP–steric equations. The same approach can be used to modify the Poisson–Boltzmann equations used widely in physical chemistry, applied mathematics, and molecular biology.(70, 155) Once modified, these new steric Poisson–Boltzmann equations are particularly useful for the study of crowded boundary layers near charged walls, including the special behavior usually attributed to Stern layers. As the length δ goes to zero, the singular integral (eq 8) can be approximated by the integral Sδci(x) cj(x)dx with Sδ ≈ δ–12+d. Hence, the energy functional (eq 1) can also be approximated bywhich gives eqs 2, 3, 6, and 7, where gij = εij(ai+aj)12Sδ for i = 1,···, N–1 and gNj = εNj(ai + aj)12cδSδ.
The PNP–steric equations (eqs 2, 3, 6, and 7) are convection–diffusion equations following the energy dissipation lawwhere μi = ∑j = 1N gijcj is the chemical potential.
Note that eqs 6 and 7 contain no singular integrals such as in eqs 4 and 5 but have extra nonlinear differential terms. These extra nonlinear terms are crucial to simulating the selectivity of ion channels that cannot be found by simulating the standard PNP equations. Hence, the (local) PNP–steric equations are more useful than the standard PNP equations and are significantly more efficient and easy to work with than the global equations of the EnVarA treatment (eqs 25) discussed in the Introduction and Discussion sections of this article. These extra local terms in differential equations 6 and 7 have global effects when the differential equations are solved. Thus, pair correlation functions described by the solutions to differential equations 4 and 5 may have properties not present in classical equilibrium analyses containing local steric forces. (Classical analyses often deal only with positions and correlations and not with solutions of differential equations containing the forces.) Detailed fitting to specific experimental data is needed to compare the solutions of the local steric differential equations (eqs 6 and 7), the solutions of the more general differential equations (eqs 4 and 5), and the actual nonlocal phenomena of experiments.
We Adopt the PNP–Steric Equations
We replace eqs 4 and 5 with more approximate eqs 6 and 7 from now on. The real 3D geometry of an ion channel shown in Figure 1 is replaced with the simple axisymmetric geometry shown in Figure 2, with eqs 2, 3, and 6 valid in Ω and eq 7 valid only in Ωf, where Ωf is the filter part of the channel and Ω Ω. The associated boundary conditions are also shown in Figure 2, with the Dirichlet boundary conditions specified for both the ionic concentration and potential at the channel inlet (left end) and outlet (right end); no-flux boundary conditions are set for both the ionic concentration and potential at the side wall of the channel. Extra no-flux boundary conditions are set for side chains JO–1/2 = 0 at the interfaces (z = α and β) between the filter and the other part of the channel because side-chain molecules are free to move inside the filter only.
figure

Figure 1. Typical geometry configuration of an ion channel. The usual time scale for an ion passing through the channel is 200 ns. Specifically, a channel passing 1 pA of current with an occupancy of 1 ion has a mean passage time of 160 ns.

figure

Figure 2. Cartoon of the configuration of an ion channel with specified boundary conditions. Ω denotes the domain of the whole channel; Ωf denotes the filter part of the channel bounded by α ≤ z ≤ β and a side wall; is the unit outward vector normal to the side wall.

This model is meant to be nearly identical to that used in the many papers using Monte Carlo methods reviewed in ref 9, particularly the key papers.(54, 90, 93-96, 99, 156) The treatment of the region outside the channel and thus buildup phenomena are different from those in refs 75, 76, and 78.
Because no-flux boundary conditions are implemented for both the ionic concentration and potential at the side walls (orthogonal to the direction of current flow), this 2D problem (eqs 2, 3, 6, and 7) can be well approximated by a reduced 1D problem along the axial direction z, with a cross-sectional area factor A(z) included as in refs 25, 43, 88, and 116−121, described succinctly in ref 119, and perhaps included most carefully in the 3D spectral element calculations of Hollerbach.(-34, 43) Of course, some phenomena cannot be reproduced well in one dimension. See Figure 7 in ref 99.
The resulting 1D equations are(9)(10)(11)and(12)with boundary and interface conditions(13)and(14)The no-flux interface conditions for the side chains guarantee that side chains are not allowed to leave the filter. Mass conservation is preserved inside the filter:Also, in eqs 11 and 12, the Einstein relation is used for both the drift current and hard-sphere-potential flux, which is μe,i = Di/kBT, μLJ,i = Di/kBT(i = 1,···, N – 1), and μV,N = DN/kBT, where μe,i and μL J,i are the electrical mobility and the mobility associated with the hard sphere potential for ionic species i and μV,N is the mobility associated with the restraining potential for glutamate. Note that Di and i do not have to be homogeneous in space, nor do the dielectric coefficients of the solution and channel protein. We have not yet studied the effects of variation in these parameters, however. Usually, Di and i are set to 1/20th of their bulk solution values inside the channel filter but are set to bulk solution values in the rest of the channel as discussed at length in the Supporting Information and body of ref 109.
Dimensionless Equations
Nondimensionalization of the governing equations is especially important in discovering the structure, such as the boundary or internal layer, of the solution of PNP-type equations in advance, and so perturbation methods,(21, 24, 25, 135, 157) including some using the powerful and rigorous methods of geometrical perturbation theory,(24, 147, 148) have been used. We follow this work and scale the dimensional variables by physically meaningful quantitieswhere s denotes all length scales and L = rmin (the narrowest radius in the channel shown in Figure 2) unless otherwise specified. Note the scaling with respect to the physical dimension and not the Debye length. The Debye length varies with concentration and concentration varies with location and conditions. Equation 8 becomes(15)where Γ = λ2/L2 and the Debye length is λ = (εkBT/cmaxe2)1/2. Γ is the reciprocal of the length of the channel in units of Debye lengths. Note that Γ can vary dramatically with location and conditions if the contents of the channel vary with location or conditions. In the channels dealt with here, the contents of the channel are “buffered” by the charge of the side chains of the protein, most clearly in calcium channels EEEE and EEEA but also for the salty DEKA channel. Such buffering is not expected in all channels (e.g., potassium channels).
Equations 911 then become(16)(17)(18)We remove all of the tilde decorations () and rewrite the dimensionless governing equations (eqs 1518) for the mixture of Na+, Ca2+, Cl, and O–1/2 as follows:(19)(20)(21)(22)Note that ε̃ijcδδ̃–9(ãi + ãj)12 in eq 17 is lumped into ij along with cδδ̃–9 and is assumed to be the same for all species. Lennard-Jones parameters ε̃ij are obtained from the literature for alike species (i = j) and computed by Kong’s rule for unlike species (i j). They do not change diffusion coefficients in our model or calculations.
Channel Wall Shape Function
The wall shape function g(z) in Figure 2 can be arbitrarily specified, for example, g(z) = ((rmazrmin)/((a/2)2p))(za/2)2p + rmin, or nondimensionalized asWe remove all of the tilde accents () and write(23)In our calculation, we choose p = 4. The geometrical parameters are typically a = 50 Å, rmin = 3.5 Å, rmax = 40 Å, and D = DNa,bulk = 1.334 × 10–5 cm2/s . We set ε = 30ε0 inside the filter; in the rest of the channel, ε = εwater = 80ε0. For a typical cmax = 100 mM, we have Debye length λ = 8.48 Å and Γ=5.87 inside the filter and λ = 13.8 Å and Γ = 15.65 outside the filter. The above Γ values are based on L = rmin = 3.5Å.
The fact that Γ is not small implies that no internal or boundary layer is expected in the radial (transverse) direction. However, we sometimes choose L = a = 50 Å and Γ = 0.0288 inside the filter and Γ = 0.07668 outside the filter. Though Γ is not as small as in semiconductor devices, an internal/boundary layer is still expected in the lateral (axial) direction. We must not forget that the nonlocal hard sphere potential term (present in ionic solutions but not in semiconductors) may produce internal layers as well. See Figure 7 in ref 99.
A noticeable problem in all PNP (and Poisson–Boltzmann) theories without finite size are internal boundary layers near all boundaries with charge (permanent or induced). Such layers are customarily removed by introducing a single distance of closest approach for all ions, however ill-defined. Of course, no single distance of closest approach can deal with ions of different diameters, a problem that Debye, Hückel, and Bjerrum were evidently quite aware of. The need for multiple distances of closest approach (different for each ion and very nonideal, depending on each other and everything else in the system) means that the complex layering phenomena seen near walls of charge can have counterparts in channels. Indeed, complex layering is expected when ionic solutions are mixtures, such as the salt water in oceans or biology, and spatial inhomogeneities are present. Reference 158 is a gateway into the enormous chemical literature on layering phenomena near walls of charge. Reference 78 is a mathematical approach. Layering in ionic solutions might be able to produce nonlinear phenomena as important as pn junctions or even pnp junctions in semiconductors.
One might argue the validity of choosing rmin = 3.5 Å and wonder if an exclusion zone adjacent to the channel side wall for an ion sphere center with the thickness of an ion radius should be given extra consideration. In 3D models, this may be necessary and requires extra care in computation because the radial exclusion zone is different for ions of different sizes and charge: see Figure 7 in ref 99. However, in the 1D continuum model studied here, which ignores such ion-specific radial effects, our single radial zone would change only the value of rmin. Because all of the governing equations are scaled to be dimensionless, the effect of changing rmin would change only the value of Γ in eq 19. That in turn would only change the distribution of the electric potential. Γ is proportional to 1/rmin2. The permanent charge concentration ρ0 in eq 19 is also proportional to 1/rmin2. The net effect is simply reduced to the amplification or shrinking of the influence of the contribution of the electric potential exerted by ion distributions. The permanent charge concentration ρ0 is generally much larger than all ion concentrations and dominates the distribution of electric potential, we imagine. We then reason that a minor change in rmin would not affect our results significantly.
Numerical Methods
Now we apply the multiblock Chebyshev pseudospectral method(153) together with the method of lines (MOL) to solve eqs 1922 with the associated boundary/interface conditions (eqs 13 and 14). These governing equations are semidiscretized in space together with boundary/interface conditions.
The resulting PNP delta representation is a set of coupled ordinary differential algebraic equations (ODAE’s). The algebraic equations come from the boundary/interface conditions that are time-independent. The resulting ODAE’s have an index of 1, which can be solved by many well-developed ODAE solvers. For example, ode15s in MATLAB is a variable-order, variable-step index-1 ODAE solver that can adjust the time step to meet the specified error tolerance and integrate with time efficiently. The numerical stability in time is automatically assured at the same time. The spatial discretization here is performed by the highly accurate Chebyshev pseudospectral method with the Chebyshev Gauss–Lobatto grid and its associated collocation derivative matrix. To cope with the computational domain of side chains being strictly within the [α, β] region and the conformation of grids, we need to use domain decomposition. We decompose the whole domain into [0, α], [α, β], and [β, a].
The extra interface conditions from this domain decomposition for ions are implemented simply by continuity of ion concentration and the associated flux. Finally, the Poisson equation for the electric potential is solved by the direct inverse at every time step, which is easy because it is only 1D.

Results


Here we consider a calcium channel (EEEE) with four glutamate side chains and eight O–1/2 particles that are free to move inside the filter, which is essentially the model introduced by Nonner and Eisenberg(53, 87-89) and used by them and their collaborators since then (reviewed in refs 9 and 72). Our goal is to demonstrate the feasibility of a PNP–steric model and the range of phenomena that can be calculated despite its local approximation. Note that the effects of a local approximation on the right-hand side of partial differential equations are not the same as the effects of a local approximation in a classical analysis of the BBGKY hierarchy. The much-needed detailed comparison with experimental results lies in the future. We are particularly interested in the effects of the steric parameters that we call εij in eq 1 and then the effects of gij as well, so we concentrate on steady-state results. Transients of the type previously reported(75) have been computed and will be reported separately.
The channel geometry used for the current 1D simulations is shown in Figure 3, and the parameters used are shown below. These parameters are not changed in the calculations (e.g., the diffusion coefficients are always the same and are not changed as the interaction (Kong) parameters are changed).
figure

Figure 3. Channel geometry. A precise specification of the geometry of our model.

Parameters


Filter radius, 3.5 Å; filter length, 10 Å; diffusion coefficients (cm2/s) in the nonfilter region, DNa+ = 1.334 × 10–5, DCa+2 = 0.792 × 10–5, DCl = 2.032 × 10–5, and DO–1/2 = 0.76 × 10–5; diffusion coefficients in the filter region, 1/20 of the above values (see Gillespie,(109) particularly the Supporting Information); ion radii, aNa+ = 0.95 Å, aCa+2 = 0.99 Å, aCl = 1.81 Å, and aO–1/2 = 1.4 Å; relative dielectric constant, 30 in the filter region and 80 in the nonfilter region.
Dimensionless restraining potential for O–1/2 inside the filter required by eq 22(24)where γ is a scaling constant that makes V reach Vmax at z = α and β.

Boundary Conditions


Voltage in both reservoirs, 100 mV; Na+ concentration in both reservoirs, 100 mM; Ca2+ concentration in both reservoirs, 0.001–10 mM. The Cl concentration in both reservoirs is chosen to keep the whole solution electrically neutral.
To measure the selectivity, as reported in the biological literature of calcium channels, we need to compute the Ca binding ratio(25)Turning to the precise description of the spherical ions, we need to specify many Lennard-Jones-style parameters. There are 10 parameters of the gij’s that must be chosen, without specific experimental data relevant to the interior of a channel. (We note that this problem is not ours alone. The same situation is faced for any atomic-scale model. No one knows how to choose force fields of molecular dynamics that are suitable for the special conditions inside an ion channel. If one follows the convention of molecular dynamics and uses force fields that depend only on the distance between two atoms, this problem is particularly serious. Note that dielectric boundary forces are almost always of great importance in confined systems such as ion channels.(159, 160) It is not likely that dielectric boundary forces acting on two ions can be well approximated as a function of only the distance between two ions.)
We use the energy well εij data from the traditional Lennard-Jones (12–6 rule) potential as a reference. From the work in refs 151 and 152. we choose εNa,NaCl,ClCa,CaO,O = 1:1:1:1.56. Kong’s combining rule(150) seems to be the best for ionic solutions,(151, 152) with σij = (a+ aj). This gives us the εij values for the rest of the cross hard-sphere potential terms.and similarlyWe can see that gCl,j (especially gCl,Cl) is much larger than the other gij terms because the size is increased so dramatically by the exponent in (ai + aj)12. These large values would make the governing equations very stiff in numerical properties and hard to integrate in time. To resolve this numerical difficulty, we remove all of the hard sphere forces involving Cl, which means that gCl,j = 0, j. This approach can be rigorously justified because Cl is usually very dilute inside the filter, as are all co-ions in ion exchangers,(161) because of the electrostatic repulsion from the highly concentrated permanent charge of eight O–1/2.
Choosing Self-Coupling Coefficients gNa,Na
The above coupling terms gij are derived as ratios, and we still need to determine actual gij values by choosing self-quantities gNa,Na. The self-coupling gNa,Na is known only from its relationship to its effective ion radius: it is proportional to (aNa + aNa)12. Larger gNa,Na values imply a stronger hard sphere potential and more pushing among particles. Smaller gNa,Na implies less interaction and pushing among particles. If we choose gNa,Na to be too small, then the finite-size effect will be trivial and (judging from previous work cited above) the correct selectivity of calcium ions will not be found. If we choose gNa,Na to be too large, then repulsion will be too strong and the profile of concentration of all species (inside the filter) will be flat. Selectivity, as biology knows it, will not be present.
A numerical experiment (Figure 4 and Table 1) shows how gNa,Na changes the Ca binding ratio. The conditions of each case are stated in the figure captions. Note that many of the cases considered below correspond to different physiological states that may have profound implications on function. Cycling between such states has been the explanation of most behavior of channels and transporters for some 60 years, since Hodgkin and Huxley (who, one notes, did not use such explanations themselves). But the states in those explanations are ad hoc, arising as inputs of models from wisdom and experimentation on macroscopic systems, not from a direct physical knowledge of channels. The states shown here in our calculations arise without human intervention or wisdom. Rather, they arise as outputs of direct self-consistent calculations. Sometimes it is better to be wise, and sometimes it is not. The choice between handcrafted traditional models of states and direct calculations of ions that are sometimes in definite states should be made, in our view, by success or failure in explaining and predicting experimental results with models. The models should be parsimonious and specific, of course, so they can be falsified, at least in principle. Otherwise, they are more poetry than science.
figure

Figure 4. Species concentration distributions with various gNa,Na values. For Vmax = 200: (a) gNa,Na = 0 and Ca2+ binding ratio = 0.60214; (b) gNa,Na = 10–4 and Ca2+ binding ratio = 0.59418; (c) gNa,Na = 10–3 and Ca2+ binding ratio = 0.75433; (d) gNa,Na = 10–2 and Ca2+ binding ratio = 0.86109; (e) gNa,Na = 10–1 and Ca2+ binding ratio = 0.82580; (f) gNa,Na = 1, Ca2+ binding ratio = 0.71644, and the symmetrical symmetric boundary conditions are [Ca2+]L = [Ca2+]R = 1 mM, [Na+]L = [Na+]R = 100 mM, L = R = 100 mV. Note the 4-fold scaling of the [O–1/2] concentration.

Table 1. Effect of Increasing εglobal on the Ca Binding Ratio with [Ca2+]L = [Ca2+]R = 1 mM [Na+]L = [Na+]R = 100 mM, L = R = 100 mM, and Vmax = 200
gNa,Na010–410–310–210–11
gNa,Cl000000
gCl,Cl000000
gNa,Ca06.41 × 10–56.41 × 10–46.41 × 10–36.42 × 10–26.42 × 10–1
gCl,Ca000000
gCa,Ca01.64 × 10–41.64 × 10–31.64 × 10–21.64 × 10–11.64
gNa,O–1/208.19 × 10–48.19 × 10–38.19 × 10–28.20 × 10–18.20
gCl,O–1/2000000
gCa,O–1/201.00 × 10–31.00 × 10–21.00 × 10–11.00341.00 × 101
gO–1/2,O–1/201.63 × 10–21.64 × 10–11.641.65 × 101.64 × 102
Ca binding ratio0.6020.5940.7540.8610.8250.72
We can summarize the observations in the following text.
(1) From eqs 21 and 22, the flux of ion species consists of diffusion, migration (i.e., electrostatic drift driven by electric potential), and particle-to-particle steric-effect interaction. The typical steric-effect flux of −iiij(∂j/∂) can be seen as a chemical potential drift exerted on ion species i by species j, in which i = j is also allowed. The steric-effect flux generally includes a flux coupling between species i and j, and this coupling is not captured in plain PNP and DFT–PNP theory, where fluxes of one species are driven only by gradients of the chemical potential of that one species. Here, everything interacts with everything else: the flux of one species is driven by gradients of the chemical potential of another species, even though we use (nearly) the same constitutive (NP) relation for transport as in classical PNP or DFT–PNP. Our chemical-potential drift term is not like electrostatic drift. The electrostatic drift can flow uphill or downhill along the electric potential, depending on the sign of the valence zi being negative or positive. The chemical potential drift, however, always flows downhill along the chemical potential unless iij is negative, which is impossible if the particle-to-particle steric-effect interaction is always repulsive (not attractive). If particles push each other away, then peaks of ion concentration (for different species) tend to separate from each other as best they can, unless frustrated or overcome by additional electrostatic force. Here, in the present case, Na+ and Ca2+ chiefly feel a strong push from O–1/2 as gNa,Na gets larger (because aO–1/2 is large and O–1/2 is kept inside the filter) as well as the electrostatic attractive force from O–1/2. This can be clearly seen in Figure 4. Note the 4-fold scaling of O–1/2 concentration. In Figure 4a, for gNa,Na = 0, the concentration profiles of O–1/2, Na+, and Ca2+ reach equilibrium (at the minimum of total energy) simply by diffusion and the electrostatic force because the extra restraining force is felt by glutamate O–1/2 only. O–1/2, Na+, and Ca2+ all form single-peak concentration profiles in the same region of the channel for the same range of z. Physically, Na+ and Ca2+ are attracted to focused O–1/2 by the electrostatic force. The attraction for Ca2+ has a larger effect than the attraction for Na+ because Ca2+ is divalent. Also, Na+ and Ca2+ at the same time repel each other by electrostatic (and steric) forces.
This complex balance of forces can produce a wide range of behavior that varies a great deal as concentrations and conditions are changed. The biological function of channels and transporters has been defined experimentally for many decades by their behavior in complex ionic mixtures of variable composition as different voltages are applied across the cell membrane. We have not yet investigated the properties of our model as the concentrations of ions are made unequal on either side of the channel, as the electrical potential is varied, or, most importantly, as different species of ions are included in ionic mixtures on both sides of the channel.
The interaction forces (without the attractive component) may be responsible for many of the single-file properties of channels. More complete descriptions of Lennard-Jones forces include an attractive component. The interaction forces with the attractive component might (conceivably) produce the phenomena that define transporters, whether they are co-transporters or counter-transporters.
(2) As gNa,Na becomes larger in Figure 4b–f, the primary force is still the electrostatic attraction of Na+ and Ca2+ by O–1/2 but is modified by the finite-size effect (hard-sphere force). This electrostatic force will make peaks of Na+ and Ca2+ occur in the selectivity filter (i.e., in the same region of the channel that contains the O–1/2 side chains). For Na+, Ca2+, and O–1/2, the hard-sphere forces between ions of the same species will make the concentration profiles flatter as gNa,Na gets larger, when particles feel the push from similar particles. Flattening is clearly seen in Figure 4b–f.
(3) As gNa,Na gets larger in Figure 4b–f, Ca2+ pushes Na+ away from the middle part of the filter and forms a depletion zone and double-peak profile for Na+ in Figure 4c–e. Depletion zones of this sort have profound effects on the selectivity of ion channels in Monte Carlo simulations. See Figure 6 in refs 99 and 54. Depletion zones have profound effects on the behavior of transistors(132-135, 162) and the selectivity of channels.(9, 54, 99) A single transistor can have qualitatively distinct properties (e.g., gain, switch, logarithm, and exponential) for different boundary electrical potentials (bias for one transistor; power supply more generally) because the different boundary potentials produce different arrangements (layering) of depletion zones. Each arrangement of layers or depletion zones makes the same transistor a different device with a different device equation, corresponding to a different reduced model for the transistor. Different reduced models are appropriate for different conditions and have different functions. The function of the depletion zones found here is not yet known nor is the pattern or effect of cycling through structures known, but the complex properties of the Ca/Na exchanger, wonderfully characterized by Hilgemann,(163-166) immediately come to mind.
The major mechanism in Figure 4c–e is still that Ca2+ is more attracted to O–1/2 than to Na+. However, with extra help from the interspecies hard sphere force between Na+ and Ca2+ (in addition to the electrical repulsive force between Na+ and Ca2+), Ca2+ is able to push Na+ out of the filter. Note that again Figure 4f is an exception because all concentration profiles in it are very flat. In Figure 4f, the interspecies repulsive forces are greatly reduced and Na+ resides in the middle of the filter along with Ca2+. Note that Ca2+ forms a single-peak concentration profile in all of these cases without splitting into two peaks. Na+, however, can form a double peak when gNa,Na increases. The splitting occurs when the electrostatic attraction by O–1/2 is large enough to survive the push of different species from O–1/2 and Na+ in addition to the repulsive electrostatic force from Na+.
The singular behavior in Figure 4f may have direct functional significance. The splitting of a single peak into a double peak can create a depletion zone that can dominate the channel behavior (even though it is very small) because it is in series with the rest of the channel. A depletion zone can block flow and create switching behavior, as it does in transistors.
The depletion zones could help create the many states of a channel, identified as activated, inactivated, slowly inactivated, blocked, and so forth.(2) The depletion zones could be responsible for many of the similar (but correlated) states identified in classical experiments on transporters. In our calculations, of course, states arise as outputs of a self-consistent calculation and as a result of theory and computation, not as handcrafted metaphors summarizing the experimental experience and wisdom of structural biology and classical channology.(2)
Despite our enthusiasm and focus on our work, it must be clearly understood that our treatment of correlations is incomplete. We leave out many of the correlation effects of more complete variational treatments(1, 75, 76) and the more subtle correlations in the BBGKY hierarchy, derived for nonequilibrium systems(67) such as PNP (without finite-sized ions) from a Langevin description of trajectories in refs 278, 167, and 168. Our treatment is mathematically fully self-consistent but physically incomplete in its treatment of correlations and of course chemical interactions as well. The classical discussion of the BBGKY hierarchy and its treatment of correlations is not directly applicable to our analysis, however. The correlation forces of the hierarchy appear as driving forces in our partial differential equations, so the results of those forces are spread through all the terms of the solution of our partial differential equations. Thus, the effects of the correlations are likely to be more widespread than the effects of classical equilibrium analysis. Only detailed fitting of theory to data will show which correlations must be included in our model to explain the experimental data. Conclusions from equilibrium analysis may not apply.
(4) From Table 1, the Ca2+ binding ratio starts from 0.60214 when the finite size effect is zero (i.e., gNa,Na = 0). Affinity for Ca2+ in the filter region shows itself, even without the finite-size effect. This is obviously because Ca2+ feels a stronger electrostatic force than Na+ because of the larger valence. The valence effect dominates even though the concentration of Ca2+ is much lower than that of Na+ in reservoirs. The fact that valence effects overwhelm concentration effects when studying divalents has been known for at least 100 years. The binding ratio of Ca2+ decreases, increases, and then decreases again as gNa,Na increases, which shows the influence of the finite-size effect. The Ca2+ binding ratio roughly reaches its maximum at gNa,Na = 0.01 with a value of 0.861. This is far larger than the 0.602 that occurs when the finite size effect is zero and helps generate the affinity of Ca (selectivity). These effects can be further seen from the computational results shown later, when gNa,Na = 0.01.
We have not yet studied the effects of concentration gradients. Note that in biological systems, gradients of Ca2+ are large and have profound effects in experiments. Calcium concentrations outside cells are typically 2 × 10–3 M, whereas those inside cells (cytoplasm) are <10–7 M. There are many compartments within cells (vesicles, mitochondria, endoplasmic reticulum, and sarcoplasmic reticulum) essential to living function that maintain distinctive concentrations of Ca2+ without which they cannot function. We anticipate complex behavior of concentration profiles of ions within the selectivity filter when our model is studied in realistic ionic mixtures. It seems unlikely that these are uninvolved in biological function, however obscure that involvement seems today and however difficult it is to discover. It seems wise to do calculations under conditions in which the systems can actually perform their natural function, and it is unwise to simulate conditions in which biological systems are known not to function properly.
Ca Binding Curve
In these calculations, we first studied how finite size effects change the Ca binding curve, and the results are shown in Table 2 and its associated Figure 5. We choose gNa,Na = 0 and an appropriate finite size effect by choosing gNa,Na = 0.01 to calculate the binding curves of Na+ and Ca2+ respectively. Vmax = 200 mV as above, boundary conditions L = R = 100 mV, and concentration of Na+ inside and outside = 100 mM.
figure

Figure 5. Binding curves corresponding to Table 2.

Table 2. Ca Binding Ratio vs [Ca2+]L = [Ca2+]R with gNa,Na = 0 (No Finite-Size Effect) and gNa,Na= 0.01 (Having a Finite-Size Effect)a
[Ca2+], mMCa binding ratio gNa,Na = 0Ca binding ratio gNa,Na = 0.01
10–69.2286 × 10–64.4525 × 10–3
10–59.2257 × 10–53.1819 × 10–2
10–49.1970 × 10–40.12233
10–38.9241 × 10–30.28671
10–27.0641 × 10–20.49926
10–10.291710.70778
10.602140.86109
100.828160.94502
1000.936610.98080
a

Vmax = 200, L = R = 100 mV, and [Na+]L = [Na+]R = 100 mM.

The results are shown in Table 2 and its associated Figure 5. The finite-size effect greatly enhances the selectivity, as observed previously in Monte Carlo simulations. The concentration profiles for different species are shown for both cases in Figures 6 and 7, respectively. Note the 2-fold scaling of the O–1/2 concentration in both figures. Figure 6 shows the increase in Ca2+ and the decrease in Na+ as the Ca2+ concentration increases (in the baths on both sides of the channel). Increases are due to the interplay of diffusion and electrostatic forces. There are no special chemical forces in our calculations. Binding forces are the output of our calculation, not the input as in so many treatments of selectivity.
figure

Figure 6. Species concentration distributions under various [Ca2+]L = [Ca2+]R with gNa,Na = 0 (no finite-size effect). Vmax = 200, L = R = 100 mV, and [Na+]L = [Na+]R = 100mM: (a) [Ca2+]L = [Ca2+]R = 10–7 M; (b) [Ca2+]L = [Ca2+]R = 10–6 M; (c) [Ca2+]L = [Ca2+]R = 10–5M; (d) [Ca2+]L = [Ca2+]R = 10–4 M; (e) [Ca2+]L = [Ca2+]R = 1 mM; (f) [Ca2+]L = [Ca2+]R = 10mM. Note the 2-fold scaling of the O–1/2 concentration.

figure

Figure 7. Species concentration distributions under various [Ca2+]L = [Ca2+]R with gNa,Na = 0.01 (with finite-size effect). Vmax 200, L = R = 100 mV, and [Na+]L = [Na+]R = 100 mM: (a) [Ca2+]L = [Ca2+]R = 10–7M, (b) [Ca2+]L = [Ca2+]R = 10–6 M, (c) [Ca2+]L = [Ca2+]R = 10–5 M, (d) [Ca2+]L = [Ca2+]R = 10–4 M, (e) [Ca2+]L = [Ca2+]R = 1 mM, (f) [Ca2+]L = [Ca2+]R = 10 mM. Note the 2-fold scaling of the O–1/2 concentration.

In our model, Na+ and Ca2+ are both attracted to confined O–1/2. They repel each other at the same time. Na+ and Ca2+ both form only single-peak concentration profiles; therefore, the depletion of Na+ in the middle of the filter occurs only at the largest peak as the Ca2+ concentration increases (on both sides).
Figure 7 shows the extra influence of the finite-size effect compared to Figure 6. The single-peak profile of Ca2+ found in all cases means the pull from O–1/2 by electrostatics survives the electrostatic repulsive force from Na+ as well as the hard-sphere pushes from both O–1/2 and Na+. The pull is so strong that there is no splitting and no double-peak profile.
However, the electrostatic pull for Na+ from O–1/2 is much smaller than its counterpart for Ca2+. The combination of the hard-sphere pushes from both O–1/2 and Ca2+ and also the electrostatic repulsive force from Ca2+ have qualitative effects. The concentration profile of Na+ changes from a single-peak to a double-peak profile as the Ca2+ concentration increases (on both sides of the channel). Also, a depletion zone of Na+ inside the filter is observed as the Ca2+ concentration increases (on both sides of the channel). A depletion zone, which arises when a peak splits in two, can have profound functional consequences because it is in series with the entire channel. A series barrier can entirely block the current flow.
DEEA Ca2+ Binding Curve
We have also computed the Ca2+ binding curve of a mutant sodium channel (DEEA) with three glutamate side chains. These are represented as six O–1/2 particles freely moving inside the filter as in previous work, consisting mostly of Monte Carlo simulations previously cited. Figure 8 shows the effect of a −4e side chain in EEEE and a −3e side chain in DEEA on the Ca2+ binding curve. Obviously, EEEE with a −4e side chain has a slightly larger affinity for Ca2+ than does DEEA with a −3e side chain. This DEEA binding curve, employing the PNP–steric model, agrees well with its counterpart in ref 75 using the PNP–LJ model and its counterpart in ref 99 using Monte Carlo simulations.
figure

Figure 8. Binding curves of EEEE (−4e) and DEEA (−3e).

DEKA Ca2+ Binding Curve
Here we compute the Ca2+ binding curve of the sodium channel (DEKA) with two glutamate side chains (four O–1/2 particles) and one lysine side chain (one NH4+ particle) free to move inside the filter. The Ca2+ binding curve is shown in Figure 9, and the associated species concentration profiles are shown in Figure 10. Note that the scaling of [O–1/2] is the same as the scaling of other concentrations in Figure 10, unlike those in Figures 4, 6, and 7. The loss of affinity for Ca2+ is obviously due to the existence of a lysine side chain with a +1e charge, though the net charge on glutamate and lysine side chains taken together is still −1e. We also computed an artificial EAAA channel with only one glutamate side chain (two O–1/2 particles) in which the net permanent charge on the filter is −1e, to correspond precisely to the experimental situation as discussed in ref 99 and references cited therein. EAAA still has a much higher affinity for Ca2+ than for Na+ (data not shown). The DEKA binding curve, employing the PNP–steric model, agrees well with its counterpart in ref 75 computed using the PNP–LJ model and also with its counterpart in ref 99 computed using Monte Carlo simulation. The concentration profiles of individual species shown in Figure 10 also resemble those in ref 99. We have not yet performed the calculations with multiple ion species needed to evaluate the selectivity of the DEKA model for K+ ions.
figure

Figure 9. Binding curves of DEKA (−1e).

figure

Figure 10. Species concentration distributions under various [Ca+2]L = [Ca+2]R with gNa,Na = 0.01 (having a finite-size effect). Vmax = −200 for glutamate side chain, Vmax = 200 for lysine side chain, L = R = 100 mV, and [Na+]L = [Na+]R = 100 mN: (a) [Ca2+]L = [Ca2+]R = 10–7M, (b) [Ca2+]L = [Ca2+]R = 10–6M, (c) [Ca2+]L = [Ca2+]R = 10–5M, (d) [Ca2+]L = [Ca2+]R = 10–4M, (e) [Ca2+]L = [Ca2+]R = 1 mM, and (f) [Ca2+]L = [Ca2+]R = 10 mM. Note that the scaling of [O–1/2] is the same as the scaling of other concentrations in this figure, unlike that in Figures 4, 6, and 7.

Discussion


Ion channel function depends on the properties of ionic solutions and ions in channels, along with the properties of the channel protein itself, so it is necessary to relate our work to previous work in each field, emphasizing the properties of ions in bulk solutions, ions in proteins, and proteins that determine the biological function in our model of channels, if not in the real world.
Relation to Classical Work on Ionic Solutions, Poisson–Boltzmann, and PNP
The limitations of Poisson–Boltzmann and PNP models of ionic solutions have been known a very long time to the physical chemistry community but seem not to be so well known to either applied mathematicians or biophysicists. Exhaustive references to the literature are in refs 52, 59, 65, and 169−171. Applied mathematicians understandably are attracted to the simplicity of the Poisson–Boltzmann–PNP equations and view them as a starting point for more realistic treatments.(70) Biophysicists(2) use the independence principle that worked so well(-172, 173) when applied to membranes in which ions flow through separated, independent protein channels.(174-177) When channels are not selective(178, 179) or two types of ions flow through one channel, as in classical ligand-gated acetylcholine channels nAChRs,(180) the independence does not apply.
The independence principle is a restatement of Kohlraush’s law of a century ago, which does not apply to bulk ionic solutions of the type found in biology (eq 3.27a,b on p 125 of ref 67). References to the physical chemistry literature include refs 21, 57−63, 66−69, 74, 129, 158, 169, 171, and 181−230. Many of these papers are strictly experimental, presenting compilations of physical chemistry data. These papers also show that ionic solutions are not well described by Poisson–Boltzmann or PNP, if they contain divalents, multivalents, or mixtures of monovalents. (Biological solutions are mixtures usually containing divalents.) Interactions are strong in all ionic solutions because they all satisfy global electroneutrality. Thus, ionic solutions are nothing like ideal, so the law of mass action, for example, does not apply as usually used with rate constants independent of concentration. Note that rate constants in complex models will depend on each other as well as on concentration because the electric field in one part of a channel (described by one rate constant) will interact with charges in another part of the channel (described by another rate constant). It is that variation of the electric field that allows Kirchoff’s current law (and its generalization to the Maxwell equations(-232-234)) to be true. The electric field is long-range and cannot be broken into independent spatial components as it is in most classical treatments.(2)
Classical Work on Channel Proteins: Permeation
Currents permeate biological membranes by flowing through channel proteins that are either open or closed. A single ionic channel controls the current by opening and closing (spontaneous gating),(5, 6, 235-237) thereby making a random telegraph signal(238) that has been studied in enormous detail for many (hundreds or thousands of) channel types(178, 179, 239, 240) using the wonderful techniques of single-channel recording, patch clamp,(6, 241-245) and bilayer reconstitution.(246)
Sadly, the structures and mechanisms that produce this gating are still mostly unknown. However, progress is at hand.(247-251) Special structures modulate spontaneous stochastic gating in most channels to produce the macroscopic gating properties of classical electrophysiology.(2, 4, 5, 252, 253) The properties of macroscopically modulated gating are complex, as is clear from the variety and number of complex schemes in the classic text(2) (chapters 18 and 19 and Figure 19.11). Some of these schemes involve nearly 100 ill-determined rate constants (Figures 18-11 and 18-12 in ref 2) and have attracted the attention of literally hundreds of investigators over many decades.(239, 240) So far, no theoretical model can explain gating and selectivity using the fundamental physics that is described (crudely) in the PNP equations, but this situation may change.
Despite our ignorance of the mechanism of gating, the phenomena of spontaneous gating is remarkably clear, one might even say crystal clear, despite the amorphous structures involved. Once the single channel is open, the current through the single channel is remarkably stable. Single-channel (mean) currents are independent of time from say 10 μs to 10 s or longer, strongly suggesting that the channel protein has only one structure over a range of (at least) 106. A glance at an MD simulation of channels or the numerical values of energies computed from Coulomb’s law or Lennard-Jones potentials suggests that the structure of the channel pores and the pore walls must be very constant indeed over these time scales. A change in radius of 3% would produce a change in current of at least (and probably much more than) 9%. Single-channel currents are routinely resolved to within 2% (and can be resolved much better(5, 242, 243, 254-256) if necessary because the stability is nearly perfect and signal-to-noise ratios are often larger than 50). The complexity and ignorance of the gating mechanism reappear when we consider the time course of the opening and closing processes themselves (e.g., on a much faster time scale(244) or in cooled systems(245)). The opening and closing processes do not in fact have well-defined time courses, and nothing seems to be known about their physical origin in either case.
Classical Work on Channel Proteins: Permeation and Selectivity
Once open, channels select between ions of different chemical types. Channels allow only some types of ions to flow, even though the different chemical types are quite similar. For example, Na+ and K+ ions differ only in diameter; Na+ and Ca2+ ions differ only in charge. Simple models do surprisingly well in dealing with the selectivity of some types of channels. This work was reviewed in the Introduction of this article and elsewhere.(7, 9, 72)
Dealing with Biological Reality
It is important that the study of permeation and selectivity be carried out in the context of specific channels using parameters known to fit a wide range of experimental data properly.(9, 16-18, 90, 94-97, 99-112)
It is a surprise, particularly to structural biologists and traditional channologists, who customarily deal with metaphors(257-259) and not quantitative fits to data that powerful results, with quantitative fits to data in important cases, can arise from the Nonner and Eisenberg models with their very simple structures. We do not know why these models work, but one reason may be that the structures are the computed consequences of the forces in the model so the structures of the channel protein and of the ionic solutions are always exactly self-consistent. Even tiny deviations in the location of side chains from their free-energy minima produce large energy and functional effects.(54, 93) Exact self-consistency between channel protein and ionic solution seems to be necessary to make reasonable models. We suspect that exact self-consistency is why some simulations fit some data so well.
If exact self-consistency is necessary to make reasonable models, then classical models in much of molecular biology(260-264) will need to be reconsidered, even in much of chemistry,(73, 74) because classical chemical and biochemical models are almost never self-consistent. They almost never calculate the electric field from the charges present, let alone deal self-consistently with boundary conditions, steric forces, or the resulting interactions of everything with everything else. Our approach in this article represents the ionic atmosphere around an ion consistently in a simplified way using the approximated LJ potential instead of the original LJ potential. Surely, we have not included all of the correlations among ions. Only detailed fitting to large amounts of data will show whether we have captured enough correlations, and captured them well enough.
It is also possible that the Nonner and Eisenberg models work well because the community of scientists working on them has recapitulated evolutionary history. Perhaps those scientists have stumbled on the adaptation that biological organisms found eons ago, as evolution selected mutations that allowed cells to live and reproduce. It is even possible (for the same reason) that simple, nearly 1D models will capture most of what we need to understand time-dependent nonequilibrium properties of channels.
Nonequilibrium Treatments
An important advantage of the methods considered here is their extension to nonequilibrium in a mathematically precise and defined way, always fully self-consistent. Other approaches depend on physical approximations that are not self-consistent. They were the best that could be done at the time but cannot substitute for self-consistent treatments, in our view. For example, the DFT–PNP method(18, 105, 107, 108, 110-112) is not self-consistent (i.e., it does not precisely satisfy sum rules of statistical mechanics(265, 266) or Poisson’s equation and boundary conditions) and apparently leaves out relaxation and dielectrophoresis terms of the (more or less) self-consistent Debye–Hückel–Onsager equation. See the classical work(267-278) and the textbook(207) (p 282, Figure 7.7). DFT–PNP indeed assumes local equilibrium, as do other approaches using combinations of simulation and PNP equations,(279-281) although it computes global flux.
It must be clearly understood that any assumption of local equilibrium is also an assumption of local zero flux. It is not clear how a system can have zero local flux and long-range substantial flux as does DFT–PNP, particularly when the system is a nanovalve connected in series to a high-impedance entry process and macroscopic baths. DFT–PNP is inconsistent in its treatment of both flux and electrostatics. Adopting models that are inherently inconsistent is dangerous because the results of calculations can depend on how the inconsistency is resolved, and that resolution may be presented sotto voce, or chosen without conscious thought.
It is striking that very successful PNP calculations in a closely related field such as computational electronics do not use inconsistent assumptions (such as local equilibrium and global nonequilibrium) and always satisfy Poisson and boundary conditions with great accuracy. Simulations in computational electronics directly solve the relevant equations and so are fully self-consistent. Otherwise, they have difficulty accounting for the function of semiconductor devices that arises from small differences in large forces. Calculations of computational electronics account for the macroscopic function using atomic-scale models.(8, 9, 11, 12, 134, 136, 149, 282-284) Most treatments of ions in water and channels have been much less successful, perhaps because they are inconsistent.
Resolving inconsistencies can be a difficult task. It took a detailed stochastic analysis (lasting many years) of a second-order Langevin equation with doubly conditioned nondifferentiable Brownian trajectories(285, 286) to resolve a similar inconsistency in an analysis of noninteracting particles. And the results of that analysis were not at all what had been anticipated, although they were pleasingly simple when interpreted correctly with the classical theory of mass action.(74) The stochastic analysis allowed one to derive the law of mass action, but with variable rate constants that were specific functions of everything in the system, as given by the analysis. It is not clear how one can evaluate or resolve the paradox of local equilibrium and global flow in DFT–PNP. One attempt is in refs 76 and Appendix C of ref 75.
Flows in Complex Mixtures
An important advantage of the methods presented here is their indifference to flow. The methods work at thermodynamic equilibrium and when flows are vigorous. Thus, calculations can be made in the nonequilibrium situations and mixed ionic solutions used nearly always in experimental work. Such calculations require some further numerical work because they involve many species of ions and ionic tracers (radioactive ions with properties identical to nonradioactive isotopes but present in trace amounts) that must be included appropriately in our Euler–Lagrange equations. Such calculations will allow the direct simulation of the experiments used historically to define single filing by ratios of unidirectional fluxes (e.g., which are estimated by the net fluxes of tracers) and most importantly to define the properties of transporters of every type, whether they transport species in the same direction or in opposite directions or in both. It is likely that some of the properties attributed to the interactions of ions with the channel protein actually occur between ions themselves.(73, 74) After all, interactions between ions had been ignored almost entirely in classical theories because ions were treated as ideal solutions(264, 287, 288) even when present in enormous concentrations.(9, 260) The effect of nonideality on the very definition of transporters remains to be investigated. This investigation can be done abstractly in general but much better if it can be done by realistic simulations of actual experimental setups using the PNP–steric model and other models that deal even more realistically with interactions.
Future Work: Inverse Problems
The distribution of permanent charge within the channel can be determined reliably from measurements of current voltage relations. Surprisingly, the inverse problem of determining the charge within the channel of the Nonner and Eisenberg model has been solved.(16-18) The sensitivity to noise and errors is small when the problem is solved by standard methods of inverse problems. The inverse problem of interest to biologists has a well-posed solution and can be used to determine the internal structure of the model channel from the kind of experimental information recorded in hundreds if not thousands of laboratories every day. The inverse problem for the PNP–steric model needs to be studied so that experiments can be designed to reveal properties of interest.
Future Work: One-Dimensional Models
The forward PNP problem has in some ways been more difficult to compute than the inverse problem because it has had to deal with the complex geometry of the channel protein. The inverse problem hides much of this geometrical complexity in effective 1D parameters. In fact, most numerical work assumes simple geometry for the ion channel and reduces the problem to one dimension. These papers all assumed that the pore diameter had some simple dependence on location, either parabolic or some kind of funnel shape that is easier to deal with analytically.(24, 25, 47, 48, 51, 88, 116, 143, 147, 148, 167, 289, 290) Although some 2D and 3D work has been reported in the past,(24, 35-40, 43-45, 49, 50, 281, 291-295) many results have not been as well converged as one might wish, and others simply were not checked as carefully as refs 26, 43, and 230 showed was necessary, as the semiconductor community had learned earlier (reviewed in ref 136; see ref 149). Very few computations have been reported using the real shape of channels, and even then the accuracy of the electrostatic treatment was not sufficient to be sure that important details were resolved (Claudio Berti, personal communication(139, 296)). Obviously, the difficulty is expected to be much greater when extending the traditional PNP equations to the modified ones of Eisenberg et al. and others.(26, 75-77, 230)
The difficulties in dealing with the full structural complexity of an ion channel should not be underestimated. The spatial resolution needed can put severe burdens on the memory bandwidth of even modern day computers. The relation of structures determined by the X-ray analysis of crystals to the spatial distribution of the mass density of each species, the spatial distribution of the (effective) diffusion coefficient (of each species), and the spatial distribution of polarization (i.e., induced dielectric) charge cannot be determined by presently available methods, but these distributions must be specified with precision in 3D calculations. The issues of temporal properties are hardly ever discussed, yet the polarization properties of electrolyte solutions change tremendously in the time range from biological function to atomic motion. It is not clear how these effects are involved in protein function or how to include them in models. It may be that these issues are less important in 1D models than in 3D models because the lower-dimensional models smooth over them in an appropriate way.
Time Dependence: Future Work
Future work needs to address each of the time-dependent phenomena to see what part of the classical properties of ion channels, studied in innumerable experimental papers, might arise from a model as simple as that used here. Obviously, many of those classical properties will involve conformational changes of the channel protein not described by our simple model. But just as obviously, those conformational changes will be coupled to ions in the channel by the electric field and probably by steric interactions as well, so everything must be analyzed together—ions, channel conformation, bathing solutions, ion flux, and current flow—as is usually the case in complex fluids flowing through complex spatial domains. Theory and simulations must allow everything to interact with everything else. They must not assume nothing interacts with nothing, as in ideal solutions.

Conclusions


Traditional PNP equations do not include the finite-size effect that is known to be significant in ionic solutions containing divalents and containing mixtures and even in pure monovalent solutions more concentrated than 50 mM. The concentration of ions in seawater, in and around cells, and inside channels is much higher than that. Therefore, classical PNP cannot describe the specific ion properties of bulk solutions such as seawater and the solutions in living systems, the plasmas of life. It cannot predict the ion-selectivity behavior of ion channels correctly. Here, we introduce the finite-size effect by treating ions and side chains as solid spheres and using hard sphere potentials to characterize this effect. Our work shows that selectivity is found in a simple 1D analysis and simulations.
Complex effects of changes in the repulsion parameters show a variety of states and depletion zones that are likely to be important in the functioning of channels and transporters. For example, the sudden appearance of a depletion zone, because of instability or a stochastic fluctuation, would surely gate a channel closed. If that gating happened on one side of a channel, then the properties on the other side would surely be affected. Everything interacts with everything else in these systems profoundly coupled by Coulomb and steric exclusion forces. If an exclusion zone moved from one side of a channel to another and then back and forth, then the channel protein could easily produce a reciprocating ping pong effect and mimic the alternating access states of transporters discovered with such wisdom and work by experimentalists who did not have the help of self-consistent models. The everything interacts with everything nature of the crowded charge environment inside a channel, or active sites,(2, 260, 297) makes such nonlinear interactions possible. It is not clear if the correlations included in our model are sufficient to produce ping pong effects or not: our model leaves out many forms of correlation, we are sad to say. It also remains to be seen whether biology actually uses such interactions at all. Alternating access could be produced in quite different ways, as most assume.
It is very important for the reader in the physical sciences to understand that complex systems of states and rates have been used by experimental biologists to characterize the function of hundreds of channels(178, 179, 239, 240) and transporters(298) studied by thousands of laboratories daily because of their medical and biological importance.
It is very important for the reader in the biological sciences to understand that an enormous wealth of living behavior could be controlled by the physical phenomena described here, as outputs of a self-consistent model and as solutions of a set of partial differential equations and boundary conditions, without invoking classical, vaguely defined effects. Those classical effects are more vitalistic than vital in many cases, in our view.
Early workers of some reputation in molecular biology, including Nobel Prize winners,(299-301) attributed the secret of life to allosteric interactions of chemical signals acting on proteins and then channels.(2) It is striking to the biologists among us that a self-consistent model of ions and side chains in channels produces strong interactions over long distances (i.e., more than 1 nm) without invoking the metaphors of vitalistic allostery. The calculations of a self-consistent variational theory of the energetics of complex fluids seem ready to replace the poetry of our ancestors.
Self-consistent theory is useful only because it can be evaluated with computers. Those computers in turn are possible only because of the successful treatment of complex physical interactions by self-consistent mathematics. It is amusing that physicists learned to use self-consistent mathematics to analyze (control and build) complex interacting systems of holes and electrons(302-304) during the same years that biologists used poetry to describe complex interacting systems of cations and anions.
Channels are nearly enzymes,(9, 297) and it is possible that the interactions described by models of the sort described here for channels underlie the complex interactions of a wide range of proteins that produce the special properties of life. Certainly, a theoretical and computational approach to biology and its molecules must allow everything to interact with everything else, instead of assuming that everything is ideal and nothing interacts with nothing.

The authors declare no competing financial interest.

Acknowledgment


The Taida Institute for Mathematical Sciences (TIMS) and National Center for Theoretical Sciences at Taipei (NCTS/Taipei) funded two meetings at which we met and developed these ideas. This paper is a result of their generous support, for which we are grateful. This work would not have been possible without the previous investigations of Dr. YunKyong Hyon. We are most grateful for his generosity in sharing all the details with us in innumerable communications and discussions. We thank Wolfgang Nonner and the referees for helpful suggestions that materially improved the article.

This research was sponsored in part by the National Science Council of Taiwan under grant no. NSC-100-2115-M-035-001 (T.L.H.). B.E. is supported by the Bard Endowed Chair of Rush University. He is particularly grateful for the administrative work he was not asked to do.

References


This article references 304 other publications.

  1. 1.
    Lin, T. C.; Eisenberg, B. To be submitted for publication, 2012.
  2. 2.
    Hille, B. Ionic Channels of Excitable Membranes, 3rd ed.; Sinauer Associates Inc.: Sunderland, England, 2001.
  3. 3.
    Ashcroft, F. M. Ion Channels and Disease; Academic Press: New York, 1999.
  4. 4.
    Neher, E.Ion Channels for Communication between and within Cells. In Nobel Lectures, Physiology or Medicine 1991–1995; Ringertz, N., Ed.; World Scientific: Singapore, 1997; p 10,

    Nobel Lecture, December 9, 1991.

  5. 5.
    Neher, E.; Sakmann, B. Nature 1976, 260, 799
  6. 6.
    Sakmann, B.; Neher, E. Single Channel Recording, 2nd ed.; Plenum: New York, 1995.
  7. 7.
    Eisenberg, B. ( 2012, available on arXiv as http://arxiv.org/abs/1206.6490.
  8. 8.
    Eisenberg, B. Fluctuations Noise Lett. 2012, 11, 76
  9. 9.
    Eisenberg, B.Crowded Charges in Ion Channels. Advances in Chemical Physics; John Wiley & Sons: New York, 2011; p 77.
  10. 10.
    Eisenberg, B. J. Phys. Chem. C 2010, 114, 20719
  11. 11.
    Eisenberg, R. S.Atomic Biology, Electrostatics and Ionic Channels. In New Developments and Theoretical Studies of Proteins; Elber, R., Ed.; World Scientific: Philadelphia, 1996; Vol. 7, p 269.
  12. 12.
    Eisenberg, R. S. J. Membr. Biol. 1996, 150, 1
  13. 13.
    Warshel, A. Proc. Natl. Acad. Sci. U.S.A. 1978, 75, 5250
  14. 14.
    Warshel, A.; Russell, S. T. Q. Rev. Biophys. 1984, 17, 283
  15. 15.
    Warshel, A.; Sharma, P. K.; Chu, Z. T.; Aqvist, J. Biochemistry 2007, 46, 1466
  16. 16.
    Arning, K. Mathematical Modelling and Simulation of Ion Channels; Johannes Kepler University Linz, 2009.
  17. 17.
    Arning, K.; Burger, M.; Eisenberg, R. S.; Engl, H. W.; He, L. PAMM 2007, 7, 1120801
  18. 18.
    Burger, M.; Eisenberg, R. S.; Engl, H. SIAM J. Appl. Math. 2007, 67, 960
  19. 19.
    Engl, H. W.; Hanke, M.; Neubauer, A. Regularization of Inverse Problems: Kluwer: Dordrecht, The Netherlands, 2000.
  20. 20.
    Kaipio, J.; Somersalo, E. Statistical and Computational Inverse Problems: Springer: New York, 2005.
  21. 21.
    Bazant, M. Z.; Thornton, K.; Ajdari, A. Phys. Rev. E 2004, 70, 021506
  22. 22.
    Coalson, R. D.; Kurnikova, M. G. IEEE Trans. Nanobiosci. 2005, 4, 81
  23. 23.
    Eisenberg, R.; Chen, D. Biophys. J. 1993, 64, A22
  24. 24.
    Liu, W.; Wang, B. J. Dynam. Differ. Equations 2010, 22, 413
  25. 25.
    Singer, A.; Norbury, J. SIAM J. Appl. Math. 2009, 70, 949
  26. 26.
    Zheng, Q.; Wei, G.-W. J. Chem. Phys. 2011, 134, 194101
  27. 27.
    Chen, D. P.; Lear, J.; Eisenberg, R. S. Biophys. J. 1997, 72, 97
  28. 28.
    Dieckmann, G. R.; Lear, J. D.; Zhong, Q.; Klein, M. L.; DeGrado, W. F.; Sharp, K. A. Biophys. J. 1999, 76, 618
  29. 29.
    Chen, D. P.; Nonner, W.; Eisenberg, R. S. Biophys. J. 1995, 68, A370
  30. 30.
    Chen, D.; Xu, L.; Tripathy, A.; Meissner, G.; Eisenberg, R. Biophys. J. 1997, 73, 1337
  31. 31.
    Allen, T. W.; Kuyucak, S.; Chung, S. H. Biophys. J. 1999, 77, 2502
  32. 32.
    Chen, D.; Xu, L.; Tripathy, A.; Meissner, G.; Eisenberg, B. Biophys. J. 1999, 76, 1346
  33. 33.
    Corry, B.; Kuyucak, S.; Chung, S. H. J. Gen. Physiol. 1999, 114, 597
  34. 34.
    Eisenberg, R. S. J. Membr. Biol. 1999, 171, 1
  35. 35.
    Hollerbach, U.; Chen, D.; Nonner, W.; Eisenberg, B. Biophys. J. 1999, 76, A205
  36. 36.
    Kurnikova, M. G.; Coalson, R. D.; Graf, P.; Nitzan, A. Biophys. J. 1999, 76, 642
  37. 37.
    Corry, B.; Kuyucak, S.; Chung, S. H. Biophys. J. 2000, 78, 2364
  38. 38.
    Graf, P.; Nitzan, A.; Kurnikova, M. G.; Coalson, R. D. J. Phys. Chem. B 2000, 104, 12324
  39. 39.
    Hollerbach, U.; Chen, D. P.; Busath, D. D.; Eisenberg, B. Langmuir 2000, 16, 5509
  40. 40.
    Moy, G.; Corry, B.; Kuyucak, S.; Chung, S. H. Biophys. J. 2000, 78, 2349
  41. 41.
    Im, W.; Roux, B. Biophys. J. 2001, 115, 4850
  42. 42.
    Edwards, S.; Corry, B.; Kuyucak, S.; Chung, S. H. Biophys. J. 2002, 83, 1348
  43. 43.
    Hollerbach, U.; Chen, D.-P.; Eisenberg, R. S. J. Comput. Sci. 2002, 16, 373
  44. 44.
    Hollerbach, U.; Eisenberg, R. Langmuir 2002, 18, 3262
  45. 45.
    Im, W.; Roux, B. J. Mol. Biol. 2002, 322, 851
  46. 46.
    Im, W.; Roux, B. J. Mol. Biol. 2002, 319, 1177
  47. 47.
    van der Straaten, T. A.; Tang, J.; Eisenberg, R. S.; Ravaioli, U.; Aluru, N. R. J. Comput. Electron. 2002, 1, 335
  48. 48.
    van der Straaten, T. A.; Tang, J. M.; Eisenberg, R. S.; Ravaioli, U.; Aluru, N.; Varma, S.; E., J. Biophys. J. 2002, 82, 207a
  49. 49.
    Corry, B.; Kuyucak, S.; Chung, S. H. Biophys. J. 2003, 84, 3594
  50. 50.
    Mamonov, A. B.; Coalson, R. D.; Nitzan, A.; Kurnikova, M. G. Biophys. J. 2003, 84, 3646
  51. 51.
    Nadler, B.; Hollerbach, U.; Eisenberg, R. S. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2003, 68, 021905
  52. 52.
    Eisenberg, B., available on arXiv as http://arxiv.org/abs/1207.4737 , 2012.
  53. 53.
    Eisenberg, B. Biophys. Chem. 2003, 100, 507
  54. 54.
    Giri, J.; Fonseca, J. E.; Boda, D.; Henderson, D.; Eisenberg, B. Phys. Biol. 2011, 8, 026004
  55. 55.
    Ganguly, P.; Mukherji, D.; Junghans, C.; van der Vegt, N. F. A. J. Chem. Theory Comput. 2012, 8, 1802
  56. 56.
    Molina, J. J.; Dufreche, J.-F.; Salanne, M.; Bernard, O.; Turq, P. J. Chem. Phys. 2011, 135, 234509
  57. 57.
    Fraenkel, D. J. Phys. Chem. B 2010, 115, 557
  58. 58.
    Kunz, W.; Neueder, R.An Attempt at an Overview. In Specific Ion Effects; Kunz, W., Ed.; World Scientific: Singapore, 2009; p 11.
  59. 59.
    Kunz, W. Specific Ion Effects; World Scientific: Singapore, 2009.
  60. 60.
    Kontogeorgis, G. M.; Folas, G. K. Thermodynamic Models for Industrial Applications: From Classical and Advanced Mixing Rules to Association Theories; Wiley: Chichester, U.K., 2009.
  61. 61.
    Abbas, Z.; Ahlberg, E.; Nordholm, S. J. Phys. Chem. B 2009, 113, 5905
  62. 62.
    Lee, L. L. Molecular Thermodynamics of Electrolyte Solutions; World Scientific: Singapore, 2008.
  63. 63.
    Jungwirth, P.; Winter, B. Annu. Rev. Phys. Chem. 2008, 59, 343
  64. 64.
    Jungwirth, P.; Finlayson-Pitts, B. J.; Tobias, D. J. Chem. Rev. 2006, 106, 1137
  65. 65.
    Ben-Naim, A. Molecular Theory of Solutions; Oxford University Press: New York, 2006.
  66. 66.
    Fawcett, W. R. Liquids, Solutions, and Interfaces: From Classical Macroscopic Descriptions to Modern Microscopic Details; Oxford University Press: New York, 2004.
  67. 67.
    Barthel, J.; Krienke, H.; Kunz, W. Physical Chemistry of Electrolyte Solutions: Modern Aspects; Springer: New York, 1998.
  68. 68.
    Durand-Vidal, S.; Turq, P.; Bernard, O.; Treiner, C.; Blum, L. Physica A 1996, 231, 123
  69. 69.
    Pitzer, K. S. Thermodynamics, 3rd ed.; McGraw Hill: New York, 1995.
  70. 70.
    Xu, Z.; Cai, W. SIAM Rev. 2011, 53, 683
  71. 71.
    Chazalviel, J.-N. Coulomb Screening by Mobile Charges; Birkhäuser: New York, 1999.
  72. 72.
    Eisenberg, B. Trans. Faraday Soc. 2012, doi: 10.1039/C2FD20066J.
  73. 73.
    Eisenberg, B., posted on arXiv.org under paper ID arXiv:1105.0184v1, 2011.
  74. 74.
    Eisenberg, B. Chem. Phys. Lett. 2011, 511, 1
  75. 75.
    Eisenberg, B.; Hyon, Y.; Liu, C. J. Chem. Phys. 2010, 133, 104104
  76. 76.
    Hyon, Y.; Eisenberg, B.; Liu, C. Commun. Math. Sci. 2011, 9, 459
  77. 77.
    Mori, Y.; Liu, C.; Eisenberg, R. S. Physica D 2011, 240, 1835
  78. 78.
    Hyon, Y.; Fonseca, J. E.; Eisenberg, B.; Liu, C. Discrete Continuous Dynamical Syst., Ser. B 2012, 17, 2725
  79. 79.
    Zhang, J.; Gong, X.; Liu, C.; Wen, W.; Sheng, P. Phys. Rev. Lett. 2008, 101, 194503
  80. 80.
    Doi, M. J. Phys. Soc. Jpn. 2009, 78, 052001
  81. 81.
    Liu, C.An Introduction of Elastic Complex Fluids: An Energetic Variational Approach. In Multi-scale Phenomena in Complex Fluids: Modeling, Analysis and Numerical Simulations; Hou, T. Y.; Liu, C.; Liu, J.-G, Eds.; World Scientific: Singapore, 2009.
  82. 82.
    Lin, F.-H.; Liu, C.; Zhang, P. Commun. Pure Appl. Math. 2007, 60, 838
  83. 83.
    Lin, F.-H.; Liu, C.; Zhang, P. Commun. Pure Appl. Math. 2005, 58, 1437
  84. 84.
    Sheng, P.; Zhang, J.; Liu, C. Prog. Theor. Phys. 2008, 131, Suppl. 175 .
  85. 85.
    Larson, R. G. The Structure and Rheology of Complex Fluids; Oxford University Press: New York, 1995.
  86. 86.
    Buyukdagli, S.; Manghi, M.; Palmeri, J. Phys. Rev. E 2010, 81, 041601
  87. 87.
    Nonner, W.; Chen, D. P.; Eisenberg, B. Biophys. J. 1998, 74, 2327
  88. 88.
    Nonner, W.; Eisenberg, B. Biophys. J. 1998, 75, 1287
  89. 89.
    Nonner, W.; Catacuzzeno, L.; Eisenberg, B. Biophys. J. 2000, 79, 1976
  90. 90.
    Boda, D.; Gillespie, D.; Nonner, W.; Henderson, D.; Eisenberg, B. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2004, 69, 046702
  91. 91.
    Boda, D.; Busath, D.; Eisenberg, B.; Henderson, D.; Nonner, W. Phys. Chem. Chem. Phys. 2002, 4, 5154
  92. 92.
    Boda, D.; Henderson, D.; Busath, D. Mol. Phys. 2002, 100, 2361
  93. 93.
    Eisenberg, B.Institute for Mathematics and its Applications; University of Minnesota; http://www.ima.umn.edu/, 2009.
  94. 94.
    Boda, D.; Giri, J.; Henderson, D.; Eisenberg, B.; Gillespie, D. J. Chem. Phys. 2011, 134, 055102
  95. 95.
    Boda, D.; Valisko, M.; Henderson, D.; Eisenberg, B.; Gillespie, D.; Nonner, W. J. Gen. Physiol. 2009, 133, 497
  96. 96.
    Boda, D.; Nonner, W.; Henderson, D.; Eisenberg, B.; Gillespie, D. Biophys. J. 2008, 94, 3486
  97. 97.
    Krauss, D.; Eisenberg, B.; Gillespie, D. Eur. Biophys. J. 2011, 40, 775
  98. 98.
    Dolphin, A. C. Br. J. Pharmacol. 2006, 147, S56
  99. 99.
    Boda, D.; Nonner, W.; Valisko, M.; Henderson, D.; Eisenberg, B.; Gillespie, D. Biophys. J. 2007, 93, 1960
  100. 100.
    Miedema, H.; Meter-Arkema, A.; Wierenga, J.; Tang, J.; Eisenberg, B.; Nonner, W.; Hektor, H.; Gillespie, D.; Meijberg, W. Biophys. J. 2004, 87, 3137
  101. 101.
    Miedema, H.; Vrouenraets, M.; Wierenga, J.; Gillespie, D.; Eisenberg, B.; Meijberg, W.; Nonner, W. Biophys. J. 2006, 91, 4392
  102. 102.
    Vrouenraets, M.; Wierenga, J.; Meijberg, W.; Miedema, H. Biophys. J. 2006, 90, 1202
  103. 103.
    Malasics, A.; Boda, D.; Valisko, M.; Henderson, D.; Gillespie, D. Biochim. Biophys. Acta 2010, 1798, 2013
  104. 104.
    Krauss, D.; Gillespie, D. Eur. Biophys. J. 2010, 39, 1513
  105. 105.
    Gillespie, D.; Giri, J.; Fill, M. Biophyis. J. 2009, 97, 2212
  106. 106.
    Gillespie, D.; Fill, M. Biophys. J. 2008, 95, 3706
  107. 107.
    Gillespie, D.; Boda, D.; He, Y.; Apel, P.; Siwy, Z. S. Biophys. J. 2008, 95, 609
  108. 108.
    Gillespie, D.; Boda, D. Biophys. J. 2008, 95, 2658
  109. 109.
    Gillespie, D. Biophys. J. 2008, 94, 1169
  110. 110.
    Gillespie, D.; Valisko, M.; Boda, D. J. Phys.: Condens. Matter 2005, 17, 6609
  111. 111.
    Gillespie, D.; Nonner, W.; Eisenberg, R. S. Phys. Rev. E 2003, 68, 0313503
  112. 112.
    Gillespie, D.; Nonner, W.; Eisenberg, R. S. J. Phys.: Condens. Matter 2002, 14, 12129
  113. 113.
    Chen, D. P.Nonequilibrium Thermodynamics of Transports in Ion Channels. In Progress of Cell Research: Towards Molecular Biophysics of Ion Channels; Sokabe, M.; Auerbach, A.; Sigworth, F., Eds.; Elsevier: Amsterdam, 1997; p 269.
  114. 114.
    Chen, D.; Xu, L.; Tripathy, A.; Meissner, G.; Eisenberg, R. Biophys. J. 1997, 72, A108
  115. 115.
    Hyon, Y.; Kwak, D. Y.; Liu, C. Discrete Continuous Dynamical Syst., Ser. A 2010, 26, 1291
  116. 116.
    Barcilon, V.; Chen, D. P.; Eisenberg, R. S. SIAM J. Appl. Math. 1992, 52, 1405
  117. 117.
    Gillespie, D.; Eisenberg, R. S. Eur. Biophys. J. 2002, 31, 454
  118. 118.
    Gillespie, D.; Nonner, W.; Eisenberg, R. S. Biophys. J. Abstr. 2002, 84, 67a
  119. 119.
    Gardner, C. L.; Nonner, W.; Eisenberg, R. S. J. Comput. Electron. 2004, 3, 25
  120. 120.
    Aguilella-Arzo, M.; Aguilella, V.; Eisenberg, R. S. Eur. Biophys. J. 2005, 34, 314
  121. 121.
    Luchinsky, D. G.; Tindjong, R.; Kaufman, I.; McClintock, P. V. E.; Eisenberg, R. S. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2009, 80, 021925
  122. 122.
    Doi, M. J. Phys.: Condens. Matter 2011, 23, 284118
  123. 123.
    Hou, T. Y.; Liu, C.; Liu, J.-G. Multi-scale Phenomena in Complex Fluids: Modeling, Analysis and Numerical Simulations; World Scientific: Singapore, 2009.
  124. 124.
    Hyon, Y.; Carrillo, J. A.; Du, Q.; Liu, C. Kinetic Relat. Models 2008, 1, 171
  125. 125.
    Hyon, Y.; Du, Q.; Liu, C. J. Comput. Theor. Nanosci. 2010, 7, 756
  126. 126.
    Ryham, R. J.An Energetic Variational Approach to Mathematical Moldeling of Charged Fluids, Charge Phases, Simulation and Well Posedness. Ph.D. Thesis, The Pennsylvania State University, 2006.
  127. 127.
    Xu, X.; Liu, C.; Qian, T. Commun. Math. Sci. 2012, 10, 1027
  128. 128.
    Ryham, R.; Cohen, F.; Eisenberg, R. S. Commun. Math. Sci. 2012, not supplied.
  129. 129.
    Dan, B.-Y.; Andelman, D.; Harries, D.; Podgornik, R. J. Phys.: Condens. Matter 2009, 21, 424106
  130. 130.
    Malasics, A.; Gillespie, D.; Boda, D. J. Chem. Phys. 2008, 128, 124102
  131. 131.
    Bruna, M.; Chapman, S. J. Phys. Rev. E 2012, 85, 011103
  132. 132.
    Streetman, B. G. Solid State Electronic Devices, 4th ed.; Prentice Hall: Englewood Cliffs, NJ, 1972.
  133. 133.
    Sze, S. M. Physics of Semiconductor Devices; John Wiley & Sons: New York, 1981.
  134. 134.
    Selberherr, S. Analysis and Simulation of Semiconductor Devices; Springer-Verlag: New York, 1984.
  135. 135.
    Markowich, P. A.; Ringhofer, C. A.; Schmeiser, C. Semiconductor Equations; Springer-Verlag: New York, 1990.
  136. 136.
    Jerome, J. W. Analysis of Charge Transport: Mathematical Theory and Approximation of Semiconductor Models; Springer-Verlag: New York, 1995.
  137. 137.
    Howe, R. T.; Sodini, C. G. Microelectronics: An Integrated Approach; Prentice Hall: Upper Saddle River, NJ, 1997.
  138. 138.
    Hess, K. Advanced Theory of Semiconductor Devices; IEEE Press: New York, 2000.
  139. 139.
    Berti, C.; Gillespie, D.; Eisenberg, R. S.; Fiegna, C. Nanoscale Res. Lett. 2012, not supplied.
  140. 140.
    Eisenberg, B. Pro. Natl. Acad. Sci. U.S.A. 2008, 105, 6211
  141. 141.
    Eisenberg, B. available on http://arxiv.org/ as q-bio/0506016v2, 2005.
  142. 142.
    Eisenberg, B. J. Comput. Electron. 2003, 2, 245
  143. 143.
    Barcilon, V.; Chen, D.-P.; Eisenberg, R. S.; Jerome, J. W. SIAM J. Appl. Math. 1997, 57, 631
  144. 144.
    Chen, D.; Eisenberg, R.; Jerome, J.; Shu, C. Biophys. J. 1995, 69, 2304
  145. 145.
    Chen, D. P.; Eisenberg, R. S. Biophys. J. 1993, 64, 1405
  146. 146.
    Chen, D. P.; Barcilon, V.; Eisenberg, R. S. Biophys. J. 1992, 61, 1372
  147. 147.
    Abaid, N.; Eisenberg, R. S.; Liu, W. SIAM J. Appl. Dynam. Syst. 2008, 7, 1507
  148. 148.
    Eisenberg, B.; Liu, W. SIAM J. Math. Anal. 2007, 38, 1932
  149. 149.
    Vasileska, D.; Goodnick, S. M.; Klimeck, G. Computational Electronics: Semiclassical and Quantum Device Modeling and Simulation; CRC Press: New York, 2010.
  150. 150.
    Kong, C. L. J. Chem. Phys. 1973, 59, 2464
  151. 151.
    Cazade, P. A.; Dweik, J.; Coasne, B.; Henn, F.; Palmeri, J. J. Phys. Chem. C 2010, 114, 12245
  152. 152.
    Lakatos, G.; Patey, G. N. J. Chem. Phys. 2007, 126, 024703
  153. 153.
    Horng, T.-L.; Teng, C.-H. J. Comput. Phys. 2012, 231, 2498
  154. 154.
    Temes, G. C.; Barcilon, V.; Marshall, F. C. Proc. IEEE 1973, 61, 196
  155. 155.
    Lee, C. C.; Hyon, Y.; Lin, T. C.; Liu, C. To be submitted for publication, 2012.
  156. 156.
    Boda, D.; Henderson, D.; Eisenberg, B.; Gillespie, D. J. Chem. Phys. 2011, 135, 064105
  157. 157.
    Singer, A.; Gillespie, D.; Norbury, J.; Eisenberg, R. S. Eur. J. Appl. Math. 2008, 19, 541
  158. 158.
    Howard, J. J.; Perkyns, J. S.; Pettitt, B. M. J. Phys. Chem. B 2010, 114, 6074
  159. 159.
    Varsos, K.; Luntz, J.; Welsh, M.; Sarabandi, K. Proc. IEEE 2011, 99, 2112
  160. 160.
    Fiedziuszko, S. J.; Hunter, I. C.; Itoh, T.; Kobayashi, Y.; Nishikawa, T.; Stitzer, S.; Wakino, K. IEEE Trans. Microwave Theory Tech. 2002, 50, 706
  161. 161.
    Helfferich, F. Ion Exchange; Dover: New York, 1995.
  162. 162.
    Shur, M. Physics of Semiconductor Devices; Prentice Hall: New York, 1990.
  163. 163.
    Collins, A.; Somlyo, A. V.; Hilgemann, D. W. J. Physiol. 1992, 454, 27
  164. 164.
    Hilgemann, D. W.; Collins, A.; Matsuoka, S. J. Gen. Physiol. 1992, 100, 933
  165. 165.
    Hilgemann, D. W.; Matsuoka, S.; Nagel, G. A.; Collins, A. J. Gen. Physiol. 1992, 100, 905
  166. 166.
    Matsuoka, S.; Hilgemann, D. W. J. Gen. Physiol. 1992, 100, 963
  167. 167.
    Schuss, Z.; Nadler, B.; Eisenberg, R. S. Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys. 2001, 64, 036116
  168. 168.
    Schuss, Z.; Nadler, B.; Singer, A.; Eisenberg, R.A PDE Formulation of Non-Equilibrium Statistical Mechanics for Ionic Permeation. In Unsolved problems of noise and fluctuations : UPoN 2002 : Third International Conference of Unsolved Problems of Noise and Fluctuations in Physics, Biology, and High Technology, AIP Conference Proceedings, Washington, DC, September 3–6, 2002; Bezrukov, S. M., Ed.; American Institute of Physics: Melville, NY, 2003; No. 665.
  169. 169.
    Hünenberger, P. H.; Reif, M. Single-Ion Solvation; RSC Publishing: Cambridge U.K., 2011.
  170. 170.
    Ben-Naim, A. Molecular Theory of Water and Aqueous Solutions Part II: The Role of Water in Protein Folding, Self-Assembly and Molecular Recognition; World Scientific: Singapore, 2011.
  171. 171.
    Fraenkel, D. Mol. Phys. 2010, 108, 1435
  172. 172.
    Hodgkin, A. L.; Huxley, A. F. J. Physiol. 1952, 116, 449
  173. 173.
    Hodgkin, A. L.; Huxley, A. F. J. Physiol. 1952, 116, 473
  174. 174.
    Mullins, L. J. J. Gen. Physiol. 1968, 52, 555
  175. 175.
    Mullins, L. J. J. Gen. Physiol. 1968, 52, 550
  176. 176.
    Moore, J. W.; Blaustein, M. P.; Anderson, N. C.; Narahashi, T. J. Gen. Physiol. 1967, 50, 1401
  177. 177.
    Moore, J. W.; Narahashi, T.; Anderson, N. C.; Blaustein, M. P.; Watanabe, A.; Tasaki, I.; Singer, I.; Lerman, L. Science 1967, 157, 220
  178. 178.
    Conley, E. C. The Ion Channel Facts Book. I. Extracellular Ligand-Gated Channels; Academic Press: New York, 1996; Vol. 1.
  179. 179.
    Conley, E. C. The Ion Channel Facts Book. II. Intracellular Ligand-Gated Channels; Academic Press: New York, 1996; Vol. 2.
  180. 180.
    Unwin, N. J. Mol. Biol. 2005, 346, 967
  181. 181.
    Kraus, C. A. Bull. Am. Math. Soc. 1938, 44, 361
  182. 182.
    Harned, H. S.; Owen, B. B. The Physical Chemistry of Electrolytic Solutions, 3rd ed.; Reinhold Publishing Corporation: New York, 1958.
  183. 183.
    Friedman, H. L. Ionic Solution Theory; Interscience Publishers: New York, 1962.
  184. 184.
    Rice, S. A.; Gray, P. Statistical Mechanics of Simple Fluids; Wiley-Interscience: New York, 1965.
  185. 185.
    Barlow, C. A., Jr.; Macdonald, J. R.Theory of Discreteness of Charge Effects in the Electrolyte Compact Double Layer. In Advances in Electrochemistry and Electrochemical Engineering; Delahay, P., Ed.; Interscience Publishers: New York, 1967; Vol. 6.
  186. 186.
    Resibois, P. M. V. Electrolyte Theory; Harper & Row: New York, 1968.
  187. 187.
    Conway, B. E. Electrochemical Data; Greenwood Press Publishers: Westport CT, 1969.
  188. 188.
    Barker, J.; Henderson, D. Rev. Mod. Phys. 1976, 48, 587
  189. 189.
    Ben-Naim, A. Inversion of the Kirkwood Buff Theory of Solutions: Application to the Water Ethanol System; American Institute of Physics: Melville, NY, 1977; Vol. 67.
  190. 190.
    Pytkowicz, R. M. Activity Coefficients in Electrolyte Solutions; CRC: Boca Raton FL, 1979; Vol. 1.
  191. 191.
    Torrie, G. M.; Valleau, A. J. Phys. Chem. 1982, 86, 3251
  192. 192.
    Hovarth, A. L. Handbook of Aqueous Electrolyte Solutions: Physical Properties, Estimation, and Correlation Methods; Ellis Horwood: New York, 1985.
  193. 193.
    Patwardhan, V. S.; Kumar, A. AIChE J. 1986, 32, 1429
  194. 194.
    Stell, G.; Joslin, C. G. Biophys. J. 1986, 50, 855
  195. 195.
    Zemaitis, J. F., Jr.; Clark, D. M.; Rafal, M.; Scrivner, N. C. Handbook of Aqueous Electrolyte Thermodynamics; Design Institute for Physical Property Data sponsored by the American Institute of Chemical Engineers: New York, 1986.
  196. 196.
    Lee, L. L. Molecular Thermodynamics of Nonideal Fluids: Butterworth-Heinemann: New York, 1988.
  197. 197.
    Pitzer, K. S. Activity Coefficients in Electrolyte Solutions; CRC Press: Boca Raton, FL, 1991.
  198. 198.
    Cabezas, H.; O’Connell, J. P. Ind. Eng. Chem. Res. 1993, 32, 2892
  199. 199.
    Patwardhan, V. S.; Kumar, A. AIChE J. 1993, 39, 711
  200. 200.
    Chhih, A.; Bernard, O.; Barthel, J. M. G.; Blum, L. Ber. Bunsen-Ges. Phys. Chem. 1994, 98, 1516
  201. 201.
    Hunenberger, P. H.; McCammon, J. A. J. Chem. Phys. 1999, 110, 1856
  202. 202.
    Baker, N. A.; Hunenberger, P. H.; McCammon, J. A. J. Chem. Phys. 2000, 113, 2510
  203. 203.
    Durand-Vidal, S.; Simonin, J.-P.; Turq, P. Electrolytes at Interfaces; Kluwer: Boston, 2000.
  204. 204.
    Kumar, P. V.; Maroncelli, M. J. Chem. Phys. 2000, 112, 5370
  205. 205.
    Heinz, T. N.; van Gunsteren, W. F.; Hunenberger, P. H. J. Chem. Phys. 2001, 115, 1125
  206. 206.
    Abbas, Z.; Gunnarsson, M.; Ahlberg, E.; Nordholm, S. J. Phys. Chem. B 2002, 106, 1403
  207. 207.
    Laidler, K. J.; Meiser, J. H.; Sanctuary, B. C. Physical Chemistry, 4th ed.; BrooksCole: Belmont, CA, 2003.
  208. 208.
    Kastenholz, M. A.; Hunenberger, P. H. J. Chem. Phys. 2006, 124, 124106
  209. 209.
    Kornyshev, A. A. J. Phys. Chem. B 2007, 111, 5545
  210. 210.
    Che, J.; Dzubiella, J.; Li, B.; McCammon, J. A. J. Phys. Chem. B 2008, 112, 3058
  211. 211.
    Henderson, D.; Boda, D. Phys. Chem. Chem. Phys. 2009, 11, 3822
  212. 212.
    Horinek, D.; Mamatkulov, S.; Netz, R. J. Chem. Phys. 2009, 130, 124507
  213. 213.
    Ibarra-Armenta, J. G.; Martin-Molina, A.; Quesada-Perez, M. Phys. Chem. Chem. Phys. 2009, 11, 309
  214. 214.
    JaneCek, J.; Netz, R. R. J. Chem. Phys. 2009, 130, 074502
  215. 215.
    Kalcher, I.; Dzubiella, J. J. Chem. Phys. 2009, 130, 134507
  216. 216.
    Li, B. Nonlinearity 2009, 22, 811
  217. 217.
    Li, B. SIAM J. Math. Anal. 2009, 40, 2536
  218. 218.
    Luo, Y.; Roux, B. J. Phys. Chem. Lett. 2009, 1, 183
  219. 219.
    Maginn, E. J. AIChE J. 2009, 55, 1304
  220. 220.
    Kalcher, I.; Schulz, J. C. F.; Dzubiella, J. Phys. Rev. Lett. 2010, 104, 097802
  221. 221.
    Kalyuzhnyi, Y. V.; Vlachy, V.; Dill, K. A. Phys. Chem. Chem. Phys. 2010, 12, 6260
  222. 222.
    Lipparini, F.; Scalmani, G.; Mennucci, B.; Cances, E.; Caricato, M.; Frisch, M. J. J. Chem. Phys. 2010, 133, 014106
  223. 223.
    Sala, J.; Guardia, E.; Marti, J. J. Chem. Phys. 2010, 132, 214505
  224. 224.
    Vincze, J.; Valisko, M.; Boda, D. J. Chem. Phys. 2010, 133, 154507
  225. 225.
    Yu, H.; Whitfield, T. W.; Harder, E.; Lamoureux, G.; Vorobyov, I.; Anisimov, V. M.; MacKerell, A. D.; Roux, B. J. Chem. Theory Comput. 2010, 6, 774
  226. 226.
    Zhang, C.; Raugei, S.; Eisenberg, B.; Carloni, P. J. Chem. Theory Comput. 2010, 6, 2167
  227. 227.
    Bazant, M. Z.; Storey, B. D.; Kornyshev, A. A. Phys. Rev. Lett. 2011, 106, 046102
  228. 228.
    Gee, M. B.; Cox, N. R.; Jiao, Y.; Bentenitis, N.; Weerasinghe, S.; Smith, P. E. J. Chem. Theory Comput. 2011, not supplied.
  229. 229.
    Xiao, T.; Song, X. J. Chem. Phys. 2011, 135, 104104
  230. 230.
    Zheng, Q.; Chen, D.; Wei, G.-W. J. Comput. Phys. 2011, 230, 5239
  231. 231.
    Boron, W.; Boulpaep, E. Medical Physiology; Saunders: New York, 2008.
  232. 232.
    Heras, J. A. Am. J. Phys. 2007, 75, 652
  233. 233.
    Heras, J. A. Am. J. Phys. 2008, 76, 101
  234. 234.
    Heras, J. A. Am. J. Phys. 2011, 79, 409
  235. 235.
    Hladky, S. B.; Haydon, D. A. Nature 1970, 225, 451
  236. 236.
    Sigworth, F. J.; Neher, E. Nature 1980, 287, 447
  237. 237.
    Hamill, O. P.; Marty, A.; Neher, E.; Sakmann, B.; Sigworth, F. J. Pflugers Arch. 1981, 391, 85
  238. 238.
    FitzHugh, R. Math. Biosci. 1983, 64, 75
  239. 239.
    Conley, E. C.; Brammar, W. J. The Ion Channel Facts Book IV: Voltage Gated Channels; Academic Press: New York, 1999.
  240. 240.
    Conley, E. C.; Brammar, W. J. The Ion Channel Facts Book III: Inward Rectifier and Intercellular Channels: Academic Press: New York, 2000.
  241. 241.
    Cherny, V. V.; Murphy, R.; Sokolov, V.; Levis, R. A.; DeCoursey, T. E. J. Gen. Physiol. 2003, 121, 615
  242. 242.
    Levis, R. A.; Rae, J. L.; Iverson, L.; Rudy, B. Methods Enzymol. 1992, 207, 66
  243. 243.
    Levis, R. A.; Rae, J. L. Methods Enzymol. 1998, 293, 218
  244. 244.
    Tang, J.; Levis, R.; Lynn, K.; Eisenberg, B. Biophys. J. 1995, 68, A145
  245. 245.
    Miodownik-Aisenberg, J. Gating Kinetics of a Calcium-Activated Potassium Channel Studied at Subzero Temperatures; University of Miami Medical School, 1995.
  246. 246.
    Miller, C. Ion Channel Reconstitution; Plenum Press: New York, 1986.
  247. 247.
    Miceli, F.; Vargas, E.; Bezanilla, F.; Taglialatela, M. Biophys. J. 2012, 102, 1372
  248. 248.
    Dekel, N.; Priest, M. F.; Parnas, H.; Parnas, I.; Bezanilla, F. Proc. Natl. Acad. Sci. U.S.A. 2012, 109, 285
  249. 249.
    Vargas, E.; Bezanilla, F.; Roux, B. Neuron 2011, 72, 713
  250. 250.
    Lacroix, J. J.; Bezanilla, F. Proc. Natl. Acad. Sci. U.S.A. 2011, 108, 6444
  251. 251.
    Lacroix, J.; Halaszovich, C. R.; Schreiber, D. N.; Leitner, M. G.; Bezanilla, F.; Oliver, D.; Villalba-Galea, C. A. J. Biol. Chem. 2011, 286, 17945
  252. 252.
    Vandenberg, C. A.; Bezanilla, F. Biophys. J. 1991, 60, 1511
  253. 253.
    Vandenberg, C. A.; Bezanilla, F. Biophys. J. 1991, 60, 1499
  254. 254.
    Rae, J. L.; Levis, R. A. Pflugers Arch. 1992, 420, 618
  255. 255.
    Levis, R. A.; Rae, J. L. Biophys. J. 1993, 65, 1666
  256. 256.
    Levis, R. A.; Rae, J. L.Technology of Patch Clamp Recording Electrodes. In Patch-Clamp Applications and Protocols; Walz, W.; Boulton, A.; Baker, G., Eds.; Humana Press: Totowa, NJ, 1995.
  257. 257.
    Doyle, D. A.; Cabral, J. M.; Pfuetzner, R. A.; Kuo, A.; Gulbis, J. M.; Cohen, S. L.; Chait, B. T.; MacKinnon, R. Science 1998, 280, 69
  258. 258.
    MacKinnon, R. Angew. Chem., Int. Ed. 2004, 43, 4265
  259. 259.
    Schmidt, D.; Cross, S. R.; MacKinnon, R. J. Mol. Biol. 2009, 390, 902
  260. 260.
    Jimenez-Morales, D.; Liang, J.; Eisenberg, B. Eur. Biophys. J. 2012, 41, 449
  261. 261.
    Lutz, S. Science 2010, 329, 285
  262. 262.
    Ringe, D.; Petsko, G. A. Science 2008, 1428
  263. 263.
    Warshel, A. J. Biol. Chem. 1998, 273, 27035
  264. 264.
    Dixon, M.; Webb, E. C. Enzymes; Academic Press: New York, 1979.
  265. 265.
    Martin, P. A. Rev. Mod. Phys. 1988, 60, 1075
  266. 266.
    Henderson, J. R.Statistical Mechanical Sum Rules. In Fundamentals of Inhomogeneous Fluids; Henderson, D., Ed.; Marcel Dekker: New York, 1992; p 23.
  267. 267.
    Fuoss, R. M.; Kraus, C. A. J. Am. Chem. Soc. 1933, 55, 2387
  268. 268.
    Fuoss, R. M. Chem. Rev. 1935, 17, 27
  269. 269.
    Fuoss, R. M.; Onsager, L. Proc. Natl. Acad. Sci. U.S.A. 1955, 41, 274
  270. 270.
    Fuoss, R. M.; Onsager, L. J. Phys. Chem. B 1958, 62, 1339
  271. 271.
    Fuoss, R. M. J. Phys. Chem. 1959, 63, 633
  272. 272.
    Fuoss, R. M.; Onsager, L. J. Phys. Chem. 1962, 66, 1722
  273. 273.
    Fuoss, R. M.; Onsager, L. J. Phys. Chem. 1962, 66, 1722
  274. 274.
    Fuoss, R. M.; Onsager, L. J. Phys. Chem. 1963, 67, 621
  275. 275.
    Fuoss, R. M.; Onsager, L. J. Phys. Chem. 1963, 67, 628
  276. 276.
    Fuoss, R. M.; Onsager, L. J. Phys. Chem. 1964, 68, 1
  277. 277.
    Fuoss, R. M.; Onsager, L.; Skinner, J. F. J. Phys. Chem. 1965, 69, 2581
  278. 278.
    Justice, J.-C.Conductance of Electrolyte Solutions. In Thermondynamic and Transport Properties of Aqueous and Molten Electrolytes; Conway, B. E.; Bockris, J. O. M.; Yaeger, E., Eds.; Comprehensive Treatise of Electrochemistry; Plenum: New York, 1983; Vol. 5, p 223.
  279. 279.
    Csanyi, E.; Boda, D.; Gillespie, D.; Kristof, T. Biochim. Biophys. Acta 2012, 1818, 592
  280. 280.
    Boda, D.; Gillespie, D. J. Chem. Theory Comput. 2012, 8, 824
  281. 281.
    Rutkai, G.; Boda, D.; Kristóf, T. J. Phys. Chem. Lett. 2010, 1, 2179
  282. 282.
    Damocles. Damocles Web Site, IBM Research. http://www.research.ibm.com/DAMOCLES/home.html, 2007.
  283. 283.
    Lundstrom, M. Fundamentals of Carrier Transport, 2nd ed.; Addison-Wesley: New York, 2000.
  284. 284.
    Jacoboni, C.; Lugli, P. The Monte Carlo Method for Semiconductor Device Simulation; Springer-Verlag: New York, 1989.
  285. 285.
    Eisenberg, R. S.; Kłosek, M. M.; Schuss, Z. J. Chem. Phys. 1995, 102, 17671780
  286. 286.
    Kłosek, M. M. J. Stat. Phys. 1995, 79, 313
  287. 287.
    Segel, I. H. Enzyme Kinetics: Behavior and Analysis of Rapid Equilibrium and Steady State Enzyme Systems; Wiley: New York, 1993.
  288. 288.
    Tosteson, D. Membrane Transport: People and Ideas; American Physiological Society: Bethesda, MD, 1989.
  289. 289.
    van der Straaten, T. A.; Tang, J. M.; Ravaioli, U.; Eisenberg, R. S.; Aluru, N. R. J. Comput. Electron. 2003, 2, 29
  290. 290.
    Eisenberg, B. Acc. Chem. Res. 1998, 31, 117
  291. 291.
    Lu, B.; McCammon, J. A. Chem. Phys. Lett. 2008, 451, 282
  292. 292.
    Vora, T.; Corry, B.; Chung, S. H. Biochim. Biophys. Acta 2006, 1758, 730
  293. 293.
    Mamonov, A. B.; Kurnikova, M. G.; Coalson, R. D. Biophys. Chem. 2006, 124, 268
  294. 294.
    Corry, B.; Vora, T.; Chung, S. H. Biochim. Biophys. Acta 2005, 1711, 72
  295. 295.
    Cardenas, A. E.; Coalson, R. D.; Kurnikova, M. G. Biophys. J. 2000, 79.
  296. 296.
    Berti, C.; Gillespie, D.; Eisenberg, B.; Furini, S.; Fiegna, C. Biophys. J. 2011, 100, 158a
  297. 297.
    Eisenberg, R. S. J. Membr. Biol. 1990, 115, 1
  298. 298.
    Griffiths, J.; Sansom, C. The Transporter Facts Book; Academic Press: New York, 1997.
  299. 299.
    Jacob, F. The Logic of Life: A History or Heredity; Princeton University Press: Princeton, NJ, 1993.
  300. 300.
    Monod, J. Chance and Necessity: An Essay on the Natural Philosophy of Modern Biology; Vintage Books: New York, 1972.
  301. 301.
    Schrödinger, E. What Is Life?; Cambridge University Press: New York, 1992.
  302. 302.
    Shockley, W. Electrons and Holes in Semiconductors, with Applications to Transistor Electronics; Van Nostrand: New York, 1950.
  303. 303.
    Van Roosbroeck, W. Bell Syst. Tech. J. 1950, 29, 560
  304. 304.
    Gummel, H. K. IEEE Trans. Electron Devices 1964, ED-11, 445

Tools

SciFinder Links

SciFinder subscribers:  Click to sign in | Not a SciFinder subscriber? Learn more at www.cas.org

Explore by:


History

  • Published In Issue September 20, 2012
  • Article ASAPSeptember 10, 2012
  • Just Accepted ManuscriptAugust 17, 2012
  • Received: May 30, 2012
    Revised: July 31, 2012

Recommend & Share

  • Share on ACS NetworkACS Network
  • Add to FacebookFacebook
  • Tweet ThisTweet This
  • Add to CiteULikeCiteULike
  • Add to NewsvineNewsvine
  • Digg ThisDigg This
  • Add to DeliciousDelicious

Related Content

Other ACS content by these authors: