Discrete-State Representation of Ion Permeation Coupled to Fast

ARTICLE pubs.acs.org/JPCA

Discrete-State Representation of Ion Permeation Coupled to Fast Gating in a Model of CLC-Chloride Channels: Analytic Estimation of the State-to-State Rate Constants Rob D. Coalson* and Mary Hongying Cheng Department of Chemistry, University of Pittsburgh, Pittsburgh, Pennsylvania 15260, United States ABSTRACT: Analytical estimation of state-to-state rate constants is carried out for a recently developed discrete state model of chloride ion motion in a CLC chloride channel (Coalson and Cheng, J. Phys. Chem. B 2010, 114, 1424). In the original presentation of this model, the same rate constants were evaluated via three-dimensional Brownian dynamics simulations. The underlying dynamical theory is an appropriate single- or multiparticle three-dimensional Smoluchowski equation. Taking advantage of approximate geometric symmetries (based on the details of the model channel geometry), well-known formulas for state-to-state transition rates are appealed to herein and adapted as necessary to the problem at hand. Rates of ionic influx from a bulk electrolyte reservoir to the nearest binding site within the channel pore are particularly challenging to compute analytically because they reflect multi-ion interactions (as opposed to single-ion dynamics). A simple empirical correction factor is added to the singleion rate constant formula in this case to account for the saturation of influx rate constants with increasing bulk Cl
concentration. Overall, the agreement between all analytically estimated rate constants is within a factor of 2 of those computed via threedimensional Brownian dynamics simulations, and often better than this. Current
concentration curves obtained using rate constants derived from these two different computational approaches agree to within 25%.

1. INTRODUCTION Recent advances in structural and physiological characterization of the CLC superfamily have provided fascinating clues about the evolutionary development and the functional mechanisms of these proteins.1
3 X-ray crystallization of the bacterial members of this family revealed a double-barreled pore with a narrow selectivity filter of 15 Å in the center of the protein.4,5 The crystallized wild-type EcClC or StClC4,5 displays two anionbinding sites: one near the intracellular entrance to the selectivity filter (Sint), and the other in the center of the filter (Scen). A third site (Sext), near the extracellular entrance to the filter, is occupied by the gating side chain of a glutamate residue (E148). In bacterial CLC mutants such as E148A or E148Q,5 all three binding sites can be occupied simultaneously by anions.6 Subsequent electrophysiological studies revealed that several CLC proteins (including the bacterial EcClC) are actually exchangers, moving 2 Cl
ions from the intracellular to the extracellular region while one Hþ is shuttled from the extracellular to the intracellular space.7
9 Other CLC family members behave as standard passive ion channels, allowing Cl
ions to flow continually down an electrochemical gradient established by a chloride concentration imbalance in intra- vs extracellular reservoirs or by a voltage differential across the membrane (as can be achieved via microelectrode techniques in a laboratory).1 In light of the exchanger characteristics of other family members, it is reasonable to suppose that those members of the CLC family which behave like passive ion channels are evolutionarily degraded variants of the original exchanger members of the same family.2 r 2011 American Chemical Society

All of these features suggest that the coupling between channel gating and ion permeation should be particularly intricate in the CLC chloride family. Indeed, CLC channels are known to be gated by changes in pH in the extracellular reservoir,5,10 changes in reservoir Cl
concentration,1,11
13 changes in membrane potential,10,14,15 and, for some members of this family, ATP binding to an cytoplasmic domain of the channels.16
18 This presents great challenges for theory and modeling, since the time scales of the relevant processes are long, and many aspects of gating involve significant conformational rearrangements of the protein. Thus, attempts to utilize all-atom molecular dynamics (MD) to directly simulate channel gating dynamics on a time scale long enough to reveal its complete mechanism are unlikely to succeed. Instead, coarse-grained techniques and models are called for. These, of course, have their own challenges, among them being to demonstrate the relevance of the coarse-grained model to allatom dynamics (as reflected by experiments). In an effort to address these issues, we developed a simplistic but fully three-dimensional (3D) Brownian dynamics (BD) model of Cl
ion permeation through a CLC-type channel coupled to opening and closing of the fast gate.19 Our BD model was based on structural and electrophysiological studies by Dutzler et al.,5 which Special Issue: David W. Pratt Festschrift Received: January 24, 2011 Revised: April 29, 2011 Published: June 21, 2011 9633

dx.doi.org/10.1021/jp200749s | J. Phys. Chem. A 2011, 115, 9633–9642

The Journal of Physical Chemistry A

ARTICLE

Figure 1. Panel A depicts the three-dimensional geometry of the multi-ion continuous space BD simulation model. Pink disks indicate the three Cl
binding sites (from left to right) Sint, Scen, and Sext. The fast gate is represented as a ball on a pivot having two stable states: one pinned back away from the pore and the other occupying the Sext binding site. Panel B shows the 12 possible states that can arise in a discrete state model based on three single occupancy binding sites. A blue disk indicates a Cl
ion in the site, an orange square indicates that the gate particle occupies the Sext binding site, and an unfilled circle indicates an empty site.

suggest the E148 residue in bacterial CLC family members (E166 in the eukaryotic analogue CLC-0) functions as the fast gate by swinging between two stable states, one of which competes with Cl
for the Sext binding site within the channel pore, and the other of which is “pinned back” from the pore into the extracellular space (thus allowing Cl
ions to pass through the pore). Although the recent eukaryotic CmCLC transporter structure from MacKinnon’s group20 indicates that in CmCLC the gating glutamate side chain occupies the Scen site instead of Sext, this does not change the essence of the modeling problem or the underlying dynamical mechanisms of the ion transport. The molecular conformational change involved here was represented in our computational model via a spherical particle attached to an anchored pivot rod and allowed to move along one angular coordinate.19 A Brownian dynamics simulation involving many mobile Naþ and Cl
ions plus one spherical gate particle attached to a pivot rod (described in detail below) strains currently extant computational resources. In fact, the simulations of gate closing presented in ref 19 could only be carried out using an artificially low activation barrier for the conformational transition between open and closed gate particle configurations. These time scales and other complexities involved do not bode well for generalizing a multiparticle continuous space model to treat other important features of the permeation/gating process, including the role and the mechanism of protonation of the fast gate, additional conformational changes introduced by an applied trans-membrane potential which presumably underlie the voltage dependence of gating in CLC channels, etc. As a complementary approach, discrete state models have long found productive use in analyzing ion channel function, e.g., ion permeation through the open state conformation of the channel, channel gating, etc.21
26 A number of interesting applications of such models to structure/function relations in the CLC family of channels and antiporters have been presented in recent papers.20,27
29 In a recent paper28 we pursued the possibility of reducing our multiparticle continuous space 3D BD model to one which

involves a small finite number of states, based on single occupancy restrictions for each of the three Cl
binding sites in the system. These states are connected by pairwise rate constants. Once numerical values have been assigned to all relevant rate constants, linear kinetic master equations are then solved to time-evolve the entire set of system state occupation probabilities, thus enabling analysis of all important dynamical processes, including ion permeation through the open pore configuration of the channel. The time-averaged rates of gate opening and closing can also be computed, simply and regardless of time scales involved, using this approach. In ref 28 we computed rate constants for each state-to-state transition in our discrete state model by performing BD simulations focused on the relevant transition. We found, encouragingly, that a kinetic master equation treatment of the dynamics of the resultant discrete state model gave good agreement with the properties of the full multiparticle continuous space 3D BD model, including the flux of Cl
ions through the open channel and the average rate of gate closing. In this approach we must still rely on time-consuming stochastic simulations for the numerical values of each rate constant. It would be useful to have a more analytic route to evaluating these rate constants. Not only would this provide insights into the relation between input parameters such as the single ion potential of mean force and diffusivity profile governing ion or gate particle motion and the corresponding rate constants but it would allow rapid scanning over a range of numerical values for these parameters. Such an analysis is the focus of the present work. In section 2 of this paper we briefly review relevant aspects of the discrete state model of Cl
permeation coupled to fast-gate opening/closing that was introduced in ref 28 and the basic equations underlying the kinetics of this model. In section 3, we analyze state-to-state transition rates associated with site-to-site hopping within the channel, using elementary Kramers theory.30 In section 4, we adapt elementary solutions of the single-particle Smoluchowski equation for a particle moving in a radially 9634

dx.doi.org/10.1021/jp200749s |J. Phys. Chem. A 2011, 115, 9633–9642

The Journal of Physical Chemistry A

ARTICLE

symmetric potential and diffusivity profile to predict rates of Cl
influx from the vestibule of our model CLC channel into binding sites located just inside the channel pore, and the analogous efflux rates from such binding sites into the abutting ionic reservoir. In addition, we suggest a simple empirical correction factor to account for saturation of the influx rates with increasing bathing solution [Cl
]. Finally, conclusions are provided in section 5.

2. A 3D MULTI-ION MODEL OF PERMEATION AND FAST GATE MOTION IN CLC CHANNELS AND ITS REDUCTION TO AN EIGHT DISCRETE STATE MODEL In ref 19, we constructed an approximate 3D model of a CLC channel permeation pore, based on X-ray structures of bacterial CLC family members4,5 and simple energetic considerations. Dynamics of ion permeation was then modeled at the Brownian dynamics level,31,32 treating sodium and chloride ions as mobile particles, water as a dielectric continuum, and the protein
membrane system as frozen in a single static configuration. In the CLC superfamily, the E148 residue of bacterial StCLC (E166 in CLC-0 and other eukaryotic homologues) is thought to act as a fast gate,2,5,11,13 by undergoing a conformational change which brings it into the channel, thus competing with a Cl
ion for one of the energetic binding sites in the channel (cf. Figure 1). In the mechanical model developed in ref 19, the glutamate gate was accounted for as a spherical particle attached via a rigid rod to a pivot point located in the channel, such that the gate particle can swing from the Sext binding site in the channel out into the external reservoir (where it does not occlude the pore). Energetically, this gate particle was subjected to a bistable potential (dependent on the angle of the pivot rod), allowing for one stable configuration in which the gate particle occupies the Sext binding site and a second stable configuration in which it is pinned back away from the pore. These minima in the single-particle effective potential seen by the gate particle are separated by an energy barrier, so that the gate particle undergoes activated transitions between the two stable states (potential energy minima). Ultimately, the gate particle was included in BD simulations of Naþ and Cl
ions, coupling to the motion of these ions through excluded volume interactions (i.e., spatial overlap between a Cl
ion and the gate particle is forbidden). In ref 28, we sought to reduce this complicated interacting many-body dynamical system to one involving a small number of discrete states connected pairwise via state-to-state rate constants. These states derive from the fact that there are three deep binding sites in the CLC well, namely, Sint, Scen, and Sext (cf. Figure 1). Our discrete state model includes states in which each of these three binding sites is either unoccupied or singly occupied by a Cl
ion. In addition, the Sext site can be occupied by the gate particle. None of the sites can be multiply occupied (e.g., via two Cl
ions, or in the case of Sext, by one Cl
ion and the gate particle). This leads to a 12 state model, as indicated in Figure 1. Given appropriate rate constants that describe transitions between these states, the time evolution of the state probabilities was assumed to be governed by a set of linear kinetic master equations. Specifically, let pi(t) denote the probability to find the ion channel system in state i. Then the time evolution of the state occupation probabilities of the N states in the system is given by the following set of linear ordinary differential equations dpi ðtÞ=dt ¼

N X j¼1

½kji pj ðtÞ
kij pi ðtÞ

;

i ¼ 1
N

ð1Þ

Figure 2. Diagram of effective free volume channel geometry. The free volume of the channel (bounded by dotted line) is obtained by rolling a sphere of Cl
ion radius 1.8 Å around the interior surface of the channel (shown as solid black line). Three filled disks indicate the three Cl
binding sites. An axis with tick marks at the bottom orients the channel axis z. The vestibule on the left connects to the intracellular reservoir, while the vestibule on the right connects to the extracellular vestibule.

where kij is the probability per unit time for a system in state i to make a transition to state j, with kii = 0, i = 1
N by construction. As noted above, there are formally 12 states in the discrete state model. However, the specific processes that were simulated in ref 19 could be modeled via a discrete state model consisting of only the eight “open channel” states (cf. Figure 1 and ref 22). Thus in ref 28, we considered transitions only within this eight-state manifold, as detailed via eq 1 with N = 8. In ref 28 calculation of these rate constants was carried out via 3D BD simulation. With numerical values of all nonzero rate constants in the model in hand, appeal to eq 1 gave satisfactory agreement with the results for both (i) Cl
ion flux through the open channel and (ii) gate closing rates which had been obtained in ref 19 via multiparticle continuous space 3D BD simulations. In the present paper we revisit the calculation of the rate constants for ion permeation when the channel is open and attempt to estimate them by appealing to underlying microscopic kinetic laws. In particular, since in our model ions execute motion in accord with the Smoluchowski equation (equivalent to the phase-space Langevin equation in the high-friction limit30), we solve the Smoluchowski equation for the relevant geometrical configurations, single ion energetics, and boundary conditions, as discussed in detail below. Before proceeding to the rate constant analysis, it behooves us to consider the effective channel dimensions which govern the motion of a Cl
ion through the channel in the independent ion limit (neglecting ion
ion repulsions). Each ion then behaves as a point particle (localized at the physical ion’s center), moving in a “free” volume which is restricted due to the finite size of the Cl
ion. Since the Cl
ion has a radius of R = 1.8 Å, the free volume of the cone and cylinder regions of the channel is significantly reduced compared to the full pore/vestibule volume of the system. A reasonable way to affix the free volume of the channel is to roll a sphere of radius 1.8 Å around the interior surface of the structure. The effective boundary surface of the channel is then defined by the location of the sphere center: the center of an ion 9635

dx.doi.org/10.1021/jp200749s |J. Phys. Chem. A 2011, 115, 9633–9642

The Journal of Physical Chemistry A

ARTICLE

Table 1. Table of Site-to-Site Hopping Rate Constants (in ns
1), Comparing Formula vs BD Simulation

Figure 3. “Bare” single Cl
ion PMF (no ions in the channel) is shown via dashed blue line. An ion pinned at the binding site Sint generates a modified single Cl
ion PMF (red line) in the Scen
Sext region of the channel.

in the physical channel cannot penetrate further than the surface defined in the manner just described. The resultant effective channel structure is depicted in Figure 2 This is the relevant geometry for single ion Smoluchowski dynamics that underlies analytic estimation of the state-to-state rate constants relevant to our discrete state model. [Note: In the present work we consider only the process of ion transport through the open pore structure of our model CLC channel. In this configuration the gate particle is pinned back away from the pore, and charge neutral (corresponding to a protonated glutamate residue, as discussed at length in ref 28). Thus its presence in the extracellular vestibule has negligible effect on ion permeation dynamics. Hence we remove it from our analytical model, which then becomes fully symmetric, geometrically, in the intracellular and extracellular vestibules: compare Figures 1 and 2.]

3. COMPUTATION OF INTERNAL (SITE-TO-SITE) HOPPING RATE CONSTANTS Here we can directly apply the Kramers formula for the rate constant of an activated process taking place in 1D in the highfriction (Smoluchowski equation) limit. The ion is actually moving in three dimensions, in the cylindrical pore region of the channel. However, since the radius of the channel is taken in our model to be constant throughout the pore, the geometrical details in the transverse directions are irrelevant: they cancel out of the final result for the rate constant, which, again, is given by the well-known Kramers formula19,30 pffiffiffiffiffiffiffiffiffiffi kw kb ð2Þ expð
βEa Þ k ¼ βD 2π Here kw is the second derivative (curvature) of the single-ion potential of mean force (PMF) evaluated at the local minimum corresponding to the binding site, while kb is the negative of the second derivative of the PMF evaluated at the local maximum; furthermore, Ea is the activation free energy for the process (the free energy of the barrier top relative to that at the bottom of the binding site well), D is the appropriate particle diffusivity, and β = 1/kBT. As an example, let us consider the rate constant k34, which describes hopping from Scen to Sext. If the channel is otherwise empty (no Cl
in Sint), then the relevant single ion PMF is shown in Figure 3 (dashed blue line). This is the same functional form for the single ion PMF that was employed in our previous 3D multi-ion BD simulations of ion permeation through the channel.19,28 (Note: This PMF includes a contribution due to

rate constants

3D BD

formula

k23

0.555

0.581

k34 k56

0.524 0.983

0.575 1.063

k67

0.310

0.283

k32

0.041

0.056

k43

0.694

0.896

k65

0.397

0.468

k76

0.136

0.147

an applied transmembrane potential of Vap =
100 mV, as was assumed in all the 3D BD simulations performed in ref 28. The same value of Vap will be assumed throughout the present work.) To calculate the hopping rate from Scen to Sext, we appeal to eq 2. In particular, we compute the curvature at Scen (kw) and the curvature at the maximum barrier height separating the Scen and Sext wells (this is
kb). Next, we need to assess the activation energy (PMF at the barrier maximum minus PMF at the position of the Scen well minimum). Finally, as in our previous studies,19,28 the relevant ion diffusion constant for Cl
site
site hopping processes within the channel is Dint (Dint = 20 Å2/ns). (Since the CLC pore is very narrow and Cl
permeation is coordinated not only by side chains but also by backbone atoms,33 Dint is expected to be significantly reduced relative to its bulk solution value of Dbulk = 200 Å2/ns, similar to previous findings for Kþ diffusivity inside the narrow Gramicidin A channel.34) Plugging these parameters into eq 2 gives k34 = 0.575/ns, in good agreement with the value of 0.524/ns obtained from 3D BD simulation. To calculate the reverse hopping rate constant, k43, we simply exchange the roles of initial and final states: the curvature at Sext is used to determine kw, and the relevant activation energy is given by the PMF at the barrier maximum minus the PMF at the Sext well minimum. The result of this analytical estimate also compares reasonably well to the value obtained via 3D BD simulations, as indicated in Table 1. Now, if there is a Cl
ion in the “third” site, it generates Coulombic interactions with the hopping ion that modify the relevant single ion PMF experienced by the hopping ion. For example, in the case of hopping from Scen to Sext (or vice versa), when a Cl
ion occupies Sint the single ion PMF of the hopping ion is altered due to the additional Coulombic repulsion generated by the ion fixed in Sint. Calculating this ion
ion interaction using a previously developed continuum electrostatics approximation35 that takes into account the effects of charge induced on the surface of the channel pore
interior water interface (due to dielectric inhomogeneity in the medium), we obtain the overall PMF curve for the motion of the ion hopping between Scen and Sext shown via the solid red line in Figure 3. Clearly, since the electrostatic force between two ions inside the channel is large, the PMF can be modified significantly by having another ion in the channel. The Kramers formula eq 2 provides direct insight into how these PMF modifications will affect the forward and backward hopping rates. For the case shown in Figure 3, the Sext well is now lower in energy than the Scen well. As a result, the forward rate constant, i.e., k56, is larger than the backward rate constant, i.e., k65, reversing the behavior found when there is no ion in the Sint site. Numerical values obtained via the Kramers formula, eq 2, are in good agreement with 3D BD 9636

dx.doi.org/10.1021/jp200749s |J. Phys. Chem. A 2011, 115, 9633–9642

The Journal of Physical Chemistry A

Figure 4. Expanded view of half of the model ion channel. In particular, the green lines indicate the side boundary of the free volume cone for a single Cl
ion moving in this region of the channel. The purple line marked Ro indicates a hypothetical cap of this cone in the intracellular region of the channel, while the purple line marked Ri denotes a hypothetical inner capping surface corresponding to the entrance to the channel pore (cylindrical section of the model channel).

simulation results (