Electrodiffusion Kinetics of Ionic Transport in a Simple Membrane

Oct 28, 2013 - We employ numerical techniques for solving time-dependent full Poisson–Nernst–Planck (PNP) equations in 2D to analyze transient beh...
0 downloads 0 Views 2MB Size
Article pubs.acs.org/JPCB

Electrodiffusion Kinetics of Ionic Transport in a Simple Membrane Channel Ivan Valent,*,† Pavol Petrovič,† Pavel Neogrády,† Igor Schreiber,‡ and Miloš Marek‡ †

Department of Physical and Theoretical Chemistry, Faculty of Natural Sciences, Comenius University, Mlynská dolina, 842 15 Bratislava, Slovak Republic ‡ Department of Chemical Engineering, Institute of Chemical Technology, Prague, Technická 5, 166 28 Prague 6, Czech Republic W Web-Enhanced Feature * S Supporting Information *

ABSTRACT: We employ numerical techniques for solving timedependent full Poisson−Nernst−Planck (PNP) equations in 2D to analyze transient behavior of a simple ion channel subject to a sudden electric potential jump across the membrane (voltage clamp). Calculated spatiotemporal profiles of the ionic concentrations and electric potential show that two principal exponential processes can be distinguished in the electrodiffusion kinetics, in agreement with original Planck’s predictions. The initial fast process corresponds to the dielectric relaxation, while the steady state is approached in a second slower exponential process attributed to the nonlinear ionic redistribution. Effects of the model parameters such as the channel length, height of the potential step, boundary concentrations, permittivity of the channel interior, and ionic mobilities on electrodiffusion kinetics are studied. Numerical solutions are used to determine spatiotemporal profiles of the electric field, ionic fluxes, and both the conductive and displacement currents. We demonstrate that the displacement current is a significant transient component of the total electric current through the channel. The presented results provide additional information about the classical voltage-clamp problem and offer further physical insights into the mechanism of electrodiffusion. The used numerical approach can be readily extended to multi-ionic models with a more structured domain geometry in 2D or 3D, and it is directly applicable to other systems, such as synthetic nanopores, nanofluidic channels, and nanopipettes.

1. INTRODUCTION Transport of ions plays an important role in many functional mechanisms in biological systems. Since the pioneering work of Luigi Galvani in research of ‘bioelectricity’, evidence for action of ions in various phenomena observed in living matter has been put forward. For example, the identification of potassium and sodium ionic currents in electrical behavior of a nerve fiber1,2 was one of the milestones in the development of the field of electrophysiology. Other examples include the signaling role of calcium ions3 in muscle contraction, synaptic excitation, hormone secretion, regulation of enzyme activity and many other biological processes.4 Transport of ions is usually localized in particular areas of cells, and in many cases it is related to a membrane function. The driving force of such ionic movement is not restricted to the concentration gradient but involves also the electric potential difference present across the membrane or electrostatic interactions between ions and a charged ion-binding macromolecule or charged membrane surface. Movement of charged particles driven by gradients of both chemical and electrical (i.e., electrochemical) potentials is termed electrodiffusion. The concept of electrodiffusion based on the Nernst−Planck equations for ionic fluxes coupled to the Poisson equation expressing the relation between the gradient of the electric field and the charge density (hence usually © 2013 American Chemical Society

referred to as the NPP or PNP theory) is widely used as a useful theoretical tool for description and modeling of various phenomena and processes in both living and nonliving systems. The PNP theory extends to a variety of fields such as membrane biophysics,5 enzyme kinetics,6 colloid and interfacial science,7 electrochemistry and electroanalytical chemistry,8 nanofluidics,9 solid-state physics,10 or materials science.11 The mathematical formulation of the PNP theory is expressed by a set of partial differential equations (PDEs). Although the formulation is straightforward, to solve these nonlinear equations is a difficult task. The Poisson equation is frequently eliminated using approximations such as the electroneutrality condition, first used by Planck,12 or the constant field assumption applied by Goldman13 to avoid complications with finding a complete solution. Also, solutions are often restricted to the steady-state case, in contrast with few attempts to solve the time-dependent PNP equations. For example, the steady-state solutions were used to model ion transport through liquid junctions,14 protein channels,15,16 and synthetic nanopores.17,18 Early efforts to analyze time-dependReceived: July 27, 2013 Revised: October 16, 2013 Published: October 28, 2013 14283

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

ent electrodiffusion systems involve ion exchange kinetics19,20 and transport of ions through liquid junctions 21 and membranes.22−24 Fundamental contributions to solving the time-dependent PNP equations also came from the field of electrochemistry, mostly related to electrodiffusion dynamics of the electrical double layer.25−27 Examples of recent applications of the PNP theory to time-dependent problems include translocation of DNA through nanopores and nanochannels,28 or electrochemical methods such as cyclic voltammetry29 and electrochemical impedance spectroscopy.30 In addition, most of the models presented in the literature are spatially 1D, whereas time-dependent solutions of the PNP equations in higher spatial geometries are scarce.31 Our motivation is a perspective of modeling electrodiffusion dynamics (Section 2) in biological processes. This aim requires implementation of numerical techniques of solving the PNP equations for complex models with spatially 2D and 3D geometry. We start with a simple 2D model of a membrane channel (Section 3) and simulate distribution of the concentration of ions and electric potential in space and time after a voltage step. The problem is formulated to be uniform in the direction of the y coordinate so that the computational technique can be tuned and tested by comparison with the known results for 1D models (Section 4). In addition to the development of the methodology, the presented results (Section 5) bring further insight into the kinetics of charge separation in a solution of strong electrolyte in a membrane under the voltage-clamp conditions. Although the 1D problem under study is classical, to our knowledge, spatiotemporal profiles of the ionic concentrations, potential, electric field, and both the faradaic and displacement currents have not yet been shown for a model formulation based on the full set of the PNP equations.

Equation 4 represents the electrodiffusion equation, also known as the Nernst−Planck equation. In general, movement of charge is accompanied by a change of the electric field due to a change in the charge density. The charge density ρ (in C m−3) is related to the electric field J by the Poisson equation: ∇·E =

∂ci = −∇·Ji ∂t

E = −∇φ

∂ci ⎡ ⎤ F = Di⎢∇2 ci + (∇ci·∇φ + ci∇2 φ)⎥ ⎣ ⎦ ∂t RT

(8)

In the absence of extra charges the charge density in an electrolyte solution is given by the concentrations of all ions in the solution. Then, with respect to eq 7, the Poisson equation may be rewritten as: ∇2 φ = −

F ε

∑ zici

(9)

i

For the nonstationary problems, time change of the electric field induces the displacement current, whose contribution to the total electric current may be significant in electrolyte solutions.22,26,33 Thus the total current density iT in a transient state of an electrodiffusion system should be considered as the sum of the faradaic (conducting) current iF and the displacement current iD in accordance with the Maxwell’s definition of “the true electric current”:34

(2)

i T = iF + iD = iF + ε

∂E ∂t

(10)

It should be noted that the term “faradaic current” used in the following text also includes the charging current in contrast with some electrochemical literature that treats the charging current as nonfaradaic. An interesting feature of the total electrodiffusion current is its spatial uniformity. This feature emerges as a consequence of the Poisson and continuity equations that can be used to show that ∇·iT = 0. Thus the faradaic and displacement currents obey the relation

(3)

where Di and ci are the diffusion coefficient and the molar concentration of ith ion, respectively, R is the universal gas constant, and T stands for the absolute temperature. Taking into account the fact that the chemical potential depends on concentration, the flux can be expressed as ⎛ ⎞ zF Ji = −Di⎜∇ci + i ci∇φ⎟ ⎝ ⎠ RT

(7)

and applying the continuity equation to the Nernst−Planck equation (eq 4) for the ion flux, the time derivative of the ion concentration is expressed as:

where the driving force is represented by a negative gradient of the corresponding potential. In the case of electrodiffusion, by using the Einstein−Smoluchowski relation between the mobility and the diffusion coefficient for the flux of ith ion Ji, we obtain:

Di ci( −∇μ ̃) RT

(6)

where t is time. Considering the relationship between the electric field and the electric potential

where μ is the chemical potential, z is the charge number, F stands for the Faraday constant, and φ is the electric potential. Usually, a mass transport is quantified by the flux J, defined as the amount of the particles (in moles) transported per unit time through unit area. The flux is related to its driving force by the Teorell formula:32

Ji =

(5)

where ε is the permittivity of the medium. A system of the Nernst−Planck equations describing the fluxes of individual particles (ions) and the Poisson field equation constitutes a mathematical framework for a deterministic description of the electrodiffusion in the continuum approximation. For the time-dependent systems, the mass conservation law takes the form of the continuity equation that relates the fluxes and concentrations:

2. THEORY The primary cause of ion transport is spatial nonuniformity of the electrochemical potential μ̃ μ ̃ = μ + zFφ (1)

flux = mobility × concentration × driving force

ρ ε

∇·iF = −∇·iD (4)

(11)

The proof for the 1D case can be found in ref 35. 14284

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

⎡ ⎤ ∂c 2 F = D2⎢∇2 c 2 − (∇c 2∇φ + c 2∇2 φ)⎥ ⎣ ⎦ ∂t RT

3. MODEL The membrane channel is represented by a rectangular domain of the length L and the height arbitrarily set to 3 nm. Although the model geometry is 2D, the problem to solve is intentionally defined as pseudo-1D, as explained below. This simplification makes it possible to compare solutions obtained by a PDE solver in 2D with the existing 1D solutions and to test reliability and applicability of the chosen numerical technique to deal with the electrodiffusion problems of higher spatial dimensionality. Thus the obtained solutions are functions of time and x coordinate only and the height of the domain has no effect on the presented results. The interior of the domain is filled with a water solution of a uniunivalent strong electrolyte without any internal or boundary fixed charges. The left (x = 0) and right (x = L) boundaries were kept at constant electrolyte concentrations c0 and cL, respectively. Initial concentrations of cation and anion were mutually equal with linear distribution along the x coordinate, conforming to the constant initial gradient across the membrane imposed by the boundary conditions for ionic concentrations. The left boundary was grounded at zero potential the entire time, while the right one was subject to a sudden (at t = 0) increase in the potential from zero to a constant value V. The lower and upper boundaries of the domain are subject to the Neumann boundary conditions, prescribing vanishing gradients of the ionic concentrations and the potential: ∂c/∂y = 0, ∂φ/∂y = 0. The initial and side boundary conditions are formulated without any gradients in the y direction, so the solution does not depend on y coordinate. Such specification of the initial and boundary conditions ensures effective one-dimensionality. A scheme of the model is shown in Figure 1.

∇2 φ = −

(14)

where the subscripts 1 and 2 are related to cation and anion, respectively, ε0 denotes the vacuum permittivity, εr is the relative permittivity of the medium, and other parameters have the meanings as explained in the previous section. The initial boundary value problem, eqs 12−14, was solved with the initial and boundary conditions as previously described using the following standard set of the parameters: D1 = 1 × 10−6,D2 = 5 × 10−7 cm2 s−1, T = 293 K, and εr = 80. The parameter values were taken from ref 16 with minor modifications. The channel length was extended to L = 10 nm, which is a typical thickness of the membrane. Standard boundary concentrations for the electrolyte were c0 = 50 mM and cL = 500 mM, and size of the potential step was normally set to V = +100 mV. In the following, the previous conditions are referred to as standard ones, and the resulting solution is referenced as the standard solution. When the effect of a parameter change on the model behavior was investigated, only the particular parameter value was varied, while the other ones were held at their standard values (unless specifically noted).

4. METHODS The equations 12−14 were solved with the use of an adaptivegrid finite-difference solver for time-dependent 2D systems of partial differential equations VLUGR2.36 The solver is coded in Fortran 77, and it implements a local uniform grid refinement method for the adaptive space discretization. The resulting system of ordinary differential equations or differential-algebraic equations obtained after the space discretization is integrated in time using the second-order two-step implicit backward differentiation formula with a variable step size. More detailed information about algorithmic aspects and the code structure of the VLUGR2 can be found in ref 36. We used a double-precision version of VLUGR2. For a standard calculation, the domain was covered by a basic grid of 10 × 3 cells in the x and y directions, respectively, with a grid refinement allowed up to level 4. The initial time step size DT = 10−11 s was used, and the space and time tolerances were usually specified as TOLS = 0.005 and TOLT = 0.003. The default option BiCGStab with ILU preconditioning was used as a linear solver. The Jacobian matrix was computed by numerical differencing. Convergence was tested by comparing the results obtained with various grid densities and time-step sizes. Calculated spatiotemporal profiles of the ion concentrations and the potential can be used to evaluate other quantities of interest by an a posteriori numerical differentiation of the obtained solution. The ion fluxes and their decomposition to diffusion and migration terms were determined according to eq 4. The ion fluxes were subsequently used to determine the faradaic current as the net charge flow through the channel. The electric-field intensity was evaluated using its relation to the electric potential (eq 7). The displacement current was determined from eq 10 using the second mixed spatiotemporal derivative of the potential. Numerical evaluations of the first and diagonal second derivatives (∂2/∂x2) were obtained using a seven-point difference (generally nonequidistant) quadraturelike scheme that suppresses the contamination from higher (up to sixth) order derivatives. The accuracy of the calculated mixed

Figure 1. Scheme of the used model of a simple membrane channel.

Time evolution of the ion concentrations and the electric potential after the potential jump was calculated by numerically solving the governing PDEs. The results are presented in the form of quantities Δc, Δφ representing departures of the calculated concentrations and electric potential from the initial values given by the initial linear distribution in the domain interior. Below we will refer to Δc and Δφ as difference concentration and potential, respectively. This model is inspired by a model of the gramicidine A channel16 with some modification and simplification. The governing equations (in the scalar form) are arising from the PNP theory as outlined above: ⎡ ⎤ ∂c1 F = D1⎢∇2 c1 + (∇c1∇φ + c1∇2 φ)⎥ ⎣ ⎦ ∂t RT

F (c1 − c 2) ε0εr

(13)

(12) 14285

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

spatiotemporal second derivatives (∂2/∂x∂t) was verified by comparing the ∂(∂/∂t)/∂x and ∂(∂/∂x)/∂t values. For the purposes of the a posteriori data processing, the solution profiles were calculated on a fixed uniform grid of 160 × 3 cells without further refinement, and the maximal time step was limited to 10−8 s. Accuracy of the obtained numerical solutions was tested by comparing the calculated steady-state profiles with the limit solutions available in the form of analytical expressions.37 The results confirmed that the numerical solution for low electrolyte concentrations approaches the “short-channel limit” corresponding to the constant field assumption with the decreasing channel length. At higher ionic concentrations and longer channels, the electroneutrality condition becomes valid and the solutions were close to the “long-channel limit”. The steadystate concentration profiles calculated by solving the PNP system with the right-hand side of the Poisson eq 14 set to zero fully reproduced the short-channel limits for both ions known as the Goldman−Hodgkin−Katz equations. The existence of an asymptotic analytical solution of the time-dependent PNP equations found by MacGillivray24 for early times and small voltage-clamp values allowed for testing the used numerical approach also with respect to temporal dynamics. The calculated time evolution of the ionic perturbation concentrations matched the theoretical asymptotes for times satisfying a validity condition for the asymptotic solution. Verification of the numerical solution reliability was also performed by comparing the left- and right-hand sides of eqs 12−14, evaluated by means of the a posteriori numerical differentiation of the obtained solution, as previously described. As an additional test, verification of spatial uniformity along the x coordinate of the total current (see Section 2) was carried out. More details of the testing of the adopted numerical technique for its ability to deal with the time-dependent PNP equations together with some illustrative results can be found in ref 38. Data processing and graphics were performed within the Gnuplot 4.4 environment. Rate constants were evaluated by nonlinear fitting with exponential functions using the Gnuplot implementation of the nonlinear least-squares Marquardt− Levenberg algorithm. Three-dimensional color plots were produced using Matlab 7.11.0.584.

Figure 2. Difference ionic concentrations (scale on left vertical axes; curves marked with ‘+’ for cation and ‘×’ for anion) and the electric potential (scale on right vertical axes; curves marked with circles) for x = 2.5 nm (panel a) and x = 7.5 nm (panel b) as functions of time after the application of a constant voltage across the membrane. The standard conditions (see Section 3) were used in the calculations. The initial and final phases of the process were fitted by exponential functions y = a(1 − e−kx) (lines 1 and 3) and y = ae−kx + c (lines 2, 4, 5). The rate constants ki (corresponding to the lines i = 1−5) obtained from the fits are as follows (in μs−1): k1 = 127, k2 = 8.36, k3 = 65.9, k4 = 9.23, and k5 = 156 (for x = 2.5 nm, panel a) and k1 = 493, k2 = 5.95, k3 = 454, k4 = 5.91, and k5 = 269 (for x = 7.5 nm, panel b).

this phase, the electric potential rapidly rises almost to its final steady-state value, closely following the exponential course also during the next phase, which may be seen as an interchange period. This phase is characterized by simultaneous occurrence of both processes demonstrated as a maximum in the cation concentration. After ∼40 ns from the voltage jump, the second, considerably slower, exponential process becomes dominant. During this period, the ionic concentrations decrease in a relaxation process to their steady-state values depending on the particular position in the channel. In contrast, the potential varies only a little in this phase, being almost fully developed by the first process. The described behavior can be explained in the framework of Planck’s theory, followed and complemented later by Cole23 for biological membranes. The first very fast process in the socalled Planck sequence may be treated as a simple “charging process” corresponding to the dielectric relaxation of Maxwell with a time constant τC = RC for the resistance R and capacity C of an element (the channel in our case). For an electrolyte solution with a nonuniform concentration, the time constant can be expressed as

5. RESULTS AND DISCUSSION 5.1. Mechanism of the Electrodiffusion Transport. Kinetic (time-dependent) PNP models using the current-step condition22,23,25,39 are more frequently used than models applying the voltage-step condition.24,40,35,41 This is somewhat surprising because the latter are more relevant to physiology, particularly regarding the voltage-clamp42 or the patch-clamp43 experimental techniques. Already, Planck12,44 anticipated that the transient behavior of an electrolyte solution after the application of a constant current may be described by two main processes with different time constants. Accordingly, kinetic behavior of our model system after an initial jump in the membrane potential also indicates two principal processes. Figure 2 demonstrates that time evolution of the ion concentration deviations Δc from the initial values (at onequarter of the channel length) can be represented by a combination of two exponential processes. The first one begins immediately after the initial excitation and lasts for ∼15 ns. During this phase, the cation concentration increases, while the concentration of anion decreases as the system departs from the local electroneutrality in a process of charge separation. In

τC =

ε LΛ

∫0

L

dx c(x)

(15)

where Λ and c(x) are the molar conductivity and concentration of the electrolyte, respectively. Approximating the concentration profile in the channel by the initial linear ionic distribution and considering the relation of the molar conductivity to an average diffusion coefficient of the ions D, after integration, we have c εRT ln L τC = 2 c0 2F D(c L − c0) (16) 14286

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

Figure 3. Spatiotemporal profiles of the ionic difference concentrations ((a) for cation (also see 3D Rotatable Image 1), (b) for anion (also see 3D Rotatable Image 2)) and the difference electric potential (c) (also see 3D Rotatable Image 3). Parameters of the calculation are the same as those in Figure 2 − the standard conditions. 3D rotatable images of Figure 3a−c in vrml format are available.

simulation in Figure 2a for cationic and anionic redistribution, respectively (lines 2 and 4). It should be noted that the rate (or time) constants of both of the principal processes moderately depend on the position in the channel (x coordinate), as it will be discussed later. The demonstrated data illustrate that the slower redistribution process is also very fast, being only less than two orders of magnitude slower than the initial dielectric relaxation. Thus, under the conditions studied here, the principal processes are not separated by a pseudo-steady state as Planck suggested for the liquid junction problem, but the phases significantly overlap in the interchange period. The duration of this interchange period is longer in the channel locations closer to the right boundary, as it can be seen in Figure 2b and will be discussed also in later subsections. Spatiotemporal profiles of the ion difference concentrations and the electric potential are shown in Figure 3. Useful information on the origin and mechanism of the analyzed electrodiffusion transport can be obtained by comparing magnitudes of the terms in eqs 12 and 13, which contribute to the total rates of time variations of the ionic concentrations. Such a term separation in eqs 12 and 13 calculated by the numerical differentiation of the solutions for cation and anion in Figure 2a is illustrated in Figure 4. For the cation (Figure 4a), the second term is almost constant, all the time representing a positive background initially set by the boundary conditions at D1F(cL − c0)V/ (RTL2) = 1.78 M μs−1. This term in eqs 12 and 13 also accounts for the initial rates of the ionic concentration changes because the Laplacians in the first and third terms are negligible

where the symbols have the meanings as previously explained. Using the standard set of parameters yields a value τC = 6.32 ns that is in surprisingly good agreement with the value of time constant 1/k5 = 6.41 ns obtained by fitting the numerical solution for the potential at x = 2.5 nm (Figure 2a, line 5). However, at the position x = 7.5 nm, closer to the right boundary with the applied potential, the value of 1/k5 = 3.71 ns (Figure 2b, line 5) is somewhat lower than the estimate. For the current-step conditions, the potential difference across the membrane is at this stage always proportional to the applied current, so the system behavior in the first period is a linear process. Also, for the voltage-step case, it can be shown that local potentials and the total current develop proportionally in this stage. The mechanism of the ion movement in the second exponential phase is, however, more complex, reflecting mutual interaction of the changing electric field and the ionic fluxes. This period is characterized by redistribution of ions after the charging process toward a new steady state, which was termed “the slower principal nonlinear process” by Planck. As Cole noted in his book45 “the process is quite complicated and neither Planck nor anyone else, until long after, gave an expression for its speed”. Cole23 estimated the time constant of the redistribution process as τD =

L2 π 2D

(17) −7

2 −1

For L = 10 nm and D = (D1+D2)/2 = 7.5 × 10 cm s , eq 17 yields τD = 135 ns. This value compares very well with the time constants 1/k2 = 120 ns and 1/k4 = 108 ns from the 14287

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

Effect of the Channel Length. Making the channel domain shorter causes amplification of the internal electric field induced by the applied voltage, which results in an increased effect on the difference ionic concentrations and electric potential in comparison with the standard conditions (Figure 5).The

Figure 4. Time dependence of the terms in eqs 12 (a) and 13 (b) obtained by numerical differentiation of the solution shown in Figure 2a.

at the beginning of the process due to the initial constant gradient of the ion concentrations and potential. The initial rates determined from the exponential fits in Figure 2a to be 2.02 and −0.907 M μs−1 for cation and anion, respectively, agree with the above value for cation and a value of −0.891 M μs−1 calculated analogously for anion as the initial value of the second term in eq 13 (Figure 4b). The electrodiffusion kinetics during the charging process (0−15 ns) is governed by the third term of eqs 12 and 13, as the electric potential is the most dynamically developing quantity in this stage of the process. Because of the large positive background of the second term of eq 12, the concentration of cations increases, albeit with decreasing rate. After 15 ns, the first (diffusion) term of eqs 12 and 13, which was negligible during the first period, comes into play, causing compensation of the negative electromigration rate represented by the sum of the second and third terms of the corresponding governing equations. While the electromigration rate of the anionic transport is increasing and stays negative during the whole process, the initially positive cationic electromigration rate becomes positive and drops below the rate of diffusion transport. This results in a maximum on the time course of the cationic concentration (Figure 2). After ∼40 ns, from the initiation of the whole process, the electromigration rate stabilizes at an almost constant level, and the transport dynamics is determined mainly by the diffusion term. Finally, when the negative rate of electromigration is fully compensated by diffusion, the steady state is reached. Figure 4 shows that the electrodiffusion temporal dynamics arise from variable balance of several processes contributing to the net transport by relatively high levels of magnitude. This appears, together with significantly different time scales of the two principal processes, to be the main source of difficulties encountered when attempting to solve the time-dependent PNP equations. 5.2. Effect of the Parameters on Behavior of the System. We investigated effects of changes in model parameters such as channel length, external potential, boundary ion concentrations, and diffusion coefficients on the kinetics of ion transport in our model system. In some cases, the calculated steady states are also discussed. The most important results are presented in the following paragraphs; because of the large amount of obtained numerical data, some of the figures are available as Supporting Information.

Figure 5. Time courses of the difference ion concentrations (a) and the difference potential (b) at x = 2.5 nm after the voltage step (100 mV) and the corresponding final steady-state profiles (c,d) for various lengths of the channel L shown as labels (in nanometers). In panels a and c, the upper lines correspond to cations and lower lines to anions. The line labeled ‘∞’ in panel d represents the steady-state profile for the electroneutrality condition. The other model parameters are set to the standard values.

maximum in the difference concentration of the cation at one-quarter of the channel length L = 5 nm is increased to almost 40 mM. The time constant of the first principal process varies very little with the channel length in accordance with eq 16 for the time constant of dielectric relaxation. For illustration, the time constants of the potential rise in Figure 5b were determined to be 4.66, 5.70, and 6.41 ns for L = 5, 7.5, and 10 nm, respectively. The effect of the channel length on the rate of the second principal process is more pronounced. The time constant of the cationic concentration drop in Figure 5a (35.9, 75.8, and 120 ns for L = 5, 7.5, and 10 nm, respectively) increases proportionally to L2, as predicted for kinetics of the redistribution process (eq 17). Such proportionality holds for the redistribution of anions as well, and it is preserved throughout the channel. Different dependence of the rates of the principal processes on the channel length suggests that for shorter channels the redistribution becomes more coupled to 14288

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

the initial charging process. Thus, for biological membranes, the principal processes are not well separated, and complex electrodiffusion interaction is to be expected. The channel length affects significantly also the steady state established after the membrane excitation. Figure 5a indicates that the steady-state concentration levels of cation and anion diverge from each other as the channel length decreases, indicating deviation from the initial electroneutrality. Spatial steady-state profiles of the ion concentration differences for the three selected channel lengths are shown in Figure 5c. It is obvious that validity of the electroneutrality condition even for the steady state is questionable because a local decrease in the anionic concentration reaches a value of almost 60 mM in some regions of the 5 nm channel. The constant field approximation fails even more significantly under these conditions, predicting the steady-state concentration differences as high as 187 mM. The deviation of the final steady state from the constant field assumption is directly demonstrated by the difference potential Δφ plotted in Figure 5d. For the short channel of 5 nm, there is still considerable difference of more than 20 mV from the value predicted by the linear potential profile. However, for longer channels of L > 10 nm, the steady-state profile of electric potential is quite satisfactorily described by the long-channel limit (the ‘∞’-line) arising from the electroneutrality condition. Effect of the External Potential. Increase in the potential step from +100 to +200 mV brings about similar qualitative effects on time course of the ion concentrations as the shortening of the channel. The maximum of the cationic concentration at x = 2.5 nm reaches 23 mV. (See Figure S1 in Supporting Information.) Variation of the boundary potential also affects the steady-state distribution of both cations and anions by shifting the ionic concentrations to lower total values. Change of polarity of the applied voltage across the membrane results in a reversal of kinetic behavior of the ions. The maximum can be observed in the time course of the anion concentration, while the cation concentration follows a sigmoid curve as shown in Figure 6a. The sigmoid course is deformed by subtle local extrema in the interchange period at channel locations closer to the right boundary (Figure 6b). Evolution of the electric potential shows a sigmoid decrease with a small minimum for higher negative potentials; see Figure 6c. The steady-state difference concentrations of ions for negative and positive voltages are equal, but the concentration profiles of cations belong to the corresponding profiles of anions and vice versa. An increase in the transmembrane voltage leads to depletion of the electrolyte from the channel manifested in Figure 6d by a shift of the final steady-state concentrations of both ions to lower values in most of the channel interior. Variation of the external potential in the range from −200 to +200 mV exhibits no significant effect on the time constants of the two principal processes. This finding is in accordance with the theoretical assumption that rate of the charging process is determined mainly by the electric resistance of the electrolyte solution in the channel. Also, kinetics of the redistribution process is governed primarily by diffusion; hence the assumed rate-determining factors are the diffusion coefficient and the length scale rather than the electric field. Effect of the Relative Permittivity. Choosing a suitable value for the permittivity of the medium in the channel interior is a matter of ongoing discussions. The PNP theory is based on a continuum electrostatic approximation using the mean charge density of the diffusing ions. Because of its “mean-field” characteristic,46 the electric forces are described on a macro-

Figure 6. Time courses of the difference ion concentrations ((a) for x = 2.5 nm, (b) for x = 7.5 nm) and the difference potential ((c) for x = 2.5 nm) after various negative voltage steps shown as labels in panel a (in millivolts). Panel d shows the corresponding final steady-state profiles of the difference ion concentrations. In panels a, b, and d, in each pair of equally colored lines the upper line corresponds to anions and the lower one corresponds to cations.

scopic level with the use of a single permittivity value averaging the dielectric properties of the medium. Such a mean-field approach certainly has its limitations because the spatial scale of real membrane channels suggests that the microscopic structure of the channel may be of significance. In this respect, the amount and structure of water molecules in the channel interior is uncertain, and hence the use of the “macroscopic” value for the relative permittivity of water is questionable. Therefore, we also examined the effect of several lower values of εr on the behavior of the studied system. A decrease in the relative permittivity strengthens Coulombic interactions of ions and consequently hinders the process of charge separation. This results in attenuation of the transient effects. For example, the maximum in the difference concentration of cations at x = 2.5 nm decreases from 13.5 to 7.51 mM if the standard value of εr is reduced twice. (See Figure S2 in the Supporting Information.) The rate of the first principal process is directly affected by the permittivity of the medium, as expected for a process of dielectric relaxation. Calculated time constants are proportional to the relative permittivity values, in agreement with eq 16. The rate of the redistribution process was found to be independent of the medium permittivity. Finally, the steady state is influenced by the permittivity of the medium. With decreasing εr, the steady-state difference ion concentrations 14289

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

approach zero as the electroneutrality condition becomes effective. Hence the steady-state potential profile for lower permittivities (εr < 20) is very close to the logarithmic function representing the long-channel limit (Figure S3 in the Supporting Information). Effect of the Boundary Concentrations. A change of the boundary concentrations toward reduction of the electrolyte concentration difference across the membrane naturally results in suppressing the transient fluctuations of the ionic concentrations and decreases the steady-state charge separation. When the electrolyte concentrations at both channel boundaries were set equal, no effect of the potential step could be detected in our numerical simulations. Reversing the concentration gradient to be opposite to the gradient of the potential causes the most profound transient effects to be shifted to positions closer to the right boundary. Under such conditions, the transient maximum is reduced, and it is exhibited by the anionic instead of the cationic concentration time course. The effect of the boundary concentration on the rates of the principal processes was examined using conditions with the left boundary concentration fixed to the standard value c0 = 50 mM, and the right boundary concentration cL varied from 100 to 500 mM. The rate constants determined for the charging process linearly depend on the concentration cL, and the time constants obtained from fitting the difference potential time courses follow the concentration dependence predicted by eq 16. The right boundary concentration has only a very weak effect on the redistribution rate constants. Effect of the Diffusion Coefficients. The effect of ionic mobilities represented by the diffusion coefficients on the studied system was investigated most comprehensively using a wide range covering two orders of magnitude of the variable parameters (D1 or D2). In addition to the standard set of fixed parameters, conditions with reversed boundary concentrations were also analyzed. The effect of increasing the diffusion coefficient of anion on the transient kinetics is shown in Figure 7. It is obvious that the lines representing time course of the anionic concentration diverge from each other very soon after the initiation of the process, in contrast with the corresponding lines for cations. This can be simply explained by different initial rates that are proportional to the diffusion coefficients of ions (varied for anion and fixed for cation). However, the rate constants of the first exponential (charging) process are very similar for both ions and show linear dependence on the diffusion coefficient of anion. Increasing the anionic mobility above the value set for cation brings about decreasing the maximum of the cationic concentration and occurrence of small local extrema in the interchange period of the anion time course. The time courses of the difference concentration for the particular ion approach the same steady-state value, demonstrating general independence of the steady state on diffusion coefficient. The time courses of the difference potential (Figure 7d) shift to the left due to the acceleration of the charging process, and a subtle maximum gradually appears. The rate constants of the second exponential (redistribution) process are again close for cations and anions. The reciprocal rate constants (i.e., time constants) depend linearly on the reciprocal value of the diffusion coefficient. The rate constants of both principal processes obtained for the center of the channel by fitting the calculated data are collected in Table 1. Similar linear dependences were found also for other channel locations, although the slopes increase with the increasing spatial coordinate.

Figure 7. Time courses of the difference ion concentrations ((a) for x = 2.5 nm, (b) for x = 5 nm, and (c) for x = 7.5 nm) and the difference potential ((d) for x = 5 nm) for various ratios of the diffusion coefficients of anion and cation D2/D1 as labeled in panel a. The diffusion coefficient of cation was fixed at the standard value. In panels a−c, in each pair of equally colored lines, the upper line corresponds to cations and the lower one corresponds to anions.

The effect of decreasing the anionic mobility for the standard and reversed boundary electrolyte concentrations can be found in the Supporting Information (Figures S4 and S5). Figure S6 in the Supporting Information illustrates the effect of increasing the cationic diffusion coefficient on the time courses of the difference ionic concentrations and electric potential. 5.3. Electric Field, Ionic Fluxes, and Electric Currents in the Channel. Time evolution of the electric field for various positions in the channel is shown in Figure 8. The vector of E is always negatively directed from the right to the left. For clarity, only the spatial derivative of the electric potential is shown. As with the potential, the most dramatic changes of the field intensity occur during the charging process. The initially uniform field of 10 mV nm−1 increases with time in the left half of the channel and simultaneously decreases in the right part of the channel. After 50 ns, the field profile is very close to the steady state with a value of ∼23 mV nm−1 near the left boundary and 4 mV nm−1 close to the right one. The steadystate profile predicted assuming the electroneutrality condition is marked in Figure 8b. Clearly, for x > 5 nm, the electroneutrality steady-state profile corresponds well with the full PNP steady-state solution; however, it differs significantly in the left region of the channel. 14290

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

Table 1. Rate Constants of the Two Principal Processes for Various Diffusion Coefficients of Anion D2 at x = 5 nm kC [μs−1], charging process −1

D2 [nm μs ] 2

10 25 50 100 200 400 1000

(+) 213 248 303 400 569 951 2061

a

± ± ± ± ± ± ±

(−) 2.6 4.0 5.9 8.9 10 19 39

175 197 234 311 463 763 1662

± ± ± ± ± ± ±

kD [μs−1], redist. process (φ)

0.66 0.92 0.88 0.79 0.64 0.78 1.6

168 183 210 262 359 608 1341

± ± ± ± ± ± ±

(+) 0.12 0.41 0.89 1.9 2.7 3.8 6.2

1.61 3.70 6.41 10.1 14.7 18.7 21.9

± ± ± ± ± ± ±

0.015 0.030 0.052 0.064 0.078 0.31 0.48

(−) 1.61 3.71 6.46 10.2 14.8 18.4 21.0

± ± ± ± ± ± ±

0.017 0.032 0.050 0.065 0.071 0.25 0.30

Data in the columns marked (+), (−), and (φ) were obtained by fitting the simulated time courses Δc1(t), Δc2(t), and Δφ(t), respectively. The fixed model parameters were set to their standard values in the simulations.

a

Figure 9. Time courses of the ionic fluxes ((a) cation and (b) anion) through the channel obtained from the standard solution. The lines correspond to various positions in the channel indicated as labels (in nanometers).

Figure 8. Time evolution of the electric field in the channel interior. The lines in panel a represent time courses at positions indicated as labels (in nanometers) for the outer lines. The inner lines correspond to x = 1, 2, ..., 8 nm (from the top to the bottom). The panel b shows spatial profiles at times marked by labels in nanoseconds. The labels ‘ss’ and ‘en’ stand for the full PNP steady state and the electroneutrality limit (red line), respectively. The displayed data were obtained by the differentiation (see Section 4) of the standard solution.

be taken to account. The displacement current is defined as the rate of change of the electric field. (See eq 10.) As we have previously shown, the electric field in the studied system does not display dramatic changes in space, but the changes occur within an extremely short time period. Hence, considerable magnitudes of the displacement current might be expected. Both current components as well as the total electric current are presented in Figure 10. As a consequence of the negative ionic net flux, the faradaic current is always negative, extending from approximately −0.5 pA nm−2 near the left channel boundary to −3 pA nm−2 at the right end. The displacement current also shows initial spatial nonuniformity varying with the position from −1.2 to +1.1 pA nm−2. Thus, the displacement current represents a significant part of the total current in the system studied here. For example, at 1 ns from the voltage jump, the displacement current accounts for >60% of the total current magnitude at the position x = 0.25 nm. Despite the spatial inhomogeneity of the faradaic and displacement currents, the total current is uniform throughout the channel (Figure 10b), as predicted by the theory. The total current density at t = 0.3 ns from the process initiation is −1.77 pA nm−2 and increases with time to the steady-state value of −1.35 pA nm−2 under the standard conditions. The corresponding change of the total current for the voltage step of +200 mV ranges from −3.33 to −2.27 pA nm−2 and from −0.994 to −0.801 pA nm−2 for the voltage step of +50 mV. These values illustrate surprisingly low time dependence of the total current

Another piece of valuable information on spatiotemporal properties of the ionic fluxes available from the a posteriori processing of the obtained solution is the time course of the cationic and anionic fluxes for the selected positions in the channel interior presented in Figure 9. Because the mobility of the cation is twice that of the anion for standard conditions, the cationic flux is dominant in contributing to the net flux that accounts for the conductive (faradaic) current through the channel. The cationic flux is always negative, initially ranging between −7 and −23 M nm μs−1, and significantly depends on the channel position. However, the spatial nonuniformity levels out after ∼20 ns when the charging process overlaps with the redistribution process. (See also Figure 2.) The anionic flux is mostly positive, except for early times in the close proximity to the left boundary. The initial values cover a range from −1 to 7 M nm μs−1 and preserve the spatial nonuniformity somewhat longer than in the case of the cationic flux. The net flux varies from −6 to −30 M nm μs−1 with the increasing x value. The transport of ions due to the ionic fluxes gives rise to charge conductance observable as the electric current. However, in dynamical systems with an electric field changing over time, another component of the total electric current, the so-called displacement current introduced by Maxwell,34 should 14291

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

scale of nanometers, where individual particles rather than average concentrations and charge density should be considered. To account for the statistical nature of transport on the subcontinuum scale, alternative computational approaches have been developed. For instance, the Lattice Boltzmann method,47 frequently used in fluid dynamics, simulates the interactions of mesoscopic particle populations using discrete velocities and positions. Computational demands of such methods are, however, considerably higher than those for solving continuum equations. Nevertheless, mean-field theories applicable to microscale problems are often useful, although with limitations, also for the nanoscale level.48 For example, the basic or extended PNP approach was successfully used for modeling steady-state properties of real ion channels.16,49,50 This suggests a cautious optimism that the development of the time-dependent PNP models can be helpful in understanding ion-channel characteristics obtained experimentally under dynamical conditions of an alternating external electric field or of a sequence of various potential steps. The studied model describes only transport kinetics of free ions in a single channel using deterministic and macroscopic approach. However, nonlinear electric properties of excitable cells are set by collective behavior of an ensemble of channels,51 each of them subject to individual stochastic gating.52,53 Thus, the whole cell current dynamics is determined not solely by the electrodiffusion transport of ions but also by the transition kinetics between discrete states of individual channels.

Figure 10. Time courses of the current density through the channel obtained from the standard solution. (a) Faradaic (the lower set of lines − red) and the displacement (the upper set − blue) currents at various positions in the channel. Outer lines of both line sets correspond to the positions marked as labels in nanometers, and the inner lines correspond to x = 1, 2, ..., 8 nm. (b) Total current representing the sum of the corresponding time courses from panel a.

due to mutual compensation of the faradaic and displacement currents. 5.4. Limitations of the Model. The presented model is one of the most simplified representations of a membrane channel. Hence, the obtained results should be interpreted and discussed mainly from the theoretical point of view with only tentative implications to real biological systems. Despite its simplicity, our model upgrades most of the existing basic timedependent PNP models by solving the full Poisson equation and considering contributions of both ions to the dynamically changing charge density. Furthermore, the model can be readily extended to multi-ionic systems with more structured domain geometry in 2D or 3D. Biological ion channels are complex objects with complicated internal structure, often forming nanodomains such as pores and selective filters. The usual property of membrane channels is the presence of surface charges on the channel walls. Our simple model does not consider any charge and hence a potential at the lateral boundaries. Moreover, the model ignores channel entrance and exit, which can be modeled as reservoirs of solution at the side boundaries in models with a higher spatial dimension than one. Because of these simplifications, electrodiffusion dynamics reported here may significantly differ from those of real systems. However, the presented simulations help us to understand fundamental mechanisms and effects of general validity extending even to more realistic systems. Our motivation is the development of a computational tool potentially capable of solving complex electrodiffusion systems such as biological ion channels or any other relevant real system, for example, fabricated nanofluidic devices. In the next step, the tested methodology can be applied to advanced models, and additional features such as boundary charges or the presence of a channel entrance/exit can be analyzed. In this respect, our contribution may be considered as an initial step providing a basic reference in a step-by-step study of effects of increasing complexity of the model toward real systems. As we have already mentioned, the continuum approximation used in the PNP theory implies another limitation because the processes in ion channels of biological membranes occur on a

6. CONCLUSIONS We employed numerical techniques for solving time-dependent PNP equations in 2D to analyze transient behavior of a simple ion channel subject to a sudden electric potential jump across the membrane (voltage-clamp). Calculated spatiotemporal profiles of the ionic concentrations and electric potential showed that two principal exponential processes can be distinguished in the electrodiffusion kinetics, in agreement with the Planck’s theory. The initial fast process corresponds to the Maxwell dielectric relaxation with a typical time constant of 5 ns for the standard conditions. The steady state is approached in a slower exponential process attributed to the nonlinear redistribution process with a standard time constant of 160 ns. Depending on various conditions and channel location, the principal processes are usually combined, forming an interchange period often accompanied by local extrema in the ionic concentrations and potential. Under the standard conditions, the steady state is established after ∼2 μs from the potential jump. Effects of the model parameters such as the channel length, height of the potential step, boundary concentrations, permittivity of the channel interior, and ionic mobilities on the electrodiffusion kinetics was studied. Quantitative effect of the parameters on the rates of the principal process was examined and corroborated by the analytical predictions. The obtained results show that neither constant field nor electroneutrality assumption is a sufficiently valid approximation under the constraints reflecting the data known for real membrane channels. Spatiotemporal profiles of the electric field, ionic fluxes, and both the conductive and displacement currents were derived from the obtained solution by means of numerical differentiation. We demonstrate that the displacement current is a significant transient component of the total electric current through the channel. 14292

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293

The Journal of Physical Chemistry B

Article

Undoubtedly, the used model is only a first step in constructing more realistic models to study time-dependent electrodiffusion phenomena that are far less understood than the steady-state processes. The simplified model allowed us to compare its predictions with the existing results and theoretical assumptions necessary for testing and tuning the numerical technique. The presented results provide data on the classical voltage-clamp problem solving the full PNP equations within a spatiotemporal domain and offer further physical insight into the mechanism of electrodiffusion. The used numerical approach is applicable to other systems such as synthetic nanopores and nanofluidic channels or nanopipettes.



(15) Barcilon, V.; Chen, D. P.; Eisenberg, R. S. SIAM J. Appl. Math. 1992, 52, 1405−1425. (16) Kurnikova, M. G.; Coalson, R. D.; Graf, P.; Nitzan, A. Biophys. J. 1999, 76, 642−656. (17) Ramirez, P.; Mafé, S.; Aguilella, V. M.; Alcaraz, A. Phys. Rev. E. 2003, 68, 011910. (18) Liu, Q.; Wang, Y.; Guo, W.; Ji, H.; Xue, J.; Ouyang, Q. Phys. Rev. E. 2007, 28, 051201. (19) Helfferich, F.; Plesset, M. S. J. Chem. Phys. 1958, 28, 418−424. (20) Conti, F.; Eisenman, G. Biophys. J. 1965, 5, 247−256. (21) Hafemann, D. R. J. Phys. Chem. 1965, 69, 4226−4231. (22) Cohen, H.; Cooley, J. W. Biophys. J. 1965, 5, 145−162. (23) Cole, K. S. Physiol. Rev. 1965, 45, 340−379. (24) MacGillivray, A. D. J. Chem. Phys. 1970, 52, 3126−3132. (25) Brumleve, T. R.; Buck, R. P. J. Electroanal. Chem. 1978, 90, 1− 13. (26) Murphy, W. D.; Manzanares, J. A.; Mafé, S.; Reiss, H. J. Phys. Chem. 1992, 96, 9983−9991. (27) Rudolph, M. J. Electroanal. Chem. 1994, 375, 89−99. (28) Das, S.; Dubsky, P.; Berg, A.; Eijkel, J. C. T. Phys. Rev. Lett. 2012, 108, 138101. (29) Wang, H.; Pilon, L. Electrochim. Acta 2012, 64, 130−139. (30) Grysakowski, B.; Jasielec, J. J.; Wierzba, B.; Sokalski, T.; Lewenstam, A.; Danielewski, M. J. Electroanal. Chem. 2011, 662, 143− 149. (31) Samson, E.; Marchand, J.; Robert, J.-L.; Bournazel, J.-P. Int. J. Numer. Methods Eng. 1999, 46, 2043−2060. (32) Teorell, T. Proc. Natl. Acad. Sci. U.S.A. 1935, 21, 152−161. (33) Zhou, S. A.; Uesaka, M. Int. J. Appl. Electromagn. Mech. 2009, 29, 25−36. (34) Heras, J. A. Am. J. Phys. 2011, 79, 409−416. (35) Arndt, R. A.; Roper, L. D. Math. Biosci. 1973, 16, 103−117. (36) Blom, J. G.; Tromper, R. A.; Verwer, J. G. ACM Trans. Math. Software 1996, 22, 302−328. (37) Keener, J.; Sneyd, J. Mathematical Physiology, 1st ed.; SpringerVerlag: New York, 1998; pp 82−87. (38) Valent, I.; Neogrády, P.; Schreiber, I.; Marek, M. J. Comput. Interdiscip. Sci. 2012, 3, 65−76. (39) Sokalski, T.; Lingenfelterl, P.; Lewenstam, A. J. Phys. Chem. B 2003, 107, 2443−2452. (40) Hägglund, J. V. J. Membr. Biol. 1972, 10, 153−170. (41) Offner, F. J. Theor. Biol. 1974, 45, 81−91. (42) Huxley, A. Trends Neurosci. 2002, 25, 553−558. (43) Neher, E.; Sakmann, B.; Steinbach, J. H. Pfluegers Arch. 1978, 375, 219−228. (44) Planck, M. Ann. Phys. Chem. 1890, 40, 561−576. (45) Cole, K. S. Membranes, Ions, and Impulses; University of California Press: London, United Kingdom, 1972; p 177. (46) Roux, B. J. Gen. Physiol. 1999, 114, 605−608. (47) McNamara, G. R.; Zanetti, G. Phys. Rev. Lett. 1988, 61, 2332− 2335. (48) Chang, H.-C.; Yeo, L. Y. Electrokinetically-Driven Microfluidics and Nanofluidics, 1st ed.; Cambridge University Press: New York, 2010. (49) Chen, D. P.; Xu, L.; Eisenberg, B.; Meissner, G. J. Phys. Chem. B 2003, 107, 9139−9145. (50) Gillespie, D.; Xu, L.; Wang, Y.; Meissner, G. J. Phys. Chem. B 2005, 109, 15598−15610. (51) Meunier, C.; Segev, I. Trends Neurosci. 2002 , 25, 558−563. (52) Sigworth, F.; Neher, E. Nature 1980, 287, 447−449. (53) Zahradníková, A.; Zahradník, I. Biophys. J. 1996, 71, 2996−3012.

ASSOCIATED CONTENT

S Supporting Information *

Additional figures showing time courses of the difference ionic concentrations and potential at three channel positions and a steady-state profile are available. The data describe effects of the potential jump, relative permittivity of the medium, and ionic diffusion coefficients on electrodiffusion kinetics of the studied system. This material is available free of charge via the Internet at http://pubs.acs.org. W Web-Enhanced Feature *

3D rotatable images of Figure 3a−c in vrml format are available in the HTML version of the paper.



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the European Science Foundation program FuncDyn, Slovak Research and Development Agency under the contract no. APVV-0291-11 and the grant GACR 203/09/2091 from the Czech Science Foundation. We thank Táňa Sebechlebská for assistance with the data processing and Roman Jartim for help with searching for the older literature.



REFERENCES

(1) Hodgkin, A. L.; Huxley, A. F. J. Physiol. 1947, 106, 341−367. (2) Hodgkin, A. L.; Katz, B. J. Physiol. 1949, 108, 37−77. (3) Rasmussen, H. Science 1970, 170, 404−412. (4) Berridge, M. J.; Lipp, P.; Bootman, M. D. Nat. Rev. Mol. Cell Biol. 2000, 1, 11−21. (5) Horng, T. L.; Lin, T. C.; Liu, C.; Eisenberg, B. J. Phys. Chem. B 2012, 116, 11422−11441. (6) Benzhuo, L.; Zhou, Z. C. Biophys. J. 2011, 100, 2475−2485. (7) Lim, J.; Whitcomb, J.; Bozd, J.; Varghese, J. J. Colloid Interface Sci. 2007, 305, 159−174. (8) de Paula, J. L.; da Cruz, J. A.; Kaminski Lenzi, E.; Evangelista, L. R. J. Electroanal. Chem. 2012, 682, 116−120. (9) Choi, Y. S.; Kim, S. J. Colloid Interface Sci. 2009, 333, 672−678. (10) Moya, A. A.; Hayas, A.; Horno, J. Solid State Ionics 2000, 130, 9−17. (11) Hosokawa, Y.; Yamada, K. .; Johannesson, B.; Nilsson, L. O. Mater. Struct. 2011, 44, 1577−1579. (12) Planck, M. Ann. Physik. Chem. 1889, 39, 161. (13) Goldman, D. J. Gen. Physiol. 1943, 27, 37−60. (14) Mafé, S.; Pellicer, J.; Aguilella, V. J. Phys. Chem. 1986, 90, 6045− 6050. 14293

dx.doi.org/10.1021/jp407492q | J. Phys. Chem. B 2013, 117, 14283−14293