Semi-quantitative model of the gating of KcsA ion channel. 1. Geometry and energetics of the gating

The aim of this series of papers is to develop the semi-quantitative theory of the gating of KcsA channel. For this purpose the available structural and electrophysiological data and the results of molecular dynamics simulations were used in the context of the concept of dynamical self-organization. In the first paper we describe the simplified model of the geometry and energetics of the gating process. This work is the first successful attempt of combining the structure and dynamics of a real protein and the general concept of dynamic self-organization.


Introduction.
The ion channels are classical objects of biophysics, which are studied both experimentally and theoretically for more than fifty years now.The principles of the ion channel functioning were established on the phenomenological level well before the structure of these molecules has become known.Starting from the now famous work of Doyle [1], where high-resolution three-dimensional structures of the KcsA channel were determined, a new era in the ion channel science has started.Rapid advance in the structural and experimental studies on the ion channels allows reexamining some theoretical developments at an entirely new level.One of such theories is the concept of Dynamic Self-Organization (DSO), which was initially formulated before the real structure of the ions channels has become known.The aim of the present series of works is to reanimate the ideas of the DSO in application to the ion channels using the best available data on the channel structure and functioning.Our goal is to show that the DSO is consistent with the structure of the KcsA potassium channel and does not contradict the traditional paradigm of ion channel gating.
Theory and Methods.The molecular dynamics (MD) simulations.The general picture of the opening of KcsA channel is now established.The opening occurs by concerted complex rotation of four transmembrane M2 helices, which leads to the widening of the intracellular part of the channel pore.In addition to global rotation the M2 helices seem to rotate around their long axes and bend near the gating region forming a kink [2][3][4][5].
In the present work the results of MD simulations of the channel opening, described in details in our previous papers [4,5], were used.In these studies the KcsA channel was driven from the equilibrated closed state to the hypothetical open state by means of targeted MD simulations.
The resulted state can only be considered as quasiopen, since the stability of this state and its similarity to the real open state are not known.Thus, we will use the term «quasi-open state» to refer to the results of MD simulations and the term «open state» to refer to the real physiological open state.We believe that the latter state is the same which emerges due to DSO.
First ~2 ns of simulation were used to calculate the channel average structure in the closed state.Last ~0.5 ns of the trajectory were used to determine the averaged structure of the quasi-open state.The calculations of the average structures were performed after aligning all heavy atoms of the selectivity filter for each frame with the corresponding atoms of the reference frame.
The geometry of the model and the structural coordinate.We model the geometry of opening of KcsA channel in a simplified way.We assume that only M2 helices move, while all other parts of the protein remain fixed.Only one of four symmetric subunits is considered.
We model the M2 helix as a rigid rod, which goes approximately through the geometrical axis of the helix.We subdivided the whole helix into 9 segments containing 4 residues each starting from the extracellular side (residue THR85).Remaining 3 residues from the intracellular side were ignored, since they exhibit large fluctuations and cannot be treated as a rigid rod.Each segment corresponds approximately to one turn of the helix, so the geometrical center of each segment provides good approximation of the position of the helix axis for this segment (Fig. 1, A).
The analysis of the MD trajectory of the channel opening shows that each M2 helix rotates approximately in one plane around the well-defined pivot.
Fig. 1, B shows the distances between the centers of segments in the closed and open states.It is clearly seen that the minimal distance along axes Y and Z is observed for segment 4, but the minimal distance along X axis is observed for the segments 2 and 3. Thus, we placed the pivot point between the centers of segments 3 and 4, which seems to be a reasonable approximation.The rod is positioned in such a way that it goes through the pivot and the center of the last 9 th segment of the helix.The centers of the last 9 th segment in the closed and quasi-open states (points T c and T o ) and the pivot point P were used to define the plane, in which the M2 helix moves during the opening.We will reference this plane as an opening plane hereafter (Fig. 2, A).We assume that the rod is confined to this plane.
The extent of opening is described by the angle j (calculated in the opening plane) between the current position of the rod and its position in the closed state.Thus, the angle j is the structural coordinate in our model.The values of the opening angle are j c = 0 o in the closed state and j o » 11.4 o in the quasi-open state.The position of the end point for arbitrary opening angle is defined as T (j).
We assume that the rod can stretch to accommodate small changes in the length of M2 helix upon opening, but can not bend.The length of the rod l for any opening angle can be calculated as linear interpolation between the lengths in the closed and quasi-open states The pore of the channel is aligned with Z axis of the global coordinate system.The center of coordinates (point 0) is placed to the center of the pore in XY plane and to the position of point T c (intracellular end of M2 helix in the closed state) on Z axis.
The permeating ions are assumed to be confined to the channel axis.The one-dimensional pore is assumed to extend from z min = -10 & A to z max = 60 & A. Treatment of the ions.The KcsA channel contains 2-4 K + ions in the physiological conditions. 2 or 3 ions are located in the selectivity filter, while the remaining ion(s) occupies the central chamber of the channel and the intracellular vestibule [1,6,7].It is clear that conventional single-particle diffusion theory is not applicable to the channel pore [8].A simplified model was proposed to describe the concerted motion of strongly interacting ions in the selectivity filter of KcsA channel [8].However, this theory is not applicable to the other parts of the pore.
We are only interested in gating motions, which are much slower than the characteristic time of ionic diffusion in the channel pore.Thus, an interaction of the channel structure with individual ions can be substituted by the interaction with quasi-equilibrium density of charges along the channel pore.
The charge density is in local equilibrium with any current conformation of the channel and changes in a quasi-static way in accord with the conformational changes (adiabatic decoupling of slow degrees of freedom [8][9][10]).The ion-ion interactions can influence the time-averaged charge distribution significantly in the regions, where the ions reside close to each other i. e. in the selectivity filter.It was shown that two ions, which occupy the selectivity filter most of the time, counterbalance deep energy well of the selectivity filter in highly dynamic way by strong ion-ion repulsion.As a result, the third permeating ion initiates the sequence of events known as knock-on barrier-less conduction [8,11].On the longer timescales this process can be reduced to the motion of a single quasi-particle, which transfers the unit charge trough the filter.This means that two ions, which occupy the filter most of the time, can be considered as a part of the «permanent» channel structure and the actual charge distribution of permeating ions is created by the third ion only [8,11].
Thus, in order to compute the quasi-equilibrium charge distribution an approximation of independent ions can be utilized in any part of the pore without significant loss of accuracy and the standard methodology of non-equilibrium Focker-Planck equations [10] where V ion is the energy profile for the permeating ions, c(j, z, C in , C ex , U) is quasi-equilibrium one-dimensional charge distribution in the channel, C in and C ex are effective ionic concentrations of internal and external solutions, U is the transmembrane potential, D is one-dimensional diffusion coefficient of the ions in the pore.This equation is subject to boundary conditions where C in 0 and C ex 0 are real internal and external ionic concentrations, p in and p ex are some coefficients, which convert real dimensional ionic A, which leads to the total number of movable charges in the channel in the range 1-2 for all studied conditions (two «permanent» ions in the selectivity filter are not considered as it is explained above).
Stationary solution of eq. ( 1) is is the concentration ratio.We assume that C in does not change and the only concentration-dependent parameter in our model is r.Expression is used to compute the interaction between the ionic charge distribution and the channel structure as it is described below.
Interaction between the structure and the ions.The charge density in the channel pore interacts with the charges localized in the movable M2 helix.In addition there are short-range Van-der-Waals interactions in the narrow region of the gate, which are responsible for steric closure of the pore in the closed state [2,3,6].We assume that there is only one effective charged group with charge q on the rod, which model the M2 helix.The charged group is positioned at the distance x q l from the intracellular end of the rod.
We also assume that there is only one effective particle, located in the narrowest part of the pore, which interacts with the ions by Van-der-Waals forces (the VDW group).The VDW group is positioned on the «stem» of length l w .The stem is oriented in such a way that it is directed precisely toward the channel axis in the closed state.The stem itself is positioned at the distance x w l from the intracellular end of the rod.
The positions of the effective groups x q and x W serve as adjustable parameters.The value of effective charge q is computed by summing up all atomic charges of the M2 helix as defined in AMBER94 force field [12] used in our MD simulations.It appears that q = 0.99.The arrangement of the effective particles is shown in Fig. 2, A. In order to define the forces, which act on the M2 helix and the permeating ions it is convenient to introduce new moving coordinate system linked to the M2 helix with the orts { r i(j), r j(j), r k(j)}, which depend on the opening angle j.The ort r i(j) is normal to the opening plane and oriented from the channel axis, the ort r j(j) lies in the opening plane and shows the instantaneous direction of motion of the end point of the helix during rotation, the ort r k(j) is directed from the current position of the end of the helix T(j) to the hinge point P (Fig. 2, B).Let us also define the position of the point G, which lie on the rod at the distance l(j)× x from its end (0 < x < 1): We also need to find the distance between the point G and any point z on the channel axis r R(j, z, x) = r G(j , x) -z × r k.This vector determines the direction of the Coulombic force acting between the ion and the charged group.
The projection of this vector on the direction of motion The projection (4) is used below for the qualitative analysis.
The Coulomb force f q , which acts between the charged group and the ion located in point z of the channel axis, is shielded by the water molecules in the pore and various polarizable groups, which lie between the ion and the effective charge.This shielding is distance-dependent and can be very complex.We approximate it with the exponential Debye-like term as it was done in our previous work [8] r r r r ( , , ) ), j x c j x j x j x = -3 where q is the effective charge of the charged group, is dimensional constant, which converts the force to k B T/ & A units, d is the shielding constant.The value of d is an adjustable parameter in our model.
The potential created by this force can be computed analytically, but the expression is too complicated to be given here.The VDW group is positioned at the point We are mainly interested in the strong VDW repulsion, which allows steric closure of the pore in the gating region, thus we neglect weak attractive term and the VDW force f w is written as where w 0 is empirical constant.The potential created by the VDW force is The total force, which acts on the M2 helix from the ion localized at point z, is The momentum of this force projected on the moving ort where m z f z j l q q q q q ( , , ) ( ( , , ) ( )) ( ) ( ) momenta of each force component.Only this projection can rotate the helix, while two other projections are counterbalanced by the reaction forces (the rod is confined to the plane) and thus neglected.One can compute the average momentum created by the quasi-equilibrium distribution of the ions in the channel given by eq. ( 2)

M
r U q w ( , , , , ) This equation ( 7) is only valid in the adiabatic approximation, which was used to calculate the distribution c.According to this approximation the influence of individual ions is averaged out and the M2 helix only «feels» the average quasi-equilibrium charge density in the channel pore.
Modeling the rest of the channel -an effective energy profile for the permeating ions.Since only the M2 helices are modeled explicitly, the rest of the channel is taken into account by introducing the effective energy profile (EEP) for the permeating ions created by all remaining parts of the channel structure.This profile includes implicitly not only the interactions with the channel walls, but also with other ions and water molecules.We presume that EEP in the selectivity filter does not change significantly upon opening and remains flat, thus all subsequent considerations are applied to the chamber, gate and intracellular part of the pore.
Let us express the effective energy profile for the ion V ion as V ion = V M2 + V EEP + U, where V M2 is the profile created by the M2 helices, V EEP is the profile created by the rest of the channel, water molecules and the ion-ion interactions, U is the transmembrane electrochemical potential.At this stage we assume that transmembrane potential is zero U = 0.The V M2 term is determined by the geometry of our model and the positions of charged and VDW groups.Using ( 5) and ( 6) we can write where the coefficient 4 accounts for existence of four identical subunits in KcsA channel.
In order to estimate V EEP one can use the following considerations.In the open state the KcsA channel conducts the ions at the rate, which is close to the diffusion limit, which means that there are no large energy barriers or deep wells on the effective energy profile.This means that the forces, which act on the permeating ion from the M2 helices, are counterbalanced almost completely by other forces, which contribute to V EEP : It is not known whether such compensation occurs in the closed state, because the channel conducts no current in this state.However, it seems to be established that the total number of ions accommodated by the channel depends on the opening state only slightly.This means that no additional deep energy wells appear in the closed state of the channel.Thus, one can assume that the rest of the channel compensates the potential profile created by the M2 helix almost completely at any stage of the opening and V EEP can be considered independent on the opening angle.
The exact shape of the effective energy profile in the presence of compensation is still unknown, but the depths of any energy wells or heights of any barriers can not exceed 1-2 k B T. It is usually believed that the chamber of the channel contains substantial density of the ions most of the time, in particular, an additional ion is usually placed into the chamber in most of MD simulations [2,6,13].The effective energy profile should have a shallow energy well in the chamber region to account for this effect.We modeled this well by an inverted Gaussian curve centered in the chamber region where s is the half-width of the energy well, A is the depth of the well, z c is the position of well.This well is assumed to be independent on the opening state as it is described above.
Structural potential for the M2 helix.Slow diffusive rotational motion of the M2 helix can be described by the stochastic Langevin equation.Using (7) where D j is diffusion coefficient of rotational motion, x(t) is the white noise, V str is so-called structural potential.Structural potential is the potential of mean force, which acts on the moving helix from the other parts of the channel but not from the ions in the pore.The computation of the potential of mean force V str for such complex moving object as the long M2 helix is not possible using standard MD protocols, thus we used simplified semi-quantitative approach to estimate V str .The rotation of the M2 helices is the dominant motion in the system during opening, thus it is possible to assume that this rotation makes the major contribution to the change of total energy E total .If this is true, one can write as a first approximation This allows to monitor the total energy of the system as a function of the opening angle instead of computing the potential of mean force explicitly, but resulting V str can only be considered as semi-quantitative approximation.
We analyzed the trajectories of opening by plotting the total energy of the system against the opening angle (Fig. 3).We approximated it as a sum of two components.The first one, V wall describes sharp «wall» at the right and similar wall at the left, which is not sampled in MD due to the direction of targeted pulling.The second linear component, V lin approximates the bottom of the potential.The final expression for V str can be written using as ) .
), j where coefficient b is identical to the corresponding coefficient in (10) and incorporates the conversion factor from dimensional MD units to dimensionless model units.We use the value b = 0017 .
k T B in all subsequent calculations.
Results and Discussion.In the present work we have established a simplified framework describing the geometry and the energetics of the KcsA channel gating based on the concept of dynamical selforganization.We used the geometry of the channel opening obtained from experiments and MD simulations.The structural potential for the motion of the M2 helix is derived in the semi-qualitative way from the MD results for single forced opening of the channel.Once the opening pathway of KcsA is determined with sufficient degree of confidence in atomic details, the all-atom model can be built easily using the principles described in this work.In the next paper the obtained data will be used to develop a closed DSO-based theory of the channel gating.
Acknowledgements.V. Kh. thanks the University of French-Comte for providing the invited professor grant, which stimulated this work.Prof. L. N. Christophorov is acknowledged for critical reading of the manuscript.

Fig. 1 .
Fig. 1.Locating the pivot point on M2 helix: A -Centers of the 4-residue segments in the closed (silver balls) and quasi-open (black balls) states superimposed into the cartoon representation of the single subunit of KcsA channel.The rod, which models the M2 helix, is shown as a solid line.The indexes of the segments are indicated; B -The distance between the centers of the segments in the closed and the quasi-open states.The connecting lines are drawn to aid the comparison only

Fig. 3 .
Fig. 3. Structural potential for the M2 helix obtained in MD simulations and the fit used in the model.Each dot corresponds to a single frame of the targeted MD trajectory.See text for details