From the potential of the mean force to the quasiparticle ’ s effective potential in narrow ion channels

We consider the selective permeation of ions through narrow water-filled channels in the presence of strong interaction between the ions. These interactions lead to highly correlated ionic motion, which can conveniently be described via the concept of a quasiparticle. Here, we connect the quasiparticle’s effective potential and the multi-ion potential of the mean force, found through molecular dynamics simulations, and we validate the method on an analytical toy model of the KcsA channel. Possible future applications of the method to the connection between molecular dynamical calculations and the experimentally measured current-voltage and current-concentration characteristics of the channel are discussed.


Introduction
The selective permeation of ions through narrow water-filled pores is an unsolved problem that continues to attract great attention.In Nature, this permeation happens through ion channels which constitute important control elements in biology [1,2].Over the last several decades their properties have become a focus of intensive research, not least on account of their role as potential pharmacological targets.That is why efficient theoretical and numerical methods of coupling a channel's structure to its properties -e.g conductivity, selectivity, and rectification -still remain a subject of active research.Although there has been significant progress in that direction, many unsolved problems [1,3] still remain.
The main challenge is the presence of strong all-to-all, not necessarily classical, interactions at short and long ranges between the relevant entities [1].The latter include dissolved ions, water molecules, and the protein backbone.The finite size and discrete charge of the ions need to be taken into account.So also do the waterion (hydration), water-protein and water-water interactions, as well as possible interactions between the ions and the lipid bilayer.Moreover, a protein's structure is flexible rather than rigid, and the molecule polarizes when a charged object approaches it closely.Molecular dynamics [4,5] encompasses all of these interactions, but at the cost of substantial computational time.One of the most important outputs of molecular dynamics (MD) is the potential of the mean force (PMF) [6].Brownian dynamics (BD) [7,8] provides a more computationally efficient approach to the problem.It describes the motion of individual ions with ion-water collisions being included implicitly as fluctuations.The strengths of the two approaches can be unified by calculating the forces between particles in BD simulations through use of the PMFs derived in MD simulations.This tandem approach has been applied successfully to a number of nanoscale systems [6,[9][10][11].Importantly, in the narrow channels with strong interaction between ions MD results in multi-ion PMFs that differ according to the number of ions in the channel.
The corresponding continuous description has evolved along the lines of Poisson-Nernst-Planck (PNP) [12].The rigorous derivation of the PNP relations for a given system requires ad hoc closing relations which have not been verified in narrow channels [13].Importantly, the discreteness of the ions and interaction between them are missed, resulting in disagreement with experimental and BD studies [14].The Fokker-Planck (FP) method faces the difficulty of the choice of correct boundary conditions [15].Thus, although being in some ways more attractive, the continuous description cannot be applied directly to strongly-interacting ions in a narrow channel.Possible solutions of the problem can utilize e.g.hierarchical FP systems [16] or the concept of quasiparticles.
The notion of a quasiparticle was introduced to describe collective ionic motion in narrow channels [17,18].This concept reduces the complexity of many-body motion to the one-dimensional motion of a single quasiparticle in an effective potential, providing a very clear physical interpretation of permeation [16,17].The idea was introduced in relation to a simple toy model, and the quasiparticle concept still needs to be connected to the structures of real channels.
Here, we extend this concept and derive an explicit relationship between the effective potential of a quasiparticle and the multi-ion PMFs related to the structure of the channel.It is validated on the toy model of ion-ion and ion-channel interactions from Refs.[17,18] to simplify the comparison between the theoretical and BD simulation results.Our approach is equally applicable to the more complex real biological structures, providing a clear route from a protein's structural properties to a quasiparticle's effective potential.

The motion of individual ions
The Streptomyces lividans KcsA is a prokaryotic pH-activated K + channel.Its highly-selective pore region -the selectivity filter (SF, left panel in Fig. 1) -contains the amino acid sequence TVGYG, which is highly conserved among both The quasiparticle's effective potential from the PMF 3 prokaryotic and eukaryotic K + voltage channels.The latter fact, the discovered KcsA crystallographic structure [19], and its ability to discriminate between K + and Na + with accuracy 1000:1 at the diffusion rate ∼ 10 8 ions/sec have made this channel a reference specimen for many studies of ion channels.
The SF of KcsA is only a few Å wide [19].Due to this narrowness, the permeating ions cannot pass each other and therefore move effectively in one dimension.This allows us to simplify ionic motion in the SF to dynamics along the channel's longitudinal axis z.The Langevin equation of motion for the k th ion of species α is [20] where m α is ion's mass, γ α is its friction coefficient which is coupled to its diffusivity D α by the Einstein relation ) represents the force acting on an ion located at r α k when N ions are present in the channel, k B is Boltzmann's constant, T represents the absolute temperature, and ξ is three-dimensional white noise of unit intensity ( ξ = 0, ξ(t + τ )ξ(t) = δ(τ )).
The force is given by where an N -particle PMF W N (r α 1 , r α 2 , . . . ) has been introduced [5, 6, 9, 21] Here C is a constant, and U all (r α 1 , . . ., r α N , X) represents the potential energy due to interactions between all atoms, i.e. permeable ions located at {r α 1 , . . ., r α N } and other atoms at X = {r β 1 , r β 2 , . . .}.The PMF in the selectivity filter of the KcsA channel is illustrated in the right panel of Fig. 1.
The use of MD-generated PMFs in studies of nanoscale systems yields several benefits.For instance, it self-consistently includes the ion's interactions with other ions in the pore, the effects of dehydration at the entrance and inside the channel, the induced charges on the pore walls, and the influence of the flexible pore structure.Using state-dependent PMFs, we take a step towards using the dynamical potential landscape [10].

Motion of quasiparticles
Importantly, the one-dimensional motion and strong ion-ion interactions allow one to describe ionic collective [19] behaviour as the motion of a "quasiparticle" (QP) [17], represented by the centre of mass q α and relative distances p α m between neighbouring ions of a given species with respective velocities v α = żα and ṗα k .Here ions are numbered from outside to inside as in [5].The inverse coordinate transformation z α k = z α k (q N , p) is given by Cramer's rule where, for clarity of notation, the set of relative distances is written as a vector [17] p = {p 1 , p 2 , . . ., p N −1 }.Generalization for multiple species is straightforward, and unless explicitly stated we omit species' indices in specifying the ion's and QP's coordinates, for clarity.It is worth noting that the QP defined above [17,18], also known as a "quasi-ion" or a "super-ion" [16], represents the centre of mass of the ions, whereas a "quasi-ion" ("permion") [12,22] includes the ion, water molecules and the protein channel.
The potential energy of the system of ions located at Z = (z 1 , z 2 , . . ., z N ) is given by the PMF W N (z 1 , z 2 , . . ., z N ).By expressing the coordinates of individual ions z k = z k (q, p) via their centre-of-mass and relative distances {q, p}, one can express the energy of the system via the parameters of the QP: W N = W N (q, p).
In the spirit of Ref. [17], we use the full Langevin equation (1) to describe the The quasiparticle's effective potential from the PMF 5 evolution of the QP's position q α in time The coefficient D * N represents the transport diffusivity [23].For interacting particles, this coefficient differs from a simple product N D between the number of ions N and diffusivity D given in Ref. [17].
Applying the chain rule and making use of the definitions ( 3)-( 4), one arrives at The above equations apply to the diffusion of a QP consisting of a fixed number of ions.However, during permeation this number varies.For instance, the knock-on mechanism, when an incoming ion causes the furthermost one to leave the channel, appears to be inherent to the KcsA potassium channel [19].The non-constant N has two important consequences.First, it implies that the set of energy landscapes W N for a QP is discrete, as we demonstrate explicitly in section 3. Secondly, the discrete changes of pore occupancy lead to spatial jumps of the QP (Fig. 2(a)).The positions of a QP before and after a jump are coupled [17].For instance, for an N -ion QP located at q N , the entry (exit) of an ion from side S relocates it to q N +1 (to q N −1 ) according to [17] ion enters: as shown in Fig. 2(b).Thirdly, the diffusivity, charge, mass, and effective potential of the quasiparticle simultaneously change their values in a jump-like manner as well.Fig. 2(a) visualizes motion during a typical simulation.Following the trajectory of single ions in the channel (black traces), one observes the centre-of-mass dynamics according to Eq. ( 7) (red).Ionic diffusion results in the corresponding diffusion of the quasiparticle.At approximately 0.5ns after the initial time, one of three ions escapes from the channel which results in the QP abruptly jumping to a new position.Diffusion of the ion at the boundary back and forth for some time gives rise to a series of jumps of the QP.These jumps cease once the ion finally leaves the neighbourhood of the channel and travels further into one of reservoirs.
The correlation coefficients between coordinates are found to be within [0.75, 0, 85].These high values confirm that the strong ion-ion interaction results in high correlation of the individual ions and their centre of mass, the quasiparticle.Individual ionic trajectories (black and grey traces) result in the motion of their corresponding QP (red).At time t ∼ 0.5ns the rightmost ion leaves the channel; two ions remain, and the QP consequently jumps from q 3 to a new position q 2 .(b) The connection between the QP coordinates before (plotted on the abscissa) and after (on the ordinate) a jump.Each BD data point describes a single jump q 3 to q 2 (a circle) or q 2 to q 3 (a cross).The set of q 3 to q 2 points fall on the steeper pair of parallel lines, each point corresponding to an individual ion's entry/exit through either the left (dashed lines) or the right (solid lines) edge of the channel, as illustrated in the inset where individual ions are grey and the larger orange circle represents the QP.The set of q 2 to q 3 points fall on the similar but shallower pair of lines as described.In each case, the lines represent the theoretical predictions of Eq. ( 8).The red square indicates the q 3 to q 2 jump shown in panel (a), when the rightmost ion exits the channel (pale grey trajectory).8) for the QP's coordinates immediately before and after a jump.Depending on the side of the entry/exit, the positions of the QP before and after form two parallel lines.The nonuniform distribution of original and final positions is due to the presence of interactions, and their centroids define which transition occurs.For instance, a 3-ion quasiparticle resides in two more probable positions (maxima of distributions in (c-d)).This means that the three ions are located closer to the inner side of the channel, and one would expect the innermost ion to leave the channel.This does happen (transition 3 → 2, red square), with the QP relocating to the centre.
We also note that in the KcsA selectivity filter K + ions reside at the binding sites formed by carbonyl oxygen atoms [19].Strong interactions with the latter define the localization of individual ions very precisely [4,5,9].The exact localization of individual ions results in a much sharper, almost discrete, distribution of quasiparticles during permeation.
The quasiparticle's effective potential from the PMF 7 2.3.Coupling PMF W N (r 1 , r 2 , . . . ) to the effective potential U ef f Equation ( 7) describes the dynamics of the position q of the quasiparticle, but still contains ion-ion distances p.To simplify Eqs.(7) further, the underlying properties of interactions in the channel must be discussed.Due to the strong interaction between ions, the relative distances p m reach their equilibrium values rapidly, while the centre of mass moves adiabatically [17].This allows one to introduce the equilibrium distribution of mutual distances p, in terms of the position q N of the QP, Utilizing the equilibrium property ( 9), we multiply both sides of Eq. ( 7) by the distribution function ( 9) and integrate over all {p k }.We also rewrite Eq. ( 5) in overdamped form.This yields the Langevin equation for the expected position of the quasiparticle q with In the integral above one has to ensure that the coordinates of individual ions x k (q, p) lie inside the channel.Equation ( 11) is the main result of this paper.It explicitly couples the PMF W N and the effective potential U ef f N for a quasiparticle whose dynamics evolves according to Eq. (10).Thus, the notion of quasiparticles reduces the many-body problem (1) to effective one-body motion.Moreover, the coupling with MD by means of a PMF allows one to introduce atomistic details of the channel structure.The corresponding probabilistic description, given by the coupled differential Chapman-Kolmogorov equations, will be published elsewhere.

Simplifying assumptions for BD modelling
In order to verify the applicability of Eq. ( 11), we run Brownian dynamics simulations in the toy model of ion-ion and ion-channel interactions proposed in Ref. [17].This model envisages N ions interacting with each other and with the channel, the energy PMF of the system being given by The energy comprises the ion-channel interaction and the ion-ion potential where the first term represents the screened Coulomb interaction [17] with shielding constant d and dielectric permittivity , and ionic charge c, while last term includes short-range repulsion between ions at small distance [24].The parameters take the following values: a = 9 Å, d = 2.8 Å, U 0 = 10.5kB T , F 0 = 2 × 10 −10 N, r 0 = 2.8 Å, = 1 inside the channel and = 80 in the bulk, D = 2 × 10 −9 m 2 /s.We consider the case of zero applied electrostatic field.
The force Eq. ( 11) acting on the QP reduces to The toy model ( 12) simplifies the comparison of the analytical calculations Eq. ( 15) with the PMF W N computed [25] in simulations from the equilibrium distributions P N via up to an arbitrary offsetting constant W of f .We consider explicitly the two distinct PMFs when there are either 2 or 3 individual K + ions in the channel.

Brownian dynamics simulation details
We consider two ionic species of opposite charge (K + and Cl -) in contrast to Ref. [17].The simulation domain (Fig. 1, left) includes the channel (narrow cylinder) of length 4nm and two reservoirs.The radius of each reservoir is 2nm, and its height is 2nm.
In the filter, a harmonic radial potential (k = 10k B T / Å2 ) is applied to ensure that ions move one-dimensionally.A harmonic repulsive field is applied when an ion approaches a domain boundary or the membrane closer than 0.3 Å.
The equations of motion of individual ions Eq. ( 1) are solved numerically using the method of Ref. [20].This offers third-order accuracy and does not suppose any restriction on the time step ∆t, which is ∆t = 0.2ps.The forces acting on individual ions are calculated from Eq. ( 2).Each simulation covered 10 µs of the ion's dynamics.
At small distances, large forces may develop due to the presence of the Coulomb and strong repulsion terms in Eq. ( 14).Straightforward application of the integration scheme will lead to so-called long jump exceptions [26].To accommodate these, an adaptive time step was used, namely, when the distance r between ions became smaller than a threshold value r 0 = 2.8 Å, the Euler scheme with the time step at the next iteration was ∆t * = ∆t (r/r 0 ) κ .As the forces grow as 1/r 10 (see Eq. ( 14)), the condition κ > 10 should be applied.We chose κ = 14 as a trade-off between accuracy and computational speed, although other choices are admissible as well.In the vast majority of cases, one adjustment suffices to unravel the ions, and the next iteration runs with the standard time step ∆t.The fraction of corrected steps was < 0.1%.
The number of ions in the simulation was maintained constant, with stochastic boundary conditions applied [27].This meant that, when an ion crosses the channel, say from the right, leaving it on the left, a leftmost ion from the left reservoir was transplanted to a randomly chosen rightmost position of the right reservoir.

Simulation results
In order to calculate the effective potential from Eq. 11, one has to transform the coordinates {z m } → {q N , p} and thus derive W N (q N , p).This method becomes attractive in the light of modern front-end studies of nanodevices which rely heavily  11) (dots and crosses), and theoretically via Eq.( 15) (solid lines), for the quasiparticles consisting of 1, 2, and 3 individual ions.U ef f 1 (solid grey line) coincides with U 0 in Eq. ( 13).This diagram demonstrates how the multi-ion PMFs W N , shown in Fig. 3, produce the 1D effective potential of the quasiparticle using formula (11).
on MD.Thus, PMF maps for many nanoscale systems have already been built.
On the top row of Fig. 3 we illustrate the PMFs of Eq.( 12), obtained from BD simulations for two (a,b) and three (c,d) ions.Between jumps, ions move along the permeation pathway (shown by a solid black line).For instance, considering a two-ion case, an entering ion pushes its neighbour, so that the latter eventually leaves the channel.Theoretical and simulation predictions match well (not shown), as they do in Figs.1-2 of Ref. [16].
Under the coordinate transformation Eq. ( 3), the PMFs from the top row in Fig. 3 are transformed into the PMFs of a quasiparticle, bottom row.The two-ion permeation process outlined above can now be reinterpreted in terms of permeation by a quasiparticle.Following the path on Fig. 3(c), one finds that the distance between ions decreases, passes through a minimum and then increases again, being accompanied by the overall displacement of the QP to the right.The spatial displacement of the QP eventually leads to charge translocation through the filter, i.e. to the electric current.
Finally, we compute the effective potentials U ef f N .As shown in Fig. 4, the implementation of formula (11) shows very good agreement between BD simulations and the theoretical predictions Eq. ( 15).One can see that the 2-ion effective potential (solid blue line) inherits the minimum from the 1-ion potential well U 0 (solid grey line), but becomes shallower and wider.This occurs due to the repulsive interac- The quasiparticle's effective potential from the PMF 11 tions which increase the ion-ion distances.Furthermore, the 3-ion effective potential (solid black line) become even wider and gains two local minima located around ±0.7nm.The increased width of the effective potential and its flatness imply that the QP is less localized.This energy landscape thus defines the spatial range over which the charge is transferred during a permeation event.

Discussion
We have derived the effective potential of a quasiparticle [17] from the potential of the mean force [6].The notion of quasiparticles allows one to reduce the collective motion of strongly interacting ions to the motion of a single particle in this effective potential.The inclusion of the PMF couples the effective potential with the energy landscape seen by the ions.Thus, the dynamics of the quasiparticle can now be connected to the structure of a real channel.
Application of the derived relation, demonstrated to a toy analytical model of interactions in the KcsA channel, enabled comparison of separate 2-and 3-ion analytical PMFs with those obtained from Brownian dynamics simulations, and thus simplified the comparison between theory and simulations.The method has direct application to modern molecular dynamics studies of nanoscale systems, in which PMFs are the primary targets for study of permeation mechanisms.
It was shown that during permeation a quasiparticle undergoes both diffusion and discrete jumps in space, motion leading to the electric current -the major experimental observable.The notion of the quasiparticle therefore couples the molecular dynamical calculations with the experimentally measured current-voltage and current-concentration characteristics of the channel.
Further work being planned includes BD simulations incorporating MDgenerated PMFs in e.g. the real KcsA channel [5,9].It will be important to see how the K + /Na + selectivity is reflected in the properties of the quasiparticle.Another direction of development is the connection with continuous methods.This suggests writing a set of coupled differential Chapman-Kolmogorov (DCK) equations for each separate occupancy state, describing diffusion and transitions of particles between states.Integration of the set of coupled DCK equations along the channel axis would provide a further link to the kinetic rate theory [28].

Fig. 1 .
Fig. 1.Left panel: Scheme of the simulation domain.The selectivity filter (narrow cylinder), reservoirs (wide cylinders) and ions of both signs (orange and blue spheres) are shown.The structure of the KcsA selectivity filter is shown in the licorice representation.Right panel: slices of the 3D PMF in the KcsA channel with the permeation pathway shown in black.Numbers 1-5 correspond to permeation landmarks from Ref. [5].Original data, also used in Ref. [5], has kindly been provided by Dr. D. Medovoy and Prof. B. Roux.

Fig. 2 .
Fig. 2. (a)Example of QP motion in Brownian dynamics simulations, illustrating a 3-ion QP that evolves into a 2-ion QP.Individual ionic trajectories (black and grey traces) result in the motion of their corresponding QP (red).At time t ∼ 0.5ns the rightmost ion leaves the channel; two ions remain, and the QP consequently jumps from q 3 to a new position q 2 .(b) The connection between the QP coordinates before (plotted on the abscissa) and after (on the ordinate) a jump.Each BD data point describes a single jump q 3 to q 2 (a circle) or q 2 to q 3 (a cross).The set of q 3 to q 2 points fall on the steeper pair of parallel lines, each point corresponding to an individual ion's entry/exit through either the left (dashed lines) or the right (solid lines) edge of the channel, as illustrated in the inset where individual ions are grey and the larger orange circle represents the QP.The set of q 2 to q 3 points fall on the similar but shallower pair of lines as described.In each case, the lines represent the theoretical predictions of Eq. (8).The red square indicates the q 3 to q 2 jump shown in panel (a), when the rightmost ion exits the channel (pale grey trajectory).(c) and (d) show the equilibrium distributions for 2-ion (orange) and 3-ion (blue) QPs.
Fig. 2. (a)Example of QP motion in Brownian dynamics simulations, illustrating a 3-ion QP that evolves into a 2-ion QP.Individual ionic trajectories (black and grey traces) result in the motion of their corresponding QP (red).At time t ∼ 0.5ns the rightmost ion leaves the channel; two ions remain, and the QP consequently jumps from q 3 to a new position q 2 .(b) The connection between the QP coordinates before (plotted on the abscissa) and after (on the ordinate) a jump.Each BD data point describes a single jump q 3 to q 2 (a circle) or q 2 to q 3 (a cross).The set of q 3 to q 2 points fall on the steeper pair of parallel lines, each point corresponding to an individual ion's entry/exit through either the left (dashed lines) or the right (solid lines) edge of the channel, as illustrated in the inset where individual ions are grey and the larger orange circle represents the QP.The set of q 2 to q 3 points fall on the similar but shallower pair of lines as described.In each case, the lines represent the theoretical predictions of Eq. (8).The red square indicates the q 3 to q 2 jump shown in panel (a), when the rightmost ion exits the channel (pale grey trajectory).(c) and (d) show the equilibrium distributions for 2-ion (orange) and 3-ion (blue) QPs.

Fig. 2 (
Fig. 2(b) represents Eq. (8) for the QP's coordinates immediately before and after a jump.Depending on the side of the entry/exit, the positions of the QP before and after form two parallel lines.The nonuniform distribution of original and final positions is due to the presence of interactions, and their centroids define which transition occurs.For instance, a 3-ion quasiparticle resides in two more probable positions (maxima of distributions in (c-d)).This means that the three ions are located closer to the inner side of the channel, and one would expect the innermost ion to leave the channel.This does happen (transition 3 → 2, red square), with the QP relocating to the centre.We also note that in the KcsA selectivity filter K + ions reside at the binding sites formed by carbonyl oxygen atoms[19].Strong interactions with the latter define the localization of individual ions very precisely[4,5,9].The exact localization of individual ions results in a much sharper, almost discrete, distribution of quasiparticles during permeation.

9 Fig. 3 .
Fig. 3. Two-dimensional (left column) and three-dimensional (right column) PMFs W N , derived from BD simulations, represented in terms of the coordinates of individual ions (top row) and quasiparticles (bottom row), the significance of the colours being indicated by the colorbar.Slices of 3D PMFs are shown.Orange sleeves trace the surface of constant potential energy 1k B T .The permeation paths are indicated by black lines.

Fig. 4 .
Fig. 4. The effective potential U ef f N calculated from the BD simulations using Eq.(11) (dots and crosses), and theoretically via Eq.(15) (solid lines), for the quasiparticles consisting of 1, 2, and 3 individual ions.U ef f