lipid bilayers provide osmotic barriers between cells and their surroundings and set the stage for operation of various transmembrane proteins such as ion channels and ion pumps that maintain optimal intra- and extracellular ionic composition . therefore , ion permeation across and partitioning into lipid bilayers is important for cellular function . understanding the mechanism of ion permeation and interactions of ions with lipid bilayers will ultimately help design new ways of delivering and maintaining at optimal levels small drug - like molecules inside cells . however , lipid bilayers are highly anisotropic and heterogeneous fluids , which makes detailed experimental studies of ion permeation and interactions with lipids difficult . on the other hand more convoluted , macroscopic properties such as ion permeation or leakage rate across membranes have been studied extensively . for example , it has been shown that alkali metal ions permeate lipid bilayers differently from halide ions , with the overall permeation rate for the cl being faster than that for na . furthermore , the permeation rate for cl , unlike na , has been found to exhibit a characteristic ph - dependence.(5 ) a number of attempts have been made to deconvolute the details of ion permeation mechanism from these macroscopic observations and to explain the differences between various ions . different mechanisms have been proposed to explain the physical basis of ion permeation across lipid bilayers . ions ( except for h and oh for which a unique mechanism has been proposed ) are believed to permeate by one or a combination of the following mechanisms : i ) solubility - diffusion , ii ) pore - mediated , iii ) flip - flop , or iv ) assisted transport . of these , the solubility - diffusion and pore - mediated mechanisms are the most popular . these two mechanisms exhibit quite different behaviors with respect to various defining parameters . the solubility - diffusion model is considered relatively crude because it does not take into account the anisotropy of lipid bilayers . therefore , more elaborate theories , such as the barrier domain model that accounts for heterogeneity and designates a specific rate limiting region in the bilayer , have been proposed.(13 ) in light of these mechanisms , halide permeation , which depends on bilayer thickness and ion size , is consistent with permeation of hydrated ions by the solubility - diffusion rather than pore - mediated mechanism.(14 ) on the other hand , alkaline ion permeation is best described by a combination of the solubility - diffusion and pore - mediated mechanisms . note that passage of hydrated ions through pores is often assumed to occur without a direct ionlipid interaction . the pore - mediated mechanism has an implicit connection with the flip - flop mechanism.(4 ) ions could take advantage of the flip - flop mechanism in cases where the energy costs for moving a single lipid and a single ion across a bilayer are comparable . however , this would require binding of ions with lipid molecules.(16 ) certain observations can not be explained by either solubility - diffusion or pore - mediated mechanisms . hence , to rationalize higher rate and ph - dependence of cl permeation , an assisted transport model involving carriercl complex formation has been postulated . following this principle , artificial carriers have been designed that successfully shuttle cl ions across membranes.(17 ) furthermore , contrary to earlier works based on continuum electrostatics theories a growing number of recent papers argue that ion permeation could in fact be assisted by ion - pairing . because of the labile nature of the interactions ,  transient interfacial ion - pairing  has been proposed to explain shuttling of hydrophilic ions across bilayer . experimental measurements of ion - pairing constants are difficult because they depend on position of the ions within the highly anisotropic bilayers.(22 ) recently , some progress has been made to probe interfacial(23 ) and cross - bilayer ion distributions.(24 ) concurrently , a number of microscopic studies have provided details of individual ion permeations across lipid bilayers on the atomic level . these latter studies have employed both free and biased molecular dynamics ( md ) simulations . biased md simulations are particularly useful because they yield potentials of mean force ( pmfs ) that allow a more detailed , quantitative characterization of permeation events . because ion - pairing along the lipid bilayer is difficult to study experimentally , we turn to pmf calculations instead . although pmf calculations have been used to study permeation of individual ions across lipid bilayers,(39 ) to the best of our knowledge , they have not yet been used to study corresponding ion - pairing . in this paper we attempt to fill this void by using accurate pmf calculations to study ion - pairing and the assisted permeation mechanism across the dmpc bilayer . can ion - pairing facilitate permeation ? what is the mechanism of the ion - pair permeation across the bilayer and is it different from that of the individual ions ? to give the reader the sense of current state of the microscopic picture of ion interactions with lipids from molecular dynamics simulations , we briefly review recent studies that are essential for understanding of the ion permeation across lipid bilayers . free md simulations have provided valuable insights into the interactions between ions and lipids at the lipidwater interface ( see ref ( 2 ) for a review ) . many groups have studied ionmembrane interactions initially using united atom models of lipids and , recently , more expensive , fully atomistic models . most of these simulations have been performed with constant surface tension imposed on the lipid bilayer and reported shrinkage of the area per lipid upon increasing salt concentration , resulting in a more ordered and thicker bilayer . interestingly , na has been found to bind to lipid headgroups , while cl has predominantly stayed in bulk water . this is particularly intriguing because the membranes studied were made of zwitterionic lipids with zero net charge , such as dmpc , dppc , and popc . the behavior of na and cl at the lipidwater interface has also been investigated under constant area per lipid conditions where no shrinkage of the bilayer would have been possible . for example , 20 ns long md simulations of a popc lipid bilayer in nacl solution by sachs et al . have suggested that both na and cl penetrate the bilayer to within 14.5 and 17.5  from the bilayer center . larger anions have been found to penetrate even deeper ( to within 12.5  from the center ) . have performed 200 ns long md simulations of a popc bilayer with the united atom gromos and the all - atom charmm force fields . these simulations confirm that na binds to the carbonyl region of the membrane , whereas cl stays in water . overall , numerous computational studies have demonstrated that na gets absorbed at the headgroup region while cl remains in solution . constant surface tension studies have also reported that the area per lipid shrinks significantly upon na binding to the bilayer . all these studies have documented long relaxation timescales and hence convergence difficulties for ionwatermembrane interactions . they have indicated that obtaining converged properties for ions in lipid bilayers requires 20 to 100 ns of simulation time . because ion permeation across a membrane is a rare event , with the free energy barrier much greater than thermally available kbt , it can not be observed with free md simulations . therefore , quantitative studies of ion permeation require use of biased simulations that allow computing accurate pmfs . we prefer real space approaches to computing pmfs such as umbrella sampling(49 ) with the weighted histogram analysis method ( wham ) , the umbrella integration method(52 ) with harmonic fourier beads ( hfb ) , or constrained simulations with adaptive biasing force ( abf ) . at present , accurate pmf calculations in lipids are relatively scarce in the literature , reflecting the convergence challenge . in particular , early attempts to calculate the pmfs for ion permeation across a lipid bilayer have identified areas of slow convergence around the middle of the bilayer.(34 ) they found that overall convergence is difficult to achieve and that the pmfs remain asymmetric even after 50 ns simulation time per window . in a pioneering work , pohorille et al.(12 ) computed individual pmfs for na and cl ion permeation through a glycerol 1-monooleate ( gmo ) bilayer modeled with united atom forcefield . the pmfs were computed using a predecessor of the abf method and were fairly short ( between 0.5 and 2 ns per window ) . these authors also suggested that pmfs computed with united atom models needed empirical adjustments . recently , voth et al . computed pmfs for na and oh ions permeation through a dmpc bilayer also employing a united atom forcefield.(9 ) they used umbrella sampling with wham and rather short md runs ( up to 1 ns ) to assemble the pmfs . , their pmfs could be compared to experimental studies without empirical adjustments . nonetheless , the permeation barrier for na was quite high , on the order of 25 kcal / mol . these authors also addressed an interesting permeation mechanism where ions on different sides of the bilayer could facilitate each other s permeation , which indirectly probed the role of ion - pairing in the permeation mechanism . overall , pmfs for ion permeation of lipid bilayers computed with biased md simulations have significantly advanced our understanding of the permeation mechanism at the atomic level . moreover , biased md simulations have identified trails of water molecules connecting the solvated ion inside the bilayer to the lipid / water interface , so - called water fingers , pockets , or funnels . somewhat different from the lipid bilayer permeation pmfs , but valuable for understanding ion permeation , was the work by hummer and co - workers.(28 ) these authors employed slabs of nanotubes as models of bilayers . the resulting nanotube layers lacked the heterogeneity and the fluidity of lipids but instead provided well - defined pores . thus , changing the effective diameter of the pore from 5 to 10  lowered the na ion permeation barrier from about 29 kcal / mol to about 1 kcal / mol . several other pmf calculations related to lipid bilayers are helpful for building a complete picture of ion permeation . these include permeation of lipid molecules , water , and h. lipid molecule pmf computation involves the removal of a lipid molecule from a preassembled lipid bilayer . such data are relevant for assessing the flip - flop mechanism and have been computed by several researchers . these studies reveal barriers as high as 19 kcal / mol for a regular dppc molecule(34 ) and about 30 kcal / mol for an h - ras anchor that has a third fatty acid chain(32 ) moving across the bilayer . thus , the free energy for moving lipid molecules with multiple fatty acid chains is about 10 kcal / mol per chain . water does not obey the solubility - diffusion mechanism and instead diffuses across the membrane via transient pores or fluctuating defects . experimentally measured activation free energies for water permeation are between 4 and 9 kcal / mol . several authors have computed pmfs for water permeation and found barriers of about 47 kcal / mol for dppc and dlpe membranes(27 ) in good agreement with experimental values . h permeation pmf computation reflects the anomalously high ( compared with all the other ions ) permeation rate of h across lipid bilayers . computing pmfs for proton permeation is somewhat more involved than for most classical ions because the proton requires quantum chemical description to allow charge delocalization . voth and co - workers computed a pmf for h permeation using a multistate variant of empirical valence bond ( evb ) model and found that h can permeate membranes much more readily than na but less efficiently than water itself . using biased simulations seems to be the ultimate way to study such rare events as ion permeation across membranes . nevertheless , there have been attempts in the literature to achieve ion permeation by alternative means , such as by creating transmembrane pores.(25 ) clearly , creating a water - filled pore that pierces the bilayer requires work that strongly depends on the pore size . interestingly , however , the presence of ions was found to greatly reduce stability of the pores . moreover , similar to the case with the lipidwater interface , na and cl permeations across the water - filled pores show some differences . thus , for small pore sizes , na and cl have similar permeation rates , but for pore sizes larger than 15  , cl permeates much faster . computed pmfs for na and cl ion permeation through metastable pores show a relatively strong binding for na and repulsion to only weak binding for cl . in the absence of the pore , the na permeation barrier has been estimated to be about 36 kcal / mol . the authors of that study concluded that na will permeate through transient pores , while cl will likely permeate by the solubility - diffusion mechanism . in certain cases , thus , some simulations have been engineered to cause ion permeation within the typical md simulation times by preparing molecular systems far from equilibrium . in this case , an unrealistic charge imbalance created using two lipid bilayers setup results in pore formation across the membrane and allows studying subsequent ionic leakage through such a pore.(44 ) although unrealistic , this process can still provide valuable information on the behavior of lipid bilayers subjected to fluctuations in ionic / charge concentrations . for example , similar to the work by marrink and co - workers , recent studies of ion permeation through pores created by nonequilibrium processes by gurtovenko et al . observed that na and cl leak at around the same rate under conditions conducive to narrow pores.(46 ) these authors also compared all - atom and united atom forcefields and found similar results . thus , md and in particular biased md simulations have dramatically advanced our understanding of the permeation mechanism for individual ions across bilayers . these studies have confirmed and considerably augmented the predictions of the continuum electrostatics modeling concerning individual ions . however , very little is known about the effect of ion - pairing and the microscopic behavior of ion - pairs during permeation across lipid bilayers . this , in part , is due to the fact that continuum electrostatics modeling ruled ion - pairing out at least in the case of bare ions . given the recent reports invoking assisted transport of ions across membranes , and the obvious disagreement with the early predictions derived using continuum models , here we explicitly test the possibility of assisted transport using biased md simulations with all - atom models of water , lipid bilayer , and ions . in particular , we assess the possibility of the  transient interfacial ion - pairing  that is essential to assisted ion shuttling across lipid bilayers . therefore , we compute high quality pmfs for the na and cl ion permeation across an all - atom model of the dimyristoyl phosphatidylcholine ( dmpc ) lipid bilayer individually , and as a nacl ion - pair free to dissociate . we would like to note that although this particular pair represents two hydrophilic ions , the mechanism of the ion - pair permeation might be more general . our simulation system comprised 54 dmpc lipids , 1747 ctip3p water molecules and a single pair of na and cl ions . thus , the total number of atoms n in the system was 11  615 . molecular dynamics simulations were performed in the npnat ensemble with t = 310 k. the area a was fixed using the area per lipid ratio of 59.7  / lipid and the normal pressure pn = 1 atm . we chose the constant area npnat rather than constant surface tension npnt ensemble for computational convenience and because they give similar results for unperturbed , fluid lipids . note that the area a and the surface tension  are thermodynamically conjugate variables and both impose constraints on the system.(82 ) ideally , we would like not to use any constraints and just perform constant pressure npt simulations ; however , lipid forcefields continue to undergo modifications to achieve this goal.(83 ) in the mean time , we feel that our choice of the ensemble for permeation studies is reasonable given that at the chosen simulation conditions the dmpc lipid is in the fluid state and that the permeation pmfs we pursue reflect the forces acting on small molecules along the z - axis whereas the lateral x - y forces average out to zero . the ultimate test of the effect of these approximations will have to be assessed once the npt - capable lipid forcefields become available . we used namd 2.6(84 ) for running md simulations and charmm c32b1 for reconstructing the pmfs using the hfb umbrella integration procedure and also for analysis . the temperature was maintained with the berendsen thermostat using a coupling parameter of 2.5 ps , whereas the pressure normal to the bilayer was maintained at 1 atm by the langevin piston method with the piston mass of 100 amu and langevin collision frequency of 50 ps . nonbonded interactions were computed using particle mesh ewald method with 8.5  real space cutoff for electrostatic interactions and the switching functions between 8.5 and 10  for the vdw . however , the multiple time - step method was employed for the electrostatic forces with full electrostatic interactions computed every other step . the nonbonded interaction list was constructed using a cutoff of 12  and updated every 20 steps . the covalent bonds involving hydrogen atoms were constrained using the shake algorithm with a tolerance in the bond length of 10  . the md simulations for the pmf calculations were initiated by dragging either ion or the center of mass of the ion - pair from the bulk on the positive half of the z - axis toward the bilayer center in increments of 0.5  . the system was equilibrated at every new position of ions for 200 ps . at all times the z - component of the center of mass of all the heavy atoms of the bilayer was restrained at the box center . during the pmf calculations for an individual ion , the other ion was kept restrained at a remote position in the corner of the simulation box with coordinates 20.0 , 16.0 , 33.5  . for each individual ion pmf , we used 78 beads covering the z - coordinate range of 5.0 to 33.5  , and for the ion - pair we used 67 beads spanning the z - coordinate range of 5.0 to 28.0  . the ion - pair range has been truncated to avoid problems arising from wrapping one of the ions to the other side of the simulation box to satisfy periodic boundary conditions that drastically changes the position of the center of mass of the ion - pair upon restart of the md simulations . in a few cases where such problems appeared , the corresponding simulations were restarted from adjacent windows following re - equilibration at the target center of mass position for 200 ps . we stress that this problem only arises due to the need to restart simulations and would not have been present should we use a single continuous trajectory to collect the mean forces or an alternative coordinate wrapping scheme . during the simulations for the ion - pair pmf , the ion - pair was allowed to dissociate and rotate to integrate the corresponding contributions out of the pmf along the center of mass coordinate . however , for the ion - pair we employed an additional harmonic boundary restraint of cylindrical symmetry that prevented ions from crossing the boundaries of the primary simulation cell in the x - y direction . 

 where di = ( ( xi  xcmref ) + ( yi  ycmref ) ) . here m and rcmref are the mass of the ion - pair and the reference position of its center of mass ; mi and ri are the mass and the position of an ion ; and d and k are the maximum allowed separation between the ions after projection onto the x - y plane and the force constant . this restraint allowed the ions to separate freely up to 16  measured after projection onto the x - y plane . although , this restraint is not perfectly orthogonal to the z component of the center of mass of the ion - pair , we do not expect it will have a significant effect on the pmf because of relatively infrequent sampling of the boundary by the ions . the force constants for the umbrella potentials used in the present pmf calculations were as follows : 0.5 kcalg weighted by the mass 32690.09828 g / mol of all the heavy atoms of the dmpc bilayer ; 10.0 kcalg weighted by the mass 22.98980 g / mol , 35.45000 g / mol and 58.43980 g / mol for the na , cl and the nacl ion - pair , respectively . to accumulate the mean forces , for each bead we performed one hundred consecutive 200 ps md runs , adding up to a total of 20 ns simulation time per bead . thus the total simulation time used to compute pmfs was 1560 ns for each of the two ions and 1340 ns for the ion - pair . the individual and ion - pair pmfs have been integrated using the ggahfb method(54 ) with 78 and 67 fourier basis functions , respectively . in agreement with previous theoretical work , we find that both ions and their ion - pair enter the bilayer solvated with a trail of water molecules connected to the lipid / water interface on one side of the bilayer ( see top panes in figure 1 ) . the trailing water molecules have been called fingers in the early work on single ion permeation across simple water / lipid interfaces . these trailing water molecules significantly perturb the lipid bilayer by dragging its polar head groups further toward the center ( see the middle and the bottom panes in figure 1 ) . representative electron densities for the dmpc bilayer and surrounding water along the na , cl , and nacl ion - pair permeation pmfs . the top row shows complete densities of the dmpc ( solid lines ) and tip3 ( dashed lines ) water molecules ; the middle row displays the densities for the carbonyl groups ( c = o , solid lines ) and positively charged choline ( chol , dashed lines ) ; the bottom row displays the density for the negatively charged phosphate groups ( po4 , solid lines ) . the polar head groups can also contribute to the first solvation shell of the ions inside the membrane . in fact , this coordination persists even at the center of the bilayer as evident from figure 1 . in contrast , cl ion does not directly coordinate either carbonyl or phosphate groups . its interactions with the choline group appear insufficient to cause choline to follow the cl into the center of the membrane . in the case of the ion - pair , in particular , na ion still coordinates carbonyls and phosphates , whereas cl ion does not . as a result , the electron density profiles depicted in figure 1 for the ion - pair look like a weighted average of the individual ion electron density profiles . when the ions reach the center of the membrane , the water fingers have to switch from one side of the lipid bilayer to the other by means of thermal fluctuations due to symmetry requirements.(96 ) this stochastic process appears to have a time scale that is on the order of nanoseconds . sometimes , when in the center of the bilayer , the first solvation shell of the ions can disconnect from the bulk , completely sealing off the ion with its nearest solvation shell in a capsule . these fluctuation events are quite rare and are relevant to forming the water defect on the other side of the bilayer . figure 2 depicts three pmfs for each of the ions and the ion - pair computed for different total simulation times . specifically , the thick , solid lines show the final 20 ns cumulative pmfs , and the thin , dashed lines refer to the first and last 10 ns pmfs . the pmfs give the following free energy barriers for the ion permeation : 21.9 kcal / mol for na , 23.6 kcal / mol for cl , and 27.6 kcal / mol for nacl . these barriers are located at 0.27  , 0.89  , and 0.26  along the z - axis , respectively . , we compare the free energies at the symmetric z - positions in the pmfs . specifically , we look at the quantity gz = g(z )  g(z ) at z = 5.0  , after 20 ns of md simulations . thus , the computed g5 difference for the nacl ion - pair is 2.4 kcal / mol ; for the na ion the difference is 0.2 kcal / mol , and cl ion shows the largest difference of 4.0 kcal / mol . comparing the pmfs for the first and the last 10 ns ( thin dashed lines ) , we observe that the g5 differences gradually decrease , further suggesting that the computed pmfs are not fully converged even after 20 ns . thus , the error bars on the free energy barriers for ion permeation of lipid bilayers appear to be fairly high and could be at least half of the g5 . pmfs for na , cl , and nacl ion - pair permeation across the dmpc bilayer . thick lines represent the cumulative pmfs computed by umbrella integration with harmonic fourier beads method using 20 ns simulations per window . thin short - dashed ( a ) and long - dashed ( b ) lines correspond to the pmfs computed using the first and the last 10 ns of the 20 ns simulation per window , respectively . although fully converged , perfectly symmetric pmfs are very difficult to achieve , the region in the pmfs that causes the largest convergence problems can be easily identified by shifting the bulk water values ( at the largest z - value ) of the pmfs to zero . a simple inspection of figure 2 identifies the region with the largest deviations in the pmfs as the middle of the bilayer . this is where switching of the water fingers from one leaflet to the other must occur by way of thermal fluctuations . these fluctuations are rather rare and hence cause severe convergence problems . nevertheless , the area of the pmfs between the maximum and the bulk water z - values ( roughly 0 < z < 33.5  ) , where the water fingers remain consistently connected to the lipid / water interface , appears converged to within 0.25 kcal / mol . it is worth noting that such convergence can be achieved already at 1020 ns of simulation time per window and is quite satisfactory for the purposes of this paper . we can use the na pmf to quantify the binding free energy for na with dmpc . in agreement with previous free md simulations , we find that na has a broad basin in the headgroup region ( between about 14 and 26  ) , resulting in a weak ( 0.4 kcal / mol ) binding at 19.1  . the fact that these values are comparable with the sampling error ( 0.25 kcal / mol ) of the pmf makes precise identification of the positions of the minima difficult . the small values of these minima ( on the order of kbt ) and the absence of a desolvation barrier explain the spontaneous na / lipid interactions observed by unbiased md simulations . unlike na , the cl ion interactions with the dmpc bilayer are purely repulsive , with cl preferentially staying in solution . although the barrier for the cl permeation is higher than that for na , the slope of the cl pmf is consistently smaller than that of na in the region between the center of the bilayer ( z = 0  ) and about z = 15  away from it ( see figure 2 ) . in the simulations of the ion - pair permeation across the bilayer , the na and cl ions are free to dissociate according to the underlying free energy . to assess the effect of the dmpc lipid bilayer on the free energy of ion - pairing , we compute interionic distance histograms for each of the 67 beads in the nacl ion - pair permeation pmf . although the limited sampling prevents us from converting this data into the corresponding nacl distance pmfs , the observed changes in the uncorrected distance distribution along the bilayer normal provide useful information . thus , figure 3 shows the top view of a 3d plot of the nacl distance distribution histograms with respect to position of the center of mass of the nacl ion - pair in the bilayer . normalized uncorrected histograms for the nacl distance distributions as a function of the bilayer permeation depth . color - coding was done on a logarithmic scale to accentuate the ssip relative to cip by exponentiation of the normalized densities to the 1/5 power . as seen in figure 3 , the ion - pair starts out in the bulk solution with a typical distance distribution : two well - defined discrete states at short distances representing contact- ( cip ) and solvent - separated - ion - pair ( ssip ) , and a broad continuum of states beyond that . as the ion - pair enters the bilayer , its distance distribution gradually loses the continuum states and shrinks down to the two discrete ion - pair states . near the center , the distributions are dominated by the cip with only a minor contribution from the ssip . thus , even in the most hydrophobic region of the bilayer , the ion - pair exhibits two - state behavior with the cip maximum around 2.52.7  and the ssip maximum around 4.44.6  . however , the latter becomes noticeably compressed as the ion - pair gets closer to the center of the bilayer . for comparison , at 298 k in water the na and cl are 2.6  and 5.1  apart in cip and ssip , respectively . even though the ions occasionally dissociate and separate by as far as 10  near the center of the lipid bilayer , most of the time they return to the bound state . another proof of the increasing binding free energy is the reduction of the distance in the ssip that agrees well with the hammond postulate.(98 ) the increase in bindng energy will be quantified in the following paragraphs . note also that in figure 3 we accentuate the ssip distribution by using logarithmic scaling in the color map . the stability of the ion - pair in the middle of the lipid bilayer is undermined by a stochastic formation of water - filled pores that pierce the lipid bilayer . formation of such pores is conducive to dissociation of the ion - pair along the bilayer normal , with ions separating to the opposite sides of the bilayer . although rare , we did observe two such ion separation events during our simulations ( at windows with z = 1.5 and z = 2.5  ) . these events also affect the pmf calculation , as once the ions escape into the bulk , the mean forces on the center of mass of the ion - pair approach zero . after these rare ion separation events have been detected , trajectories have been restarted from the nearby windows where ions remained associated . to assess the orientation and direction for dissociation of the ion - pair as a function of the bilayer penetration depth , we have computed the angle distribution between the clna vector and the z - axis shown in figure 4 . clearly , orientation of the ion - pair in the lipid bilayer is strongly anisotropic and changes with respect to the bilayer normal as the penetration depth increases . normalized histograms for the angle of the cl - na vector about the z - axis ( naclz ) in degrees as a function of the bilayer permeation depth . interestingly , near the center , the ion - pair precesses about the bilayer normal , with cl positioned closer to the center than na . crossing the center of the bilayer into the other leaflet , the ion - pair orientation undergoes a sharp transition to maintain cl proximity to the bilayer center as is seen from the abrupt change in the angle distribution . similarly , as the ion - pair moves further than 15  away from the center of the bilayer , the ions swap places and the horizontal component of the dissociation becomes significantly more important as confirmed by the dense population of the area near 90. however , because at such long distances cl is expelled from the bilayer into the bulk while na remains bound , the angle distribution is skewed toward higher angles . thus another more gradual transition occurs around 15  that causes the ion - pair dipole to flip . simple continuum electrostatic arguments rule out conventional ion - pairing as a possible mechanism of ion permeation , at least in the case of ions permeating without their water solvation shells . however , compelling evidence suggests that certain charged carriers facilitate ion transport via the ion - pairing mechanism . one of the requirements of such facilitated permeation is  transient interfacial ion - pairing  . using the three pmfs computed in this study , we are able to assess for the first time the change in the ion - pairing free energy as a function of the ion - pair permeation depth . to compute the change in the nacl binding free energy as a function of the position of the ion - pair inside the bilayer , we use the detailed pmfs for individual ions and the ion - pair as follows first , we find the average positions of the individual ions for each of the beads in the ion - pair pmf from the corresponding trajectories . then for the determined average positions we find the corresponding free energy values using the individual ion pmfs . the final estimate of the change in the binding free energy between na and cl along the membrane normal with respect to bulk water is shown in figure 5 . except for a rather broad , short hump at the headgroup region , the magnitude of the binding free energy increases gradually upon going toward the center of the bilayer . the observed hump is due to preferential na binding at the headgroup region and indicates that the association of the two ions in that region is not as strong as in water . going over the hump , the electrostatic interactions between the ions become stronger and the magnitude of the binding energy inside the nonpolar region of the bilayer increases rapidly . as inferred from figure 5 , at the center of the bilayer , the absolute change in the ion binding free energy reaches a maximum of 17.1 kcal / mol near the bilayer center . if , for simplicity , we estimate the change in the binding free energy from the maxima of the three pmfs , we get a 17.9 kcal / mol absolute value ( 27.6  21.9  23.6 kcal / mol ) . evidently , ion - pairing in the bilayer reduces the barrier for the ion - pair permeation from the sum of the individual na and cl permeation barriers of 45.5 kcal / mol ( 21.9 + 23.6 kcal / mol ) to 27.6 kcal / mol ( 45.5  17.9 kcal / mol ) . however , in the present case the ion - pair permeation barrier is still larger than the barrier for the individual ions . change in the binding free energy in kcal / mol between na and cl ions relative to bulk water as a function of the bilayer permeation depth . averaged positions of individual ions are shown on the right axis , whereas the change in binding free energy is shown on the left axis . the vertical and horizontal black lines pass through values of zero of the three axes . to clarify our calculations above and to further explore if pairing of ions of opposite charge can facilitate ion permeation under different circumstances , consider a thermodynamic cycle for an ion - pair transfer from polar into apolar solvent shown in figure 6 . we label one of the ions a and the other b. we can write down the following equation for the free energy of ion binding in each phase . illustration of the thermodynamic cycle for a coupled ion transfer from polar to apolar phase . for simplicity the overlap of the solvation shells in the cip state is shown in green . here the superscripts a and p refer to apolar and polar phase , respectively . defining the transfer free energies from polar to apolar solvent for the individual ions as 

 the transfer free energy for the ion - pair then is 

 where the gabp - a = gabagabp is the change in binding free energy upon ion - pair transfer from polar to apolar solvent . thus , the free energy for ion - pair transfer from polar to apolar environment equals the difference between the binding energies of the ions in the respective phases plus the energy of transfer of individual ions from polar to apolar phase . with a. for two hydrophilic ions , the binding energy is expected to be stronger in the apolar phase than in the polar phase . this leads to negative gabp - a , just like we have seen with na and cl moving from water into the bilayer . however , the individual ion transfer free energies are expected to be of positive sign and relatively large magnitude . therefore , the only case in which the hydrophilic ions can benefit from ion - pairing is when the change in the binding energy outweighs one of the individual ion s transfer free energy . according to our calculations , b. for two hydrophobic ions , the binding free energy is expected to be stronger in the polar phase . this would result in a positive gabp - a and negative individual ion transfer free energies . this case is otherwise very similar to case a. c. in the case of two ions , one of which is hydrophobic and the other is hydrophilic , the gabp - a can be close to zero . in addition , the free energies for the individual ion transfer can have opposite signs and can potentially cancel each other out to some degree . therefore , ion - pairing can be beneficial for ion permeation in this case . in a paper describing br and tba permeation across the water / nitrobenzene interface , the pmfs provided for the individual ion transfer have opposite sign and indeed cancel each other to some degree . previous md simulations have demonstrated that na partitions preferentially into the headgroup region of bilayers . in contrast , cl ion has been found to remain in solution . those studies have further indicated that quantifying the na binding to the lipid bilayer by free md simulations requires runs of tens to hundreds of nanoseconds . in agreement with these observations , the present work identifies a broad and shallow na binding well within the headgroup region of the dmpc bilayer . the largest binding free energy of na to the dmpc bilayer is 0.4 kcal / mol at 19.1  from the bilayer center ( which is very close to the maximum of the choline electron density depicted in the middle row in figure 1 ) . our calculations find no desolvation penalty for the na ion to enter the bilayer and confirm that cl ion interaction with the dmpc bilayer is strictly repulsive . in agreement with previous simulations and theoretical predictions , furthermore , they maintain the connection to the lipid / water interface on one side of the bilayer via trails of water molecules or water fingers . therefore , it is not surprising that the ion - pair also drags water into the bilayer up to its very center . na coordinates dmpc molecules directly by oxygen atoms of their carbonyl and phosphate groups . on the other hand , cl only weakly attracts the positively charged choline groups of the dmpc , which seems insufficient to drag them close to the bilayer center . specific interactions of the na with the dmpc molecules may explain lower permeation barrier compared to cl ion and preferential orientation of the ion - pair . it is worth mentioning that incorporating additional carbonyl groups into lipid molecules has been shown to significantly lower the barrier for na permeation across composite lipid bilayers by unknown mechanism.(99 ) the two observed transitions in the ion - pair orientation , particularly the one resulting in the cl positioned closer to the center of the bilayer than na , might seem surprising , given that the individual cl pmf is higher than that of na throughout . however , a closer look at the corresponding individual ion pmfs suggests that the two ions have significantly different slopes in that region ( see figure 2 ) . specifically , the cl ion has the lesser slope of the two ions despite the overall higher free energy . thus , it is more favorable to have the cl ion in front of the na ion on the approach of the ion - pair to the center of the bilayer . note also that the na ion directly coordinates the dmpc molecules by the phosphate and/or carbonyl oxygens ( as is seen from the left column in figure 1 ) in that interval and thus drags the lipid molecules along with trailing waters behind like a parachute . when the water finger finally switches to the other side of the bilayer , cl and na swap places to allow cl remain closer to the bilayer center . as mentioned above , na interacts with the phosphate and carbonyl groups more strongly than cl does with the choline groups . this is likely related to the size of the ions , as the interaction between the smaller na ion and phosphate is more effective than the corresponding interaction between the larger cl and choline . the orientation of the ion - pair or its dipole also responds to the local dipole of the membrane near the headgroup region . the second much sharper flip happens at the center of the membrane where there is no such dipole because of the symmetry of the bilayer.(82 ) therefore , it seems more likely that the orientation of the ion - pair near the center of the bilayer is due to the preferential solvation of the na rather than cl by the trailing water . in particular , the specific interactions of the na with the phosphate and the carbonyl groups that are stronger than those with water alone retain the na bound to the membrane in the headgroup region . on the other hand , the cl does not bind effectively to the membrane because of weaker interactions with the choline groups and is expelled from the membrane into bulk water . quantifying the interactions between the ions and the lipid bilayer by means of pmfs is important for understanding the ion permeation mechanism . we stress that in computing pmfs for permeation of molecules across lipid bilayers , it is essential to identify both the plateau region and the barrier region of the pmfs . as has been demonstrated in the present paper and also shown in previous papers , the turnover point or the maximum in the pmf is often shifted away from the center of the bilayer . setting the flat , bulk water region of the pmfs to zero clearly identifies the turnover point and reveals the well - converged portions of the pmfs . the offset in the z - position of the barrier likely indicates a systematic error due to limited sampling . this procedure also exposes the intrinsically problematic area of pmf convergence near the membrane center , where the trails of water or water fingers switch from one side of the bilayer to the other by means of fluctuations . we did not find experimental free energy barriers for na and cl permeation in pure dmpc bilayers.(3 ) however , permeation barriers for na and cl measured in phosphatidylserine were found to be around 23.5  2.7 and 20.8  0.4 kcal / mol at 300 k , respectively.(4 ) in contrast to this data and some computational work,(12 ) we find that the barrier for cl permeation through dmpc lipids is 1.7 kcal / mol higher than that for na . note that although the sampling error in our work is largest for cl , the pmf is well converged up to the barrier . furthermore , behavior of both na and cl ions along the ion - pair pmf , such as the relative position of ions and the orientation of their dipole during permeation , is in agreement with the pmfs of the individual ions . we therefore believe the relative order for the barriers of the na and cl permeation to be reliable . the present ion permeation and ion - pairing study can only be considered as the limit of the nonpolarizable force field . however , the effect of polarizability on the pmf of ion permeation has been demonstrated to be relatively small.(37 ) furthermore , considering experimental permeation rates for the ions , the nacl ion - pair is unlikely to undergo a charge transfer inside the membrane . parsegian also suggested that such scenario would be very unlikely because nacl remains ionized even in the gas phase.(19 ) these factors provide additional confidence in the pmfs computed in the present work . it is worth noting that measurement of the electric parameters for zwitterionic phosphatidylcholine membranes found that na diffuses across the bilayer 4.6 times faster than cl , in sharp contrast to the isotopic measurements that predicted cl to diffuse 232 times faster than na.(3 ) this discrepancy has been reconciled by invoking an additional permeation path for the cl ion , namely the  exchange diffusion  or assisted transport path.(3 ) that study proposed a charged carrier forming an uncharged complex with cl but did not exclude the possibility of  a small trace molecule acting as a carrier  . the fact that the rate of cl ( but not na ) permeation is inversely proportional to ph further supports a more complex mechanism of cl permeation than the passive , intrinsic permeation studied here . the thermodynamic cycle discussed above suggests that ion - pairing could enhance ion transport under certain conditions . one particular scenario that seems to be operating in many reported cases is pairing of ions of opposite solvent philicity . hence , as has been done by other authors,(5 ) we can tentatively suggest a solvated proton as a counterion for cl permeation for the following reasons . first , because the permeation barrier is near the center of the bilayer , we do not expect an interfering effect from a possible protonation of the lipid headgroups . in the present study we have employed an all - atom dmpc lipid model to avoid the uncertainty with the free energy profiles associated with the united atom models . in particular , we have computed pmfs for dmpc bilayer permeation by na and cl ions individually and as a nacl ion - pair ( free to dissociate ) . in agreement with previous studies , we have identified the middle of the bilayer as the most challenging region for converging the pmfs , whereas the major part of the pmf between the bulk and the bilayer center demonstrates relatively fast convergence to within 0.25 kcal / mol . the difficulties in converging the pmfs near the center arise from the stochastic nature of fluctuations that are responsible for switching the trailing water fingers from one leaflet of the bilayer to the other to satisfy symmetry requirements . the computed pmfs provide intrinsic , unassisted free energy barriers to permeation of na and cl . the fact that the free energy barrier to permeation reduces in the order nacl ion - pair , cl ion and na ion demonstrates that there is no benefit from ion - pairing in this case . furthermore , the ordering of the barriers for the individual ions agrees with electrochemical measurements but not the isotopic measurements of the ion permeation rate . thus , our calculations indirectly support the existence of an assisted ion transport across the bilayer in the case of cl ion . one of the most likely candidates for the cl ionophore would be the hydrated h ion , which would also explain the dependence of the cl permeation rate on ph . using the data from the three pmfs we have been able , for the first time , to quantify the change in the binding free energy between na and cl ions as a function of the permeation depth of their ion - pair center of mass . we have demonstrated that the na and cl ions in the ion - pair , although still capable to dissociate in the membrane , bind much tighter so that the ion - pair primarily exists as a contact ion - pair ( cip ) in the region of hydrophobic tails . the solvent - separated ion - pair ( ssip ) is populated only slightly and its nacl distance appears to be compressed near the center of the bilayer compared to that in the bulk water , in agreement with increasing binding strength . the orientation of the ion - pair inside the bilayer is strongly anisotropic and changes with the permeation depth . in particular , in the hydrophobic tail region of the bilayer the dipole of the ion - pair precesses along the membrane normal with cl ion positioned closer to the center of the bilayer than na ion . reaching the center of the bilayer causes the largest increase in the magnitude of the binding free energy between na and cl up to 17.9 kcal / mol . the results of this work further advance our understanding of ion permeation mechanisms and are of general significance for the growing body of work on assisted ion transport .
understanding the mechanism of ion permeation across lipid bilayers is key to controlling osmotic pressure and developing new ways of delivering charged , drug - like molecules inside cells . 
 recent reports suggest ion - pairing as the mechanism to lower the free energy barrier for the ion permeation in disagreement with predictions from the simple electrostatic models . in this paper 
 we quantify the effect of ion - pairing or charge quenching on the permeation of na+ and cl ions across dmpc lipid bilayer by computing the corresponding potentials of mean force ( pmfs ) using fully atomistic molecular dynamics simulations . 
 we find that the free energy barrier to permeation reduces in the order na+cl ion - pair ( 27.6 kcal / mol ) > cl ( 23.6 kcal / mol ) > na+ ( 21.9 kcal / mol ) . 
 furthermore , with the help of these pmfs we derive the change in the binding free energy between the na+ and cl with respect to that in water as a function of the bilayer permeation depth . 
 despite the fact that the bilayer boosts the na+cl ion binding free energy by as high as 17.9 kcal / mol near its center , ion - pairing between such hydrophilic ions as na+ and cl does not assist their permeation . 
 however , based on a simple thermodynamic cycle , we suggest that ion - pairing between ions of opposite charge and solvent philicity could enhance ion permeation . 
 comparison of the computed permeation barriers for na+ and cl ions with available experimental data supports this notion . 
 this work establishes general computational methodology to address ion - pairing in fluid anisotropic media and details the ion permeation mechanism on atomic level .