J. Phys. Chem. 1995,99, 12030-12040
12030
Integral Weak Diffusion and Diffusion Approximations Applied to Ion Transport through Biological Ion Channels Alexander Syganow and Eberhard von Kitzing* Abteilung Zellphysiologie, Max-Planck-Institut f i r medizinische Forschung, Postfach 10 38 20, 0-69028 Heidelberg, Germany Received: January 9, 1995; In Final Form: April 17, 1995@
In this article a theory is presented to calculate integral properties of biological ion channels (like currentvoltage and conductance-concentration relations). The qualitative form of these relations predicted by the theory agrees well with data measured in experiments. For instance, the saturation of the channel conductance with increasing external ion concentration is predicted for a class of ion channels (as, for instance, found for the gramicidin A,' acetylcholine receptors,*s3NMDA? and sarcoplasmic reticulum channels5). In contrast to commonly used approaches such as the Eyring rate theory, this method is directly related to physical parameters of the ion channel such as the channel length and diameter, dielectric constant, ionic mobility, and minimal ionic concentration inside the channel. The theory starts from Nernst-Planck and Poisson equations. Using the method of phase trajectory (as proposed by Schottky) and the regional approximation, rather general expressions can be derived for integral channel quantities in the drift limit (IVl > kBT/eo)in the presence of multiple ionic species. The theory predicts two typical types of conductance-concentration relations found experimentally: a monotone saturating conductance and a maximum in the conductance. The realized type of relation depends on the minimal ionic concentration inside the channel. In the present form the theory is restricted to narrow ion channels where the length exceeds its diameter. The ions are assumed to behave like structureless point charges at not too high ionic concentration.
1. Introduction Theoretical approaches for understanding ion transport through biological channels requires (a) a model concerning the threedimensional structure of the channel and (b) a physical mechanism of how the ion moves through the channel. In early studies of ion permeation through membranes, Goldman: and later Hodgkin and Katz,' assumed that a single ion would pass the channel independent of the presence of others. The ionic density in the channel was expected to be constant along the channel and so low that ion-ion interactions could be neglected. In addition, the external voltage was assumed to drop linearly along the channel. These assumptions led to the well-known Goldman-Hodgkin-Katz (GHK) current equation (for an introduction see ref 8). Among experimentalists in the field of ion channels in biological membranes a very popular method to describe ion permeation is the Eyring or Kramer's rate theory.8-'0 In such an approach, the structural information is the number of binding sites, their relative position in the external electric potential, the respective ion-binding energies, and the height and position of the barriers between the different sites and the bulk solutions. Eyring rate theory is used to describe the motion of ions between the sites. Usually, the energy wells and barriers are assumed to be constant. They do not depend on the external voltage or ion concentration. Under this assumption these models largely neglect long-range electrostatic ion-ion and ion-channel interactions. The external potential is assumed to drop linearly across the membrane. One major advantage of Eyring rate theory is that it easily explains the saturation of the channel conductance with increasing ionic strength.
* To whom correspondence should be addressed. Tel.: +49 6221 486 467; FAX: f 4 9 6221 486 459. Internet:
[email protected] Abstract published in Advance ACS Absrrucrs, July 15, 1995. @
0022-365419512099-12030$09.0010
An alternative approach to model ion permeation is diffusion theories. In this case, the structural information is encoded in an external free energy profile which may result either from pure speculations or from three-dimensional models of the considered channel (for instance, the gramicidin channel,' ' - I 3 the acetylcholine receptor channel,I4-l9 or other proteins assumed to allow for ion permeation20-22). This approach has the disadvantage that an atomic structure is available only for the gramicidin ~ h a n n e l . l I - ~Up ~ to now it has generally not been possible to decide whether a differencebetween experiment and theory results from shortcomings in the theory or in the structural model. Either diffusion theories start from the Fokker-Planck equation10s23-28 and derive the Nemst-Planck equations from there, or the Nemst-Planck equation itself is taken as the starting The various approaches using the NemstPlanck equation differ in dealing with long-range electrostatic ion-ion interactions. If they are neglected, one essentially obtains a generalization of the GHK a p p r ~ x i m a t i o n . ~For .~~ calculating the mean electrostatic forces between ions two approaches were proposed in literature: (a) the PoissonBoltzmann30*32.35-38 and (b) the Vlassov approach.24,26.39-42 Because diffusion models can account for the physical properties of long-range ion-ion interactions, hybrid models containing a diffusion and an Eyring rate part have been p r o p o ~ e d . ~For ',~~ a general survey of theories used to describe the properties of ion channels see the essay of EisenbergS4 To study the distribution of charge carriers or to calculate current-voltage relations of semiconductors, solutions for the Nemst-Planck equation together with Poisson's equation have been studied since the 1930s (for reviews see refs 45-47). A commonly used procedure to obtain approximate solutions is the method of regional ~ p p r o x i m a t i o n .Using ~ ~ explicit models for the charge carrier distribution, integral characteristics such as resistance or current-voltage relations may be obtained in 0 1995 American Chemical Society
Ion Transport through Biological Ion Channels
J. Phys. Chem., Vol. 99, No. 31, 1995 12031
analytical or numerical form. In this way, a solution for the charge carrier distribution and current-voltage relation assuming constant electric field along the material was found in 1941.48 This solution was later rediscovered for describing the properties of ion channels6%’and is now known in electrophysiology as the GHK theory.* An approach which does not require specific models for the charge carrier distribution was proposed by Schottky49and has only recently been successfully ~ t i l i z e d . ~ O -This ~ ~ method is designated the phase trajectory approach. Because it does not require a priori assumptions about the charge carrier distribution, this method is sometimes labeled “model free”. In this article, the methods of regional approximation and phase trajectory are used to calculate integral properties of biological ion channels using rather general assumptions. Although these theories were originally concerned with conductivities in semiconductors, many results also apply to understanding ion permeation through ion channels in biological membranes .
2. Theory The starting point for describing integral properties of ion channels is the Nemst-Planck equation:8.’0
j =p[qEn
-k ~ m n ]
(1)
where j is the electric current density, p is the ionic mobility, q is the ionic charge, E is the electric field due to fixed charges and mobile ions, n is the ion number density (this corresponds to ionic concentration in the more biologically oriented literature), kBT is Boltzmann’s constant times the absolute temperature, and v is the nabla operator. In the following, only the one-dimensional form of the Nemst-Planck equation is used. The electric field is determined by the Poisson equation:
where 6 is the dielectric constant and e is the electric charge density of the ions and possible other charges in the system. This approach requires that the ions behave like structureless point charges and that their motion be govemed by many molecular collisions to justify the Nemst approach. The mean field approximation breaks down at too large ionic number densities and charge densities. Such a situation is designated as space-charge limited. These two equation represent a set of nonlinearly coupled equations designated here as the Nemst-Planck- Vlassov approa~h.~ Only ~ , ~rarely are analytic solutions known for this system. In general, it is simpler to obtain integral properties of this equation than a detailed local description. In the case of ion channels, only integral quantities can be obtained experimentally. It is therefore of interest to study these integral properties theoretically. 2.1. The Model System. Figure 1 shows the model system used in this article. There is a membrane containing the ion channel. Both sides of the membrane contact an aqueous electrolyte reservoir. The dielectric constant of the electrolyte is close to 80, whereas that of the membrane is much lower. The membrane acts as a barrier for the ions because it reduces their Bom energy of h y d r a t i ~ n . ~ ~ It - ~ is ’ assumed that the channel is one-dimensional; that is, its diameter d is smaller than its length L: dz >j o 11
for the diffusion domain
+ E8 *dnl+’”
for the drift domain
(1 1) (12) (13)
where j o is the current density at the boundary between a closeto-equilibrium system bl > j o and is given by
. Jo=
Figure 3. Number density of positive ions n shown along the channel axis for positive outside voltage (continuous line) and the equilibrium density (dashed line). One can distinguish the following three regions: (a) the equilibrium domain at the side toward the external solution (eq 11) with the typical phase trajectory of K,, = KQ (eq 8), (b) the drift domain inside the channel (eq 13) with the typical phase trajectory of K E (eq lo), (c) the diffusion domain at the side toward the internal solution (eq 12) with the typical phase trajectory of K , = K D (eq 9) (see text).
= eo(n - n~))?’ In such approaches the problem was always how to obtain physically reasonable approximations for the phase trajectory K across the boundaries of the regions (i.e., between equilibrium and drift domain; see Figure 2). 2.2.2. The Voltage Drop Across the Membrane. Experimentally, the external voltage V is applied across membranes containing ion channels, and the resulting ionic current is measured.8 From eq 4 the external voltage V as a function of the current density j can be obtained:
kBT
dn +-ln-
ni
(14)
n 2 , Q , ] 112
E
It should be noted that j o is a function of the local ionic density n and consequently a function of the coordinatez. The boundary between a close-to-equilibrium situation and a steady-state situation occurs at different j o values at different places in the channel. With the exception of the article of to the knowledge of the authors, the concept of regional approximations has only been used in the context of approximating the charge density Q(n)by means of an explicit function (i.e. e ( n )
where no and n, are the ion number densities outside and inside the cell, respectively (see Figure 1). The integral term on the right-hand side of eq 14 is the voltage drop across the ohmic resistance of the channel, whereas the second one is the wellknown Nernst potential which results from possible different ionic number densities on both sides of the membrane. The external voltage V as given in eq 14 is reformulated to isolate contributions from different physical origins. First the single integral is divided into two parts, integrating from the minimal ionic density nE to the solution densities ni and no:
J. Phys. Chem., Vol. 99, No. 31, I995 12033
Ion Transport through Biological Ion Channels
The dependence of K on j and n is left out for the readability of the equations. As can be seen from Figure 2, the sign of the phase trajectory yields K < 0 for nE In 5 no and K > 0 for nE In Ini. Replacing 1 / using ~ eq 6, the voltage becomes
The first integral X I in this equation is reformulated via integration by parts:
which together with eq 7 gives
With the inverse phase trajectory 1 / from ~ eq 6, the channel length L becomes
The first term on the right-hand side of this equation is the pure "drift contribution" LE of the channel length L. Similar to the external voltage V the apparent length L of the channel is split into the drift length LE and a rest 15%:
LE would equal ZE - zo, if the drift domain would tend to zo. Because this is not the case, the often small corrections are given by L* including the interval ZE to zi (see Figure 2):
with 6 = n!!Q
e"
X2
is reevaluated using integration by parts:
with K from eq 7 Now the extemal voltage has obtained the form
Here VEV) is the so-called drijl ~ o l t a g e . ~It' would give the exact voltage between zo and ZE if the drift condition 13 holds through this range. The contributions in the equilibrium domain close to the outer membrane wall (see Figure 2 ) and in the diffusion domain close to the inner membrane wall, however, are small at large voltages (\VI > kBT/eo) for long channels. Consequently, V E ~is) a good approximationof the total voltage. The corrections required to obtain exact results are given by V4j):
The physical interpretation of the various terms in V* are given below. It should be noted that the extemal voltage as given in eq 15 is still the exact voltage according to the Nemst-PlanckVlassov approach (eqs 1 and 2). 2.2.3. The Apparent Length of the ion Channel. Using the relation between the ion number density n and the channel coordinate z (see eq 4) leads to a relation between the channel length L and an integral property of the ion channel:
leading to the following expression for L*:
Again, eqs 18 and 19 represent an exact result from the NemstPlanck-Vlassov approach. 2.2.4. The Weak integral Dimsion Limit. In this section the extemal voltage Vis calculated in the limit of large voltages IVI >> kBT/q. In this case the drift approximation dominates the integral quantities such as extemal voltage V for long ion channels ( L >> L*). This limit is designated the dri9 limit. For short channels the situation depends on the system parameters. In Figures 2 and 3 the three different domains are indicated. For zo I / K ~< z < zi the equilibrium condition 11 is violated, ~ ~ resulting in j >>jo. In the range of zo < z < zo 1 / the system deviates only slightly from the equilibrium situation (see
+
+
Syganow and von Kitzing
12034 J. Phys. Chem., Vol. 99, No. 31, 1995 I
is the drift voltage between ZE and small and may often be neglected.
small voltage Ohmic law
This term is generally
is the Ohmic voltage drop between zo and zo
, 0-
zi.
I1
+11~~.
V .diffusion:regime
drilt
[email protected] with weak diffusion
Figure 4. Current-voltage relation. For relatively small energy barriers two domains can be discovered in the current-voltage curve: the diffusion domain for small voltages (IVl < k&"/eo) and the drift domain ([VI h kBT/eo), including contributions from weak diffusion. The linear ohmic law for large voltages is given by the line with short dashes. The diffusion voltage VD obtained from the large voltage limit should not be mistaken for the reversal potential VR, which is obtained from the low-voltage approximation (see text). The superlinear behavior of the current-voltage curve corresponds to the situation shown in Figures 2 and 3, where the equilibrium ion number density inside the channel is assumed to be below the extemal number density, nc < ns. In the other case, where the limiting density nc exceeds the solution density ns, a sublinear current-voltage relation is predicted.
+
Figure 3). In the intermediate range zo 1 1 > nE 2 nc the conductance approaches %*. 2.3. Channels Containing Multiple Ionic Species. Up to this point only the monopolar situation has been considered where the electric current is assumed to be dominated by the motion of a single ionic species. In the following section, the results are generalized to the case of multiple ionic species in
Syganow and von Kitzing
12036 J. Phys. Chem., Vol. 99, No. 31, 1995
The voltage V , the length L, and the resistance R are given by
$16
4
V=j j-do f3K
+ kBT -j eo
4 -do 0
10:
Here the voltage is composed of the ohmic part proportional to the current density j and a part which is designated the reversal potential VRin electrophysiologya and is defined to be the voltage at zero total current density j :
5:
r
-
.
2
4
6
8
10
*
n,lnc
Figure 5. Ion channel conductance in the weak integral diffusion drift limit. The curves are drawn for the case that the minimal ionic density nE equals the long channel limiting ionic density nc. For nE > nc, the actual values change, but not the qualitative behavior of the channel conductance. The ratio IclL (see text) ranges from 2 for the curve with the lowest maximum to 3 in steps of 115. The position of the maximum is independent of the ratio IclL and given by no x 4 nc. In the limit of large extemal ionic number densities no >> nc the conductivity ratio Z$Z$ approaches unity.
solution, which allows us to study more realistic cases. Let us consider K ionic species with number densities nu, mobilities pa, and charges qa. The total specific conductivity a(z) inside the channel and the total charge density e(z)are given by K
42)
= C l q , l ~ a n a ( Z ) and r=1
K
e(z>= [Cq,n,(z)
- efl(z)l
r=l
where N(z) is the number density of fixed charges inside the channel (Le. charged amino acid side chains). The multispecies Nernst-Planck equation now takes the form K
kBT d b j = DE - -- with I3 = &&,n,(z) eo dz i= I
According to the mean value theorem for integrals, there exists some intermediate value q5 for #(a) for which the reversal potential yields
This result is interesting because it shows that the reversal potential obtained from GHK theory8holds under rather general tonditions. For many practical applications q5 are close to unity. 4 becomes ion number density dependent if the current density results from non-negligible contributions of ionic species with different signs of the charge q, or if there is considerable fluxflux coupling. 2.3.1. The Integral Weak Diffusion Limit. With these equations the contributions to the apparent channel length in the integral weak diffusion limit are given by
where eo is the absolute value of the electronic charge and 6 is the specific conductivity of the positive ions a, minus the specific conductivity of the negative ions a,. In the following, the same terms are calculated as for the monopolar approximation. But now, only the major results are given. The phase trajectory is defined by (see eq 3)
(39)
and the different voltage contributions:
The respective differential equation yields (33)
From this equation expressions for the phase trajectory K and its inverse K-' are evaluated:
(34) The phase trajectory in the three domains yields
2.3.2. The Channel Conductance as a Function of Extemal Ion Number Density. The conductance of a channel as a function of the extemal ion number density depends critically on the minimal specific conductivity UE in the channel in relation to the specific conductivity ac at which the charge density becomes zero. As shown below, one can distinguish two important limiting situations: (a) UE >> ac and (b) OE 2 OC. At low extemal voltages 1 VI OC; for OE k ac the channels are long. Together with the ohmic length (see eq 39) the screening length ~ D His one of the characteristic quantities determining the mechanism of ion conduction through the channel. 2.3.3. The Ionic Conductance of Short Channels. The equilibrium phase trajectory KQ at the boundary between the ion channel and solution is given by eq 35. For symmetric electrolytes the internal specific conductivity ai equals the external specific conductivity uo ss ai; the length L and the ohmic resistance R are obtained using eq 32: L=f’dz=2p -0
- fU2
-0
At large external ion number densities (ao>> OE) the arctan of the length L tends to n/2, leading to the following approximate relation for the charge density at minimal specific conductivity:
provided that @(a)is largely constant in the region of low specific conductivity. This equation is used to estimate the maximal specific conductivity a ,, in the channel. Because of UE z-ac (condition for a short channel), the charge density is approximately linear in the specific conductivity: @(a)% a/p. Consequently for qh = 1 the maximal specific conductivity becomes
The existence of a maximal specific conductivity amax in the channel is responsible for the saturation of the channel conductance found in many biological ion channels.8 Dividing the resistance R by the length L allows the yet unknown @E and @E to be eliminated, yielding the conductance Z for short channels as a function of external ion number densities:
dz=
, ]’”/ ;[-+-,[ -Emax
As long as the integral conductivity is dominated by a single ionic species, #(a) is nearly constant, and thus @(a) @(uE) @E. Consequently, L and R simplify to
*
The inner integral in the equation for L and R can be estimated
z=
1
1
1]112
with
arctan
Oo aE
This form of equation is interesting because it predicts the often found saturation of the ionic conductivity with increasing ion number density.8 For the case that ffE and a, are comparable, the function behaves like a Michaelis-Menten curve describing the activity of enzymes and often used to describe the saturation for ion channek8
12038 .IPhys. . Chem., Vol. 99, No. 31, I995 In fact, eq 46 should be considered as the general solution for the case that the minimal specific conductivity U E is much larger than the value ac where @(a) equals zero. In this case I$ = 1 and e becomes a linear function of the specific conductivity. In the case of large extemal ion number densities, the following limiting law is obtained:
2.3.4. The Ionic Conductance of Long Channels. If the solution contains multiple ionic species, a limiting drift currentvoltage relation can be derived very similar to the one obtained for the monopolar approximation 29 with
U E and Y E are still functions of j . In the limit of very long channels, however, the minimal specific conductivity converges, UE ac, which makes the result independent of the current density. Also, the more general current-voltage relation as given in eq 30 is obtained with
-
Thus, the long channel conductance as a function of the extemal specific conductivity as shown in Figure 5 holds not only for the monopolar approximation but also in general.
3. Discussion Methods known from theoretical studies of semiconductors to solve the Nemst-Planck equation together with Poisson’s equation (Le. the Nernst-Planck-Vlassov approach) are used and extended in this article to derive integral properties of ion channels, as they can be obtained experimentally from measurements of current-voltage relations at different ion number densities.* Due to the complexity of the theory, this article concentrates only on the theoretical part of the problem. An application to experimental data is in preparation (Syganow, A.; Villarroel, A.; von Kitzing, E. In preparation). In this article, the approach is extended beyond the level commonly used in solid-state physics. The general formulation including multiple ionic species and the derivation of general properties of ion channels without assuming certain model charge distributions are, to the author’s knowledge, new approaches. 3.1. The Calculation of Integral Properties of Ion Channels Using the Nernst-Planck-Vlassov Approach. Two concepts taken from theories used successfully to describe properties of semiconductors (for reviews see refs 45-47) are utilized in this article to derive approximate analytical solutions for integral properties from the Nemst-Planck-Vlassov theory. The phase trajectory approach was originally proposed by Sch0ttky.4~ It transforms the Nemst-Planck equation (1) and Poisson’s equation (2) into two other differential equations (7 and 6) where the electric field E is replaced by the phase trajectory K . The differential equations as a function of the channel coordinate z are reformulated into two differential equations depending on the ion number density n (6 and 7) for a single dominating ionic species or in the case of multiple ions in solution on the specific conductivity a (eq 34).
Syganow and von Kitzing The newly formed equation for the phase trajectory, which cannot be solved in general, consists of three terms. In many situations a single term is considerably smaller than the other two. This idea leads to the concept of regional approximation, where three domains are considered: (a) the equilibrium domain, where the ion number density gradient is compensated for by respectively strong fields and the current density is assumed to be small, (b) the diffusion domain, where the electric field term is neglected and the current density is determined by ion number density gradients, and (c) the drip domain, where the ionic flow is determined by a strong electric field and the contribution due to the density gradient is small. Using the phase trajectory (7), its inverse (6), and multiple integrations by part, the integral properties of the ion channel, such as the total voltage drop V (15 and 16), and the apparent channel length L (18 and 19) are reformulated. Now the integrals no longer depend on the phase trajectory K, and the result contains the phase trajectory only at the boundaries ~i and K ~ save , for a single integral depending on K ~ . The results of the regional approximation are used to derive general integral properties of the channel, such as currentvoltage relations (30) and the channel conductance as a function of external ion number densities (31 and 46). Here, these derivations where done under the assumption of large extemal voltage (IVl > kBT/eo). In this case the influence of the drift domain dominates the integral properties of the channel, but this is not necessarily true for the local properties such as the ionic number density n along the channel. Because the phase trajectory K is small in the drift domain, the integrals depending on K* are neglected. Some correction due to contributions of the diffusion domains are, however, included. A considerable advantage of the approach presented in this article over comparable methods using the Nemst-PlanckVlassov equations to study ion transport through biological channels24.25.27.28.3’-33 is that the phase trajectory approach allows one to derive approximate but analytical formulas for many integral properties of the channels. These results are not exact; however, they provide a guide to the behavior of the channels under many conditions. Numerical integration schemes may be used to test the validity of the formulas in intermediate regions. 3.2. Integral Properties for Biological Ion Channels. The apparent channel length L, calculated from this approach, may be compared with the physical length of known ion channels. The behavior of the channel is determined by the ohmic length LQ (see eqs 26 and 39) and the internal screening length ~ D H (see eq 43). For LQ > L, the GHK approximation becomes valid for the channel. In the case of L > LQ,but L = 21DH, the channel is “short” and its conductance is determined by diffusion, whereas the current-voltage relation for /VI > ~ B T / eo has drift characteristics. For long channels, L >> IDH, all the relations are described by the drift approximation. It should be noted that LQ depends on the external ionic number density no, and changes in this parameter may change the conduction mechanism. There are two characteristic ionic number densities inside the channel, nE and nc. The equilibrium number density nc is related to vanishing charge density in the channel (due to the ions and possible charge contributions from protein side chains) and would be found in an infinitely long channel. The second characteristic ionic density is the minimal number density nE, which results from diffusion of ions from the extemal solution where the ionic number density is assumed to be larger (nE < no and nE < ni). For short channels this minimal number density nE may be considerably larger than the equilibrium number
Ion Transport through Biological Ion Channels
J. Phys. Chem., Vol. 99, No. 31, 1995 12039
a n a x
A
::f 0.6
# 20
40
60
80
100
odaE
Figure 6. Short ion channel conductance as a function of the external specific conductivity a, (see eq 46).
density, nc
