Dynamic Theory of Membrane Potentials - The Journal of Physical

Aug 3, 2010 - New physical insights emerge. ... All simulations were programmed in C++ and run on a desktop computer with two quad core Intel Xeon 2.2...
0 downloads 0 Views 2MB Size
Download PDF
J. Phys. Chem. B 2010, 114, 10763–10773

10763

Dynamic Theory of Membrane Potentials Kristopher R. Ward, Edmund J. F. Dickinson, and Richard G. Compton* Department of Chemistry, Physical and Theoretical Chemistry Laboratory, Oxford UniVersity, South Parks Road, Oxford OX1 3QZ, United Kingdom ReceiVed: March 23, 2010; ReVised Manuscript ReceiVed: July 18, 2010

An accurate understanding of the dynamics of membrane potential formation underpins modern electrophysiology and much of cell biochemistry. Computer simulations using a Nernst-Planck-Poisson (NPP) finite difference method are used to model the dynamic evolution of a series of membrane systems in which two reservoirs of electrolyte solution are separated by a thin membrane which is impermeable to selected species. Two specific examples are considered in detail. The first (“type 1”) is the case in which the solutions are monophasic but of unequal concentration, and the second (“type 2”) is the case in which the solutions are of equal concentrations but different phase, with a common impermeant ion (a bi-ionic membrane). The validity of the Goldman equation for membrane potential, as applied to each case, is investigated. The type 1 case is shown to reach a steady state, and strong agreement with the Donnan equation for potential difference is observed. For the type 2 case, it is shown that the potential difference consists of two separable components: a localized, Donnan-type potential that reaches a pseudosteady state and a dynamically expanding diffuse component, with properties similar to those of a liquid junction potential, that does not reach a steady state but rather discharges at constant potential difference. This is contrary to the classical interpretation of a static diffuse layer, due to Planck, Henderson, and Goldman. 1. Introduction Understanding membranes underpins much of modern chemistry and biochemistry. In this paper, we develop a new approach for studying membrane potentials, generating new physical insights into the phenomenon. 1.1. Membrane Applications and Classical Understanding. A membrane is a layer of material that separates two solution phases. Often membranes are partially permeable or fully impermeable to at least one dissolved component. Such membranes find a wide range of applications in chemistry, particularly in the field of ion-selective electrodes (ISEs),1 and are common throughout biochemistry and the biological sciences.2 The most familiar example that occurs in nature is that of a cell membrane, which is vital to the existence of all complex organic life.2 Cell membranes surround all cells, separating the interior from the exterior and, of direct relevance to this investigation, control the movement of substances in and out of cells. Ion channels and pumps embedded in the cell membrane control the ionic concentrations on either side, thus generating or maintaining a potential difference across the membrane, which may both provide power and allow for the transmission of signals between different parts of the cell.2 Modern biochemistry therefore depends heavily on accurately understanding both the magnitude and dynamics of membrane potentials, that is, potential differences due to inhomogeneity between the two sides of the membrane. With solution of the Nernst-Planck-Poisson equations in full becoming increasingly computationally feasible, we intend to review this problem, to identify the extent to which the traditional chemical understanding of these systems remains valid. The steady-state potential difference of membrane systems may be calculated using the Goldman equation3 * Corresponding author. Fax: +44 (0) 1865 275410. Phone: +44 (0) 1865 275413. E-mail: [email protected].

∆φMem

RT ) ln F

(

∑ uq cq ,o + ∑ uq cq ,i +

+

q

-

-

q

∑ uq cq ,i + ∑ uq cq ,o +

+

q

+

q

-

)

(1.1)

where φ is the membrane potential, uq( is the mobility of species q(, and cq(,o and cq(,i are the extracellular and intercellular ionic concentrations, respectively. This equation is derived using several simplifying assumptions detailed below and may generally be considered valid for the case of a thick membrane. However, in the limit of a thin (compared to the Debye length) membrane, these assumptions are now largely seen as being unphysical, as recently expressed by Perram et al.4 and Dickinson et al.5 The first of these is the assumption of constant electric field, i.e., that the electric field is invariant across the membrane.3 For a potential difference to arise, there must somewhere be a nonzero charge density, which is not consistent with a constant electric field, as shown by the Poisson equation:

∂2φ F )2  ∂x s0

(1.2)

(where ∂2φ/∂x2 is the divergence of the potential, F is the charge density, and s and 0 are the permittivities of the solvent and free space, respectively). A second common assumption is that of local electroneutrality, i.e., that the charge density is zero, which, in a linear geometry, is entirely equivalent to the condition of constant electric field, and is therefore equally inconsistent with the Poisson equation when F ) 0. A further assumption in the derivation of the Goldman equation is that the external surface charge of the membrane is

10.1021/jp102599j  2010 American Chemical Society Published on Web 08/03/2010

10764

J. Phys. Chem. B, Vol. 114, No. 33, 2010

negligible.3 This leads to the use of a Dirichlet, or “constant concentration” boundary condition: the concentrations of all species at the membrane boundary are artificially held constant, as if the external reservoirs of solution were being constantly stirred. In this paper, the validity of the Goldman equation as applied to a thin membrane system and the physical meaningfulness of its underlying assumptions (particularly the condition of constant concentration) will be considered, taking advantage of modern computational approaches and recent advances in the study of membrane potentials and liquid junctions.5 1.2. Nature of This Work. This work considers the lower limit of membrane width: a membrane that is infinitesimally thin. It is assumed that the membrane is ideally permselective, i.e., that the flux of all impermeant species through the membrane is zero and that the flux of all permeant species is unaffected. Additionally, the effects of osmosis are ignored. Two specific examples are considered in detail. The first is the case in which the solutions are monophasic but of unequal concentration, i.e., systems of the type

AX(C ) C*) L

|

X

AX(C ) C*) R

where the symbol “|X” is used to indicate a membrane that is impermeable to X- and C*L and C*R are the initial bulk concentrations on the left- and right-hand sides, respectively. The second example is systems of the type

AX(C ) C*)

|

BX(C ) C*)

X

i.e., that of a bi-ionic membrane, in which the solutions are of equal concentrations but different compositions, with one common, impermeant ion. 1.3. Past Computational Studies. Since their first application to the closely related problem of the liquid junction potential by Nernst and Planck10-12 in the late nineteenth century, the Nernst-Planck (NP) equation for flux and the Nernst-PlanckPoisson (NPP) electrodiffusion equation set have been used in numerous studies, both analytical and numerical. However, it was not until the advent of digital computers that any real attempt could be made to solve the NPP set (a set of coupled nonlinear partial differential equations) without recourse to drastic simplifications of the model (e.g., constant electric field,3 electroneutrality10). Two of the earliest attempts at numerical simulation were made, both in the same year: one by Cohen and Cooley,13 who, significantly, studied not only the steady-state condition of the system but also the time dependence; the other by Hafemann,14 who was the first to avoid the approximation of constrained diffusion, i.e., the use of Dirichlet boundary conditions, by use of an expanding grid of spatial points, albeit for a free liquid junction system. Despite this, the use of finitely positioned boundary conditions, which necessarily restrict the potential difference to the area directly surrounding the membrane, has dominated,6-9 following classical analysis from Planck10 onward, though the recent work of Perram et al.4 has argued this to be inconsistent with the law of conservation of mass, and several alternatives have been proposed.5,15 Indeed, for the free liquid junction case, finite Dirichlet boundaries have been shown to be unnecessary to derive a constant potential difference.5,16 In this work, the

Ward et al. boundaries used are of Dirichlet type but positioned such that they vastly exceed the maximum extent of the diffusion layer, so the area under investigation is essentially infinite and thus unperturbed by the imposition of boundary conditions; approximations such as electroneutrality and constant electric field are not used. The majority of the work in this field has also focused only on the steady-state potential difference of the system7,15,16 with little or no regard for how this steady state is attained. Although important information may still be gained from this approach, a more complete understanding comes from examining systems in a dynamic context, which few works5,17 have attempted. This Article is intended to account for this oversight, by providing a complete, dynamic description of the systems under investigation, detailing the evolution of the potential, concentrations, and the electric field without undue approximation. New physical insights emerge. 2. Theoretical Model 2.1. The Membrane Model. We assume two different solutions of binary monovalent electrolyte to be separated by a semipermeable, infinitesimal membrane which completely blocks the transport of one or more selected species but has no effect on the mobility of the remaining species. After some time, t ) 0, the movement of permeant species through the membrane is permitted, and the evolving concentration profiles of all species, as well as their dependent properties (potential difference, electric field, etc.), may be observed. As the solutions are assumed homogeneous within all planes parallel to the membrane, mass transport is neglected in all axes except for that perpendicular to the membrane (defined as x) and is treated as one-dimensional and linear. The flux of each species q at any point in the solution is described by the Nernst-Planck equation:

Jq ) -Dq

(

∂Cq zqF ∂φ + C ∂x RT q ∂x

)

(2.1)

where Jq, Dq, and zq are the x component of the flux vector, the diffusion coefficient, and the charge, respectively, for species q, φ is the potential, and F, R, and T have their usual meanings. The first term is the diffusional contribution, and the second is the migrational one. This excepts the impermeant species at the location of the membrane, where Jq ≡ 0. By conservation of mass, the space-time evolution of Cq is given by

∂Cq ∂Jq )∂t ∂x

(2.2)

Consequently,

(

∂Cq z qF ∂ ∂2Cq ∂φ + Cq ) Dq 2 ∂t RT ∂x ∂x ∂x

(

))

(2.3)

The potential at any point must further satisfy the Poisson equation:

∂2φ F + s0 ∂x2

∑ zqCq ) 0 q

(2.4)

Dynamic Theory of Membrane Potentials

J. Phys. Chem. B, Vol. 114, No. 33, 2010 10765

2.2. Normalization. By use of the following normalization conditions, we can render the Nernst-Planck-Poisson (NPP) equation set dimensionless:

Cq C*Q,R

cq )

θ)

F φ RT

Dq′ )

Dq DQ

where concentration and diffusion coefficients are normalized against those of a standard species (Q) on the right-hand side of the junction. By defining the dimensionless time and space coordinates as

X ) kx

τ ) k2DQt

where

F2C*Q,R RTs0

k2 )

Figure 1. Evolution of the potential difference (∆θMem), maximum electric field (ξmax), and junction extent (χ98) of a type 1 membrane system with c*L /c*R ) 0.1 and DA′ /DX′ ) 1.

we obtain the dimensionless NPP equation set:

∂2cq 2

∂X

+ zq

∂θ ∂ 1 ∂cq c )0 ∂X q ∂X Dq′ ∂τ

( )

(2.5)

∂cq ∂θ + zqcq )0 ∂X ∂X

for all species q and

∂2θ + ∂X2

∑ zqcq ) 0

(2.6)

q

Using a typical diffusion coefficient of 1 × 10-9 m2 s-1 and concentration of 10 mM gives dimensionless units as approximately 1τ ) 20 ns and 1X ) 5 nm. 2.3. Outer Boundary Conditions. The outer simulation boundaries in space are taken at a distance where they will vastly exceed the diffusion layer throughout the simulation time. The maximum value of X used in the simulation, Xmax, is hence defined as

Xmax ) 6√D'τmax

(2.7)

where D′ is the greatest dimensionless diffusion coefficient in the system. The root-mean-square displacement of molecules, and therefore the mean thickness of the diffusion layer, is (2Dτ)1/2,18 so Xmax is chosen to greatly exceed it. We assume the concentrations cq(X,τ) of the ions at the boundaries are constant and equal to their bulk values, i.e., for all τ:

X f +∞

achieved in practice by setting the flux of all impermeant species to be zero on either side of the membrane, i.e.,

cq ) c*q,R

X f -∞

cq ) c*q,L

As there is no exchange of material at the boundaries, the simulation space is a Gaussian box of zero enclosed electric charge. Therefore, the electric field is zero at the boundaries, again at all τ:

X f (∞

∂θ )0 ∂X

The approximately infinitesimal membrane is modeled as existing at the center of the simulation space, which may be

(2.8)

2.4. Numerical Methods. The one-dimensional simulation space is divided into 2n + 1 points labeled X-n to Xn with X0 ) 0 as the central point. The membrane potential, ∆θMem, is extracted as θn - θ-n. The simulation uses expanding space and time grids such that there is a dense region of regular space points close to the membrane, to be denser at low X, with the same being true for time at low τ. The simulation uses a fully implicit centrally differenced finite difference discretization scheme. The set of nonlinear simultaneous equations for all chemical species and the potential at every space step may then be solved using the iterative Newton-Raphson method. Further details of the numerical methods employed, including the Newton-Raphson method, may be found in the Supporting Information. Simulation parameters that gave the most favorable compromise between accuracy and run time were determined by a detailed convergence study, which led to a runtime of approximately 4 h for a typical dynamic membrane potential difference evolution. All simulations were programmed in C++ and run on a desktop computer with two quad core Intel Xeon 2.27 GHz processors and 2.25 GB of RAM. 3. Theoretical Results and Discussion 3.1. Type 1. Here, a “type 1” membrane system is defined as a system of the type

AX(C ) C*) L

|

X

AX(C ) C*) R

Two solutions of electrolyte of the same composition (A+X-) but different concentrations are separated by an infinitesimal membrane that completely blocks the passage of X- ions. Figure 1 shows the time evolution of the membrane potential, ∆θMem, the maximum electric field, ξmax (where ξ ≡ -∂θ/∂X), and the

10766

J. Phys. Chem. B, Vol. 114, No. 33, 2010

Ward et al.

extent of the nonzero electric field, χ, for a typical system, where χ is defined as

(

χ)Xθ)

99 1 -Xθ) ∆θ ∆θ 100 Mem 100 Mem

) (

) (3.1)

It can be seen from Figure 1 that the trends in time evolution of the membrane potential (∆θMem), the maximum electric field (ξmax), and the extent of the junction (χ) all tend to a steady state. 3.1.1. Short Time BehaWior. From τ ) 0, the movement of permeant species through the membrane is permitted. At short time, the magnitude of the electric field is near zero and so the movement of the charged ions is dominated by diffusion. Diffusion of A+ down its concentration gradient proceeds until some time where the separation of charges, arising since X- is prevented from diffusing similarly, becomes significant enough for a substantial electric field to develop. This is typically on the order of 1 ns (a dimensionless time of approximately 0.1). The short time (τ f 0) behavior is not discussed by classical theory, but as a result of analysis under an approximation of zero migrational feedback (mathematical details of which may be found in the Supporting Information), it may be shown that the maximum electric field, which always occurs at the location of the membrane, is



ξmax ) (c*L - 1)

D'A · √τ π

(3.2)

since (DX′ )1/2 ≡ 1 and c*R ≡ 1, and that the membrane potential is

∆θMem ) (c*L - 1)D'Aτ

(3.3)

which was consistent with the behavior observed from simulation in all cases (figures are available in the Supporting Information). It was also observed that the initial rate of growth of χ is determined only by the rate of diffusion of the species and is independent of the initial concentrations, as expected. 3.1.2. Steady State. The accumulation of positive charge on the left-hand side of the membrane leads to a concomitant buildup of negative charge on the right-hand side (the formation of a pair of electrical double layers19), due to the action of the positive electric field on X-. The membrane “charges” as the double layers develop in time, with the electric field at the membrane increasing up until the point where the total flux of A+ at the membrane is zero. Hereafter, the diffusion into the region of lower concentration is precisely balanced by the opposite migrational flux, as envisaged classically:

(

jA ) -DA

)

∂cA ∂θ )0 + zAcA ∂X ∂X

(3.4)

In this way, the system tends toward a steady state. This is in complete contrast with the permeant case, i.e., that of a liquid junction, where the constant potential difference is not associated with a constant electric field but which may be characterized as spontaneously “charging” and “discharging” in time with the electric field relaxing to zero as τ f ∞, alongside a correspond-

Figure 2. Electric field profiles at various times for a type 1 membrane system with c*L /c*R ) 0.1 and DA′ /DX′ ) 1.

ing increase in the extent of the junction.5 This is impossible in the type 1 membrane case, as there is no possible “discharge” mechanism, since discharging would require A+ to move back against its concentration gradient, in which case a steady state would no longer be approached. Figure 2 shows the evolution of the electric field profile of the same system, clearly showing that the field tends toward a steady state at long time. Two important features of this figure are the obvious asymmetry in the electric field and the discontinuous charge separation at the membrane, as more clearly displayed in the corresponding charge density plot shown in the inset. The introduction of a semipermeable membrane necessarily introduces this discontinuity, as the charge density is the sum of continuous (A+) and discontinuous (X-) concentration profiles. Donnan20 predicted that, as the system reaches equilibrium (at long time), the membrane potential should simply be the difference in the initial electrochemical potential of the two components:

∆φMem ) φL - φR )

RT RT ln C*L ln C*R F F

(3.5)

In dimensionless terms, this is simply

()

∆θMem ) ln

c*L c*R

(3.6)

which was indeed found to be the case. Note that the Goldman equation (eq 1.1) reduces to the Donnan equation (eq 3.6) for the type 1 case. Figure 3 shows a comparison between the values of potential difference predicted by this equation and the long time results of simulations of a variety of initial concentration ratios. Strong agreement (< (0.9% difference between eq 3.6 and simulation) was seen in all cases (on the time scale studied), with a monotonic decrease in percentage error for concentration ratios closer to unity, which may be attributed to a more rapid achievement of steady state in these cases. By mathematical analysis of the system at steady state, details of which may be found in the Supporting Information, it is possible to analytically determine the maximum electric field at steady state, which was found to be in strong agreement with the results of simulation.

Dynamic Theory of Membrane Potentials

J. Phys. Chem. B, Vol. 114, No. 33, 2010 10767

Figure 3. Comparison of simulated limiting membrane potentials with those predicted by the Donnan equation for a type 1 membrane system for a range of concentration ratios, c*L /c*R, with DA′ /DX′ ) 1.

A study of the effects of the rate of diffusion of each ion was performed by varying the diffusion coefficient ratio, DA′ / DX′ , from 0.5 to 3 while keeping the concentration ratio, c*L /c*R constant. It was observed that the ratio DA′ /DX′ is an exclusively dynamic parameter and thus has absolutely no effect on the steady-state configuration of the membrane system, merely determining the rate at which steady state is achieved. This is contrary to the behavior at short time where the rate of growth is dependent only on the relative rates of diffusion. Additional plots showing the effects of varying DA′ /DX′ and c*L /c*R on the system dynamics and the rate of attainment of steady state are given in the Supporting Information. 3.1.3. Perturbation Due to Unequal Diffusion Coefficients. If the diffusion coefficients of the two ions are not equal, it is observed that the extent of the nonzero electric field continues to increase indefinitely (so long as the method of calculating χ is sensitive enough to nonelectroneutrality). Directly adjacent to the membrane, on the left-hand side, migration of X- due to the electric field generated by diffusion of A+ through the membrane leads to a slight depletion (of X-) relative to the bulk concentration, with the excess ions accumulating slightly further to the left. The strength of the electric field prevents the concentration of X- from relaxing back to its initial equilibrium position, and so it diffuses out into the bulk solution, as does a small amount of A+ that has diffused from the right. From Figure 4, it can be seen that, at a typical distance of about 200 nm (about 20 dimensionless units in this case) from the membrane, the concentration profiles of both species are coincident. At this distance, the electric field remains small and the flux is dominated by diffusion:

| ∂c∂X | . |z c ∂X∂θ | q

q q

(3.7)

Therefore, when the diffusion coefficients are unequal, one species diffuses into the bulk faster than the other and a small charge separation is observed. If D′A ) D′X, however, no further charge separation is observed, since the two species diffuse at the same rate. The same argument also applies on the righthand side of the membrane. The magnitude of the electric field due to this effect is negligible compared to that at the membrane. Furthermore, this diffuse charge separation accounts for significantly less than 1% of the total system potential difference for all systems studied, and may therefore be ignored for most purposes.

Figure 4. Concentration profiles for a type 1 membrane system with c*L /c*R ) 0.1 and DA′ /DX′ ) 4.

Figure 5. Evolution of the potential difference (∆θMem), maximum electric field (ξmax), and junction extent (χ) of a type 2 membrane system with DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

3.2. Type 2 Membrane Systems. In the type 2 membrane case, a common ion is shared between each of the two reservoirs and both solutions have the same concentration. Here, we shall consider systems in which only the common anion is blocked by the membrane, i.e., systems of the type

AX(C ) C*)

|

BX(C ) C*) X

The investigation was conducted by varying the diffusion coefficient ratios DA′ /DX′ and DB′ /DX′ through a range of values from 0.33 to 3. Most of the figures are based on the exemplar system, DA′ /DX′ ) 3; DB′ /DX′ ) 0.33, but the conclusions drawn are equally applicable to every case investigated. Figure 5 displays the evolution of the dynamic properties of such a system. 3.2.1. Comparison with Existing Models. The steady-state membrane potential is often examined with reference to the Goldman equation3 for this system:

( )

∆θGold ) ln

DA′ DB′

(3.8)

which is based on the unphysical approximation of constant (i.e., bulk) concentration at the membrane boundary (a Dirichlet

10768

J. Phys. Chem. B, Vol. 114, No. 33, 2010

Ward et al.

Figure 6. Concentration profiles of all species in the area close to the membrane at long time for a type 2 system: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

Figure 7. Comparison of Goldman and Henderson models for ∆θMem with simulated results for a range of type 2 systems.

boundary condition). The Goldman equation considers only the permeant species, and the electric field is assumed to be linear between either side of a membrane of finite thickness. Since concentrations outside this membrane are considered uniform under this model, the only source of potential difference must necessarily arise from the interior of the membrane. However, it is obvious from a plot of the simulated concentration profiles at long time (Figure 6) that the Dirichlet boundary conditions are a gross simplification, if not indeed altogether inaccurate, as the concentration inhomogeneity extends well beyond the vicinity of the membrane. Another relevant model is the Henderson equation21 for the fully permeable system:

(

∆θHend ) ln

DA′ + DX′ DB′ + DX′

)

(3.9)

which is used as an approximation to a steady-state free liquid junction potential, i.e., unconstrained diffusion of all ions. The same type of boundary conditions apply in the derivation, but this time at some arbitrarily positioned outer boundaries of a “junction layer”. The utility of eqs 3.8 and 3.9 as applied to a membrane system have been investigated to some extent by Perram et al.,4 and the Henderson equation has been discussed by Dickinson et al.5 and Ward et al.22 A comparison of the simulated ∆θMem with that predicted by both the Henderson and Goldman equations is shown in Figure 7. It can be seen that the simulated result lies somewhere between the two models in all cases, with the scatter from linearity qualitatively resembling a less pronounced version of the Henderson case. A table displaying the simulated potential differences along with those predicted by the Goldman and Henderson models for every case examined is available in the Supporting Information 3.2.2. OWerall BehaWiorsTwo-Component Potential. Figure 5 displays the evolution of the dynamic properties of a typical type 2 system. At short time, both the field and potential difference are monotonically increasing, accompanied by a growing extent of the region of nonzero electric field. At long time, it appears from this plot that, while the system potential difference and maximum electric field (which always occurs at the location of the membrane) both tend toward a steady state, the spatial extent of the nonzero field continues to increase

Figure 8. Electric field profiles at various times for a type 2 membrane system with DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

monotonically. That an expanding layer of constant charge separation cannot maintain a steady-state potential difference has been discussed in the work of Dickinson et al.,5 and this behavior would appear, at first sight, to be in violation of Maxwell’s equations. However, a closer examination will reveal that this plot disguises much of the behavior. In Figures 8 and 9, it can be seen that in the central spatial region, close to the membrane, the potential profiles are very nearly coincident at all long times; i.e., the electric field in this area is static on the time scale studied. Farther from the membrane, however, the potential “spreads out”: the electric field decays in time. It is therefore clear that the total system potential difference arises from two separate components. One is a local component, analogous to the Donnan-type potential observed in the type 1 case, for which the electric field is unchanging (on the time scale studied) and which remains close to the membrane ((25 nm for this typical case). The other is a diffuse component where the charge separation, and hence electric field, is dynamically relaxing. In Figure 5, the maximum electric field is associated only with the former, whereas the extent of nonzero electric field is associated only with the latter. We may write

∆θMem ) ∆θLocal + ∆θDiffuse

(3.10)

Dynamic Theory of Membrane Potentials

J. Phys. Chem. B, Vol. 114, No. 33, 2010 10769

Figure 9. Potential profiles at various times for a type 2 membrane system with DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3. The indicated region in the center may be considered to arise from the Donnan-type potential, whereas the other regions arise from the LJP-type potential.

Figure 10. Evolution of the system potential difference for a free liquid junction and a membrane system with identical chemical composition: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

Of particular note is that the electric field in the area outside of the central, Donnan-type region contributes significantly to the overall system potential difference. 3.2.3. Short Time BehaWior. Mathematical analysis conducted in an analogous fashion to that of the type 1 case revealed that at short time

ξmax )

√DA′ - √DB′ √π

· √τ

(3.11)

Figure 11. Evolution of the maximum electric field for a free liquid junction and a membrane system with identical chemical composition: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

Figure 12. Total charge flux (jA + jB - jX) at various times for a type 2 system: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

in Figures 10 and 11, comparisons of the evolution of system potential and maximum electric field, respectively), which at short time behaves identically. 3.2.4. Intermediate to Long Time. As the pair of double layers forms, the developing electric field significantly alters the behavior with respect to that of a free liquid junction. The flux profiles of each of the positively charged species (A+ and B+) resemble slightly asymmetric bell curves with the greatest flux at the location of the membrane. Despite the fact that these species may have quite different diffusion coefficients, the total positive flux, defined as

j+ ) jA + jB

and

∆θMem ) (DA′ - DB′ ) · τ

(3.12)

which was consistent with the behavior observed from simulation in all cases (figures are available in the Supporting Information). Since the impermeant species is of initially uniform concentration, its distribution is not perturbed as long as the environment remains approximately electroneutral. The membrane then has no effect on the short time behavior of a type 2 system, and the above equations are identical to those derived for the case of a type 2 free liquid junction5 (as shown

(3.13)

(which is approximately equal to the total charge flux shown in Figure 12, since the flux of X- is small) reaches a local minimum at the membrane. The increasing electric field leads to migration that increases the magnitude of the overall flux of the slower species, as it is in the same direction as its diffusion, but decreases that of the faster species, as it opposes its diffusion (Figure 13). Since the electric field is largest at the membrane, due to the presence of the double layers, and since the diffusive flux is initially largest here, the equilibration of flux is more rapid here than further from the membrane. At the time when the electric field at the membrane reaches its maximum, the

10770

J. Phys. Chem. B, Vol. 114, No. 33, 2010

Ward et al.

Figure 13. Direction and approximate magnitude of the diffusional and migrational fluxes for the positive species, A and B, at an intermediate time where DA′ > DB′ . A larger arrow indicates a larger flux.

Figure 14. Evolution of electric field and positive flux (jA + jB) for a type 2 membrane system: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3. The dotted line indicates the time of greatest electric field.

Figure 15. Evolution of flux at the membrane of species A and B for a type 2 membrane system with DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3. It can be seen that the magnitude of the flux of B becomes slightly greater than that of A.

Figure 16. Flux profiles of species X- at logarithmic times for a type 2 system: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

flux of B+ becomes slightly greater than that of A+ (as can be seen in Figures 14 and 15) but then relaxes slowly to a condition of equality:

||

jB f1 jA

as τ f ∞

(3.14)

The flux of X- (shown in Figure 16) does not contribute significantly to the total charge flux, though it too tends to zero close to the membrane (and is zero by definition at the membrane). For the faster positive species (A+), the diffusive flux decreases with time as the concentration gradient becomes shallower, and the migrational flux increases (in the opposite direction) with time as the electric field develops, until both become balanced. For the slower positive species (B+), both flux components are initially in the same direction, leading to an accumulation on the opposite side of the membrane to the reservoir in which it started. This in turn causes an inversion in its concentration gradient local to the membrane (Figure 17), which results in a diffusive component that opposes the migratory one.

Figure 17. Concentration profiles of species B+ at logarithmic times for a type 2 system: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

In the case of the negative species (X-), migration leads to a local accumulation on the side of A+ and depletion on the opposite side, generating a pair of concentration gradients and thus opposing diffusive fluxes, as shown in Figure 18. Therefore, at long time, the diffusive flux of each species at and close to the membrane is precisely balanced by its own opposite migrational flux:

Dynamic Theory of Membrane Potentials

∂cq ∂θ ) -zqcq ∂X ∂X

for all species, q

J. Phys. Chem. B, Vol. 114, No. 33, 2010 10771

(3.15)

and a local pseudosteady state is achieved. 3.2.5. The Impossibility of Discharge. The short time electric field induces two adjacent double layers at the membrane, of opposite charge separation. On one side of the membrane, the negative species is accumulated and the positive species are depleted with respect to the diffusing concentration gradient; on the other, the negative species is accumulated and positive species are depleted. Since each double layer has a nonzero local charge, it induces Coulombic forces which maintain the opposite double layer, and so the two double layers exist in a synergistic pseudosteady state. When jB > jA (eq 3.14), a mechanism is provided for some of the total positive charge to return to the left-hand side, which partially discharges the system and explains the behavior of the diffuse component. However, the only means of relaxation of the local component is the diffusion of anions through the membrane, where its concentration gradient is extremely large, but this diffusion is formally forbidden by the definition of the system, i.e., the assumption that at the membrane, the anion is impermeant. In a free liquid junction, diffusion of anions through the membrane is continuous and zones of accumulation and depletion discharge alongside the discharge associated with the change of direction of the net positive flux, as is confirmed by detailed examination of the data from simulations reported in Dickinson et al.5 Therefore, even though the negative species is exposed to an externally induced electric field at short time in both systems, no local potential of the Donnan type is observed in a free liquid junction, whereas, under the same conditions, the presence of a blocking membrane is both a necessary and sufficient condition for such a localized pseudosteady-state potential to occur. 3.2.6. The Nature of the “Local” Layer. Although its concentration profiles remain constant at long time, the “local layer” is a dynamic property and is not in agreement with any equilibrium for the system. It can be demonstrated analytically that a finite steady state cannot be attained for the type 2 case, as it can for the type 1 case (as detailed in the Supporting Information). As is clear from Figure 9, however, this component of the potential difference has nearly constant magnitude and is also confined to a constant region of space, so it is distinct from the “dynamic” layer for which the extent of nonzero electric field grows with the diffusion layers.

Figure 18. Concentration profile of species X- at τ ) 0.001 for a type 2 system: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

Figure 19. Concentration profiles of all species in the area close to the membrane at very long time, τ ) 106, for a type 2 system: DA′ /DX′ ) 3 and DB′ /DX′ ) 1/3.

To determine quite how fixed the layer remains, further simulations for the exemplar system above were run to very long time (τ ) 106, corresponding to the order of 10 ms for typical concentrations). No change was observed in the concentration profiles very close to the junction (|X| < (1), where the two double layers at the membrane remain unperturbed, although the concentration profiles outside this region relax, and change considerably in magnitude over logarithmic time as the diffusion layers grow. There is no evidence whatsoever of the double layers relaxing as the external concentration gradients do so, and the magnitude of the junction potential due to the local layer is also unaffected. The concentration profiles in a range close to the membrane (|X| < (100) are shown for an unconstrained system at τ ) 106 in Figure 19, with evident similarity to Figure 6. It is clear that the local layer shares many properties with Planck’s junction layer of finite thickness,12 as it is constant in size, its outer concentrations become constant, and the concentration profiles in its interior are at a local pseudosteady state. The local layer differs markedly from Planck’s junction layer, however, in that the concentrations at the edge of the local layer are not equivalent to the bulk concentrations, and the junction potential of the system as a whole has a significant dynamic component which is due to charge separation external to this layer. Since the dynamic component of the junction requires a continual flux of A+ and B+ through the membrane, the local layer cannot be precisely at steady state; at long time, it is most accurate to say that it is not observably distinct from a steady state. It is therefore our extrapolation that if the system evolves indefinitely, this local component will remain, although it is not possible to determine this directly as simulations cannot be run to infinite time. What is clear is that the local potential remains constant at least up until a real time of 10 ms for typical values of diffusion coefficient and concentration, when the diffusion layer will extend at least ≈5 µm from the membrane in either direction. In a real system, it is extremely likely that some boundary (convectional or otherwise) will be encountered before this point, thus collapsing the system and the potential difference. This is only conditional, however, on the provision of infinite solution to diffuse through the membrane. The results of the simulation with finite Neumann boundaries (e.g., the walls of a container) are shown in Figures 20 and 21, which are the evolution of the system potential difference and the evolution of the concentration of species X-, respectively. It

10772

J. Phys. Chem. B, Vol. 114, No. 33, 2010

Figure 20. Evolution of membrane potential, ∆θMem, of a system of infinite extent and a constrained system of total extent X ) 20 (≈100 nm).

Figure 21. Evolution of the concentration profile of species X- for a constrained system of total extent X ) 20 (≈100 nm).

can clearly be seen that the imposition of an insulating boundary to the space results in the eventual total collapse of the local layer along with its potential difference, with the system becoming homogeneous at long, but finite, time. A more expansive system leads to a greater time before collapse as the collapse only begins when the diffusion layer reaches the imposed boundary. Note that the diffuse component of the potential difference also collapses if infinite solution is not available, as was found in a related system by Perram et al.4 3.2.7. OWerWiew. The superficially paradoxical behavior (Figure 5) is explained by the fact that one component of the electric field discharges via a growing diffuse layer, accounting for the increasing extent of the nonzero electric field, while a second, localized component accounts for the constant magnitude of maximum field observed at the membrane. Consequently, while a system of infinite extent cannot reach an overall steady state in finite time, a local pseudosteady state in the region close to the membrane can be attained relatively rapidly. On simulation of a system of finite extent (but with zero flux rather than constant concentration boundary conditions), the steadystate potential difference is zero as expected. Crucially, for the infinite system, the dynamic, discharging component nonetheless generates a constant potential difference across the membrane, since the region of nonzero electric field is expanding and so the distance of charge separation is increasing at a rate that balances its decreasing magnitude. This

Ward et al.

Figure 22. Approximate local and diffuse components of ∆θMem with lines of best fit for a range of type 2 systems.

thesis was presented and justified by Dickinson et al.5 as the mechanism for steady-state potential difference at a free, totally permeable liquid junction. 3.2.8. Separating the Components. By making the reasonable approximation that, at very long time, the central peak in the electric field corresponds entirely to the local component of the potential, it is a simple task to determine its spatial extent. This was not found to vary significantly with the diffusion coefficient ratio and was approximately 50 nm in dimensional units for a concentration of X- of 1 mM (a typical value). It is therefore possible to resolve values to within a few % for the potential difference resulting from each component. Figure 22 shows the values for the local (∆θLocal) and diffuse (∆θDiffuse) components, respectively. It may be expected that the diffuse component should behave in a similar manner to the Henderson model, whereas the local component should be closer in nature to the Goldman model. It can be seen (Figure 22) that the values of the diffuse component show greater scatter from the line of best fit than do those of the local component which correlates with the behavior of the respective models, as shown in Figure 7. From this, it can be seen that each component comprises a similar proportion of the potential regardless of the relative diffusion coefficients. 4. Conclusions A detailed view of the dynamic evolution of both “type 1” and “type 2” infinitesimal membrane systems has been presented. The time evolution of both systems, and thus their configurations at long time, has been elucidated with reference to the mass transport of charged ions. It has been confirmed that, for the type 1 system, the Goldman equation (which is equivalent to the Donnan equation in this case) provides an excellent approximation of the steady-state value of the membrane potential. For the type 2 system, however, it has been shown that, contrary to classical understanding, the potential difference consists of two components. One is a localized, Donnan-type layer that remains at pseudosteady state very close to the membrane; the other is a dynamically expanding diffuse component, with properties similar to those of a free liquid junction potential, that discharges at constant potential difference, with the contributing electric field extending increasingly far from the membrane. It has therefore been shown that the Goldman equation does not adequately describe the τ f ∞

Dynamic Theory of Membrane Potentials membrane potential for a thin membrane in the type 2 case; the approximation of a confined steady-state membrane layer3,10 is misleading and may introduce inaccuracies, especially where such a membrane is considered in the context of a more complex system. Acknowledgment. E.J.F.D. thanks St John’s College, Oxford, for funding. Supporting Information Available: (1) Details of the numerical methods used to solve the Nernst-Planck-Poisson equations dynamically; (2) details of the mathematics of the iterative Newton-Raphson method used to numerically solve the nonlinear equations in this paper; (3) mathematical analysis of a type 1 system as τ f 0; (4) mathematical analysis of a type 1 system as τ f ∞; (5) an analytical classification of membrane and liquid junction systems; (6) plots showing agreement between the short and long time mathematical analyses and the results of simulation for type 1; (7) plots showing further information about the dynamic and long time behavior of type 1 systems; (8) additional plots for the type 2 case; (9) a table of the simulated membrane potentials along with those predicted by the Goldman and Henderson equations for a range of systems. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Koryta, J.; Sˇtulı´k, K. Ion SelectiVe Electrodes, 2nd ed.; Cambridge University Press: Cambridge, U.K., 2009.

J. Phys. Chem. B, Vol. 114, No. 33, 2010 10773 (2) Berg, J. M.; Tymoczko, J. L.; Stryer, L. Biochemistry, 6th ed.; W. H. Freeman & Co.: New York, 2006. (3) Goldman, D. E. J. Gen. Physiol. 1943, 27, 37–60. (4) Perram, J. W.; Stiles, P. J. Phys. Chem. Chem. Phys. 2006, 8, 4200– 4213. (5) Dickinson, E. J. F.; Freitag, L.; Compton, R. G. J. Phys. Chem. B 2010, 114, 187–197. (6) Morf, W. E. Anal. Chem. 1977, 49, 810–813. (7) Brumleve, T. R.; Buck, R. P. J. Electroanal. Chem. 1978, 90, 1– 31. (8) Mafe´, S.; , J. P.; Aguilella, V. M. J. Phys. Chem. 1986, 90, 6045– 6050. (9) Sokalski, T.; Lewenstam, A. Electrochem. Commun. 2001, 3, 107– 112. (10) Planck, M. Wied. Ann. 1890, 39, 161–186. (11) Nernst, W. H. Z. Phys. Chem. 1889, 4, 165. (12) Planck, M. Wied. Ann. 1890, 40, 561–576. (13) Cohen, H.; Cooley, J. W. Biophys. J. 1965, 5, 145–162. (14) Hafemann, D. R. J. Phys. Chem. 1965, 69, 4226–4231. (15) Filipek, R.; Szyszkiewicz-Warzecha, K.; Bozek, B.; Danielewski, M.; Lewenstam, A. Defect Diffus. Forum 2009, 283-286, 487–493. (16) Jackson, J. S. J. Phys. Chem. 1974, 78, 2060–2064. (17) Sokalski, T.; Lingenfelter, P.; Lewenstam, A. J. Phys. Chem. B 2003, 107, 2443–2452. (18) Compton, R. G.; Banks, C. E. Understanding Voltammetry; World Scientific: Singapore, 2007. (19) Bard, A. J.; Faulkner, E. A. Electrochemical Methods: Fundamentals and Applications, 2nd ed.; John Wiley & Sons: New York, 2001. (20) Donnan, F. G.; Guggenheim, E. A. Z. Phys. Chem. 1932, A162, 346–360. (21) Henderson, P. Z. Phys. Chem. 1907, 59, 118–127. (22) Ward, K. R.; Dickinson, E. J. F.; Compton, R. G. J. Phys. Chem. B 2010, 114, 4521–4528.

JP102599J