Validating Affinities for IonâLipid Association from Simulation against

ARTICLE pubs.acs.org/JPCA

Validating Affinities for IonLipid Association from Simulation against Experiment Benjamin Klasczyk and Volker Knecht* Max Planck Institute of Colloids and Interfaces, Science Park Golm, 14424 Potsdam, Germany ABSTRACT: Understanding biological membranes at physiological conditions requires comprehension of the interaction of lipid bilayers with sodium and potassium ions. These cations are adsorbed at palmitoyl-oleoyl-phosphatidylcholine (POPC) bilayers as indicated from previous studies. Here we compare the affinity of Na+ and K+ for POPC in molecular dynamics (MD) simulations with recent data from electrophoresis experiments and isothermal calorimetry (ITC) at neutral pH. NaCl and KCl were described using GROMOS or parameters matching solution activities on the basis of KirkwoodBuff theory (KBFF), and K+ was also described using parameters by Dang et al., all in conjunction with the Berger parameters for the lipids and the SPC water model. Apparent binding constants of GROMOS-Na+ and KBFF-K+ are the same within error and in good agreement with values from ITC. Although these force fields yield the same number of bound ions per number of lipids for Na+ and K+, they give a larger number of Na+ ions per surface area compared to K+, in agreement with the electrophoresis experiments, because Na+ causes a stronger reduction in the area per lipid than K+. The intrinsic binding constants, on the other hand, are reproduced by Dang-K+ but overestimated by GROMOS-Na+ and KBFF-K+. That no ion force field reproduces the intrinsic and the apparent binding constant simultaneously arises from the fact that in MD simulations, implicitly meant to mimic neutral pH, pure PC is usually modeled with zero surface charge. In contrast, POPC at neutral conditions in experiment carries a low but significant negative surface charge and is uncharged only at acidic pH as indicated from electrophoretic mobilities. Implications for future simulation and experimental studies are discussed.

1. INTRODUCTION Biological membranes are formed by a complex mixture of a variety of lipids and proteins and surrounded by aqueous media containing various types of ions. The concentration of these ions may be enhanced or reduced in the vicinity of the membrane if the lipids are charged, which depends only on the valency of the ions and is modeled by the self-consistent mean-field description provided by PoissonBoltzmann theory.1 In addition, ions may be adsorbed at a membrane due to specific molecular interactions that are dependent on the size of the ions, leading to ion-specific effects. Such effects can be studied theoretically using molecular dynamics (MD) computer simulations. A prerequisite for a correct description of ion adsorption is to reproduce the membrane affinity of the ions quantitatively. However, the membrane affinities of ions depend on the force field chosen.2,3 Hence, the evaluation of force fields in terms of membrane affinity by quantitative comparison to experiment is crucial. An important example is the interaction of NaCl and KCl with a palmitoyl-oleoylphosphatidylcholine (POPC) bilayer in aqueous solution. Lipids containing zwitterionic phosphatidylcholine (PC) headgroups are the most common lipids in eukaryotic cells. In particular, POPC is highly abundant. Biological membranes are surrounded by aqueous solution that contains the cations sodium and potassium and the anion chloride. Sodium is the most abundant cation in the extracellular medium where it occurs in concentrations of 145 mM. Potassium is the most common cation in the intracellular medium where it is found at concentrations of r 2011 American Chemical Society

100150 mM. Chloride occurs at concentrations of 117 mM in the extracellular medium and concentrations of 430 mM in the intracellular medium. The Na+/K+ asymmetry between the extraor intracellular medium is maintained by ion pumps at significant energy cost; this indicates that the difference between Na+ and K+ is biologically highly relevant. This relevance has been attributed to ion-specific interactions with proteins,4,5 but lipids could also play a role. Interactions of ions with lipid membranes can be monitored from the electrophoretic mobilities of liposomes in solutions of varying ion concentrations. Eisenberg et al. studied the interaction of NaCl with vesicles formed from mixtures of PC and negatively charged phosphatidylserine (PS). Although the addition of NaCl rendered electrophoretic mobilities of the vesicles less negative if significant fractions of PS were present, no effect could be resolved for pure PC.6 Tatulian, on the other hand, did see an effect of NaCl on the electrophoretic mobilities of PC vesicles in the presence of divalent cations.7 Klasczyk et al. observed that addition of NaCl rendered the electrophoretic mobilities of pure PC vesicles in the absence of divalent cations more positive.8 The effects seen for monocomponent PC vesicles were more than an order of magnitude smaller than those observed for PC/PS mixtures and required a higher resolution Received: March 29, 2011 Revised: July 28, 2011 Published: August 22, 2011 10587

dx.doi.org/10.1021/jp202928u | J. Phys. Chem. A 2011, 115, 10587–10595

The Journal of Physical Chemistry A

ARTICLE

Table 1. Force Field Parameters for Na+ and K+ Ions Used in This Worka ion

force

ε

field

(kJ/mol)

C6

C12 (kJ nm12/mol)

0.258 7.31425 105 0.245 2.76825 104

2.15718 108 5.98691 108

GROMOS21 0.566 104 0.645 1.63787 105

1.18384 106

Na+ GROMOS21 0.062 KBFF27 0.302 K+

σ

(nm) (kJ nm6/mol)

KBFF

26

22

Dang

2.185 0.4184

3

0.303 6.81362 10

5.30991 106

4

2.67955 107

0.305 3.34832 10

KBFF denotes a force field derived from a KirkwoodBuff analysis of aqueous NaCl or KCl solution. a

of the experimental setup, this is probably why they were discovered much later when the technical equipment was more advanced. MD simulations2,914 indicate that PC headgroups bind Na+ or K+ whereas Cl forms a diffuse counterion layer in the water phase. For POPC, sodium is found to bind at the ester oxygens.9 The simulations indicate that the coordination of sodium by multiple lipids leads to a lateral condensation of the bilayer, implying an increase in the membrane thickness and reduction of the area per lipid.2,9 The number of lipids bound to a sodium ion varies between 1 and 4; i.e., ionlipid complexes with different stochiometries coexist in a bilayer.11 Addition of sodium fluoride to the subphase of dipalmitoylphosphatidylcholine (DPPC) monolayers in the liquid expanded phase at constant area per DPPC molecule leads to an increase in the surface pressure of the monolayers.15 In agreement with the coexistence of Na+lipid complexes with different stochiometries observed in MD simulations, describing the surface pressure increment as a function of NaF concentration in the concentration range between 0 and 1 M assuming a 1:1, 1:2, or 1:3 ion lipid stochiometry yields poor fits. Better fits are obtained using a stepwise binding model considering ionlipid complexes with increasing numbers of lipids, which corresponds more closely to the idea of a mixed stochiometry.15 Nevertheless, data for the heat release arising from the binding of sodium or potassium ions to POPC for NaCl or KCl concentrations between 0 and 0.1 M from recent isothermal calorimetry (ITC) experiments can be fitted well assuming a 1:1 ionlipid stochiometry, probably due to the more narrow concentration range. In these experiments, the membrane affinity of NaCl and KCl in terms of the ratio of the number of bound ions and the number of lipids at given bulk concentration of the ions appears to be the same within the error.8 Rapid solution exchange on a solid supported membrane, on the other hand, indicates that PC interacts more strongly with Na+ than with K+.16 Electrophoresis experiments suggest that the number of cations per area POPC at given NaCl or KCl bulk concentration is larger for Na+ than for K+.8 In MD simulations using various force fields for the lipids, the water, and the ions, consistently show that Cl does not specifically interact with POPC and only forms a diffuse counterion layer if the corresponding cation is adsorbed at the membrane; however, the affinity of the cations is found to depend on the force field.2,914,17 Here, in any case, the interaction is predicted to be weaker for K+ than for Na+.2,1114,17 Previous simulations of the interaction of NaCl or KCl with POPC were conducted using the Berger force field for the lipids,18 the SPC model for the water,19,20 and GROMOS parameters for NaCl.21 The

potassium, however, was treated using different force fields, including GROMOS and parameters by Dang et al.22 The parameters for GROMOS-Na+, GROMOS-K+, and Dang-K+ are given in Table 1. In this paper, we compare affinities of Na+ and K+ for POPC membranes in MD simulations with the values from the ITC experiments. Here we employ the Berger force field for the lipids, the SPC model for the water, and GROMOS for NaCl or KCl, as well as the parameters by Dang et al. for K+. Four other parameter sets for K+ in SPC water are available in the literature.2326 One of these parameter sets was optimized to match solution activities of aqueous KCl solutions.26 As this K+ force field was constructed using KirkwoodBuff theory, it shall be denoted as KBFF-K+. As an accurate description of cosolvent activities is claimed27 to be particularly important for cosolventsolute interactions, we chose this force field as another example for our study. In addition, we employed parameters for NaCl matching solution activities in conjunction with SPC/E water.27 These ion parameters were also derived using KirkwooodBuff theory; hence, this ion force field shall be called KBFF-Na+. The parameters for KBFF-K+ and KBFFNa+ are also listed in Table 1.

2. APPARENT AND INTRINSIC BINDING CONSTANTS In the ITC experiments,8 NaCl or KCl solutions at 500 mM concentration (titrant) were injected into a solution containing POPC vesicles, and the heat release per injection was measured as a function of NaCl or KCl concentration. Hence, a so-called apparent binding constant Kapp for a 1:1 ionlipid stochiometry defined as Kapp ¼

cblip cion ðclip cblip Þ

ð1Þ

was determined via a fitting procedure. Here, clip denotes the total concentration of lipids available for binding, cblip is the concentration of bound lipids bound to an ion, and cion is the total ion concentration. The assumption of a 1:1 ionlipid stochiometry differs from the mixed stochiometry indicated from MD simulations11 and experimental data for expanded PC monolayers for a concentration range of 01 M.15 However, the assumed 1:1 ionlipid stochiometry yields a good fit to the ITC data for the more narrow concentration range of 0100 mM and, hence, is expected to yield a reasonable measure for the effective binding affinity in this range of concentrations. The concentration of lipids available for binding was determined considering only the outer leaflet of the lipid vesicles, corresponding to 55% of the total lipid concentration calculated for a 56 nm vesicle radius as measured from dynamic light scattering. The apparent binding constants for the interaction of ions with POPC obtained from that work were 1.25 ( 0.05 M1 for NaCl and 1.17 ( 0.13 M1 for KCl.8 The apparent binding constant for the interaction of ions with a membrane arises from two effects. First, if the membrane is charged, the concentration of unbound ions at the surface will be different from that in the bulk; the concentration of counterions will be enhanced compared to the bulk concentration whereas co-ions will be depleted at the surface. This effect only depends on the valency of the ions in the solution and is described by the double layer theory.28 Second, local interactions determine an equilibrium between unbound ions at the surface and surface-bound ions. 10588

dx.doi.org/10.1021/jp202928u |J. Phys. Chem. A 2011, 115, 10587–10595


ARTICLE

This equilibrium is ion-specific and determined by the intrinsic binding constant of the ions, Kint, defined as Kint ¼

cblip cion, s ðclip cblip Þ

ð2Þ

Here, cion,s denotes the ion concentration at the surface of the membrane. In this terminology, we will call an ion to be surface active if Kint > 0 and surface-inactive if Kint = 0. In the case of NaCl and KCl solutions, Na+ and K+ are considered to be surface active and Cl is considered to be surface-inactive, as suggested by simulation studies.2,914 The surface concentration of a given ion may be estimated from the ζ potential inferred from electrophoresis experiments. The ζ potential is the electrostatic potential at the shear plane, which is considered to be located at most 0.2 nm away from the membrane surface.6 If ζ is the experimental ζ potential, the surface concentration cion,s of cations with valency 1 is cion, s ¼ cion expð eζ=kB TÞ

ð3Þ

Here, kB denotes the Boltzmann constant, T is the absolute temperature, and cion is ion concentration in the bulk, which in experiments is approximately equal to the total ion concentration. Experimentally, the intrinsic binding constant is determined by fitting eq 2 to the released heats, yielding an effective intrinsic binding constant for the range of surface concentrations of the surface-active ions considered in this fit. Klasczyk et al., though, did not provide intrinsic binding constants for Na+ and K+ interacting with POPC.8 The ζ potentials of POPC vesicles in the absence of alkali chlorides at neutral pH are found to be negative.8 The ζ potentials of POPC vesicles become less negative with decreasing pH and change sign at pH 4.29 Possible causes of the negative ζ potentials of PC could be adsorption of hydroxide29 or impurities.30 Addition of NaCl or KCl renders the ζ potentials less negative but does not invert their sign.8 The surface charge of POPC vesicles in experiment is thus always negative. In contrast, POPC bilayers in pure water in MD simulations have zero surface charge, and addition of NaCl leads to a positive surface charge.

3. METHODS 3.1. Simulations. POPC bilayers in sodium or potassium chloride solutions were studied under periodic boundary conditions using classical molecular dynamics simulations. Initial configurations of the POPC bilayers were obtained from freely available coordinates of a bilayer consisting of 128 POPC lipids by Tieleman et al.31 This system was energy minimized using 50 000 steps of steepest descent first with flexible and, subsequently, constrained bond lengths. In the next step the bilayer was hydrated such that the final bilayer was surrounded by 50605120 water molecules depending on the ion concentration. Ten or thirty sodium or potassium chloride ion pairs were added by replacing water molecules at energetically favorable positions. This was followed by an energy minimization using 50 000 steps of steepest descent, first with flexible and, subsequently, fixed bond lengths. The water was described using the rigid simple point charge (SPC) model.19,20,32 The force field for POPC lipids was taken from Berger et al.18 where hydrogens in hydrocarbon groups are treated implicitly using united atoms. Different force fields for Na+ and K+ were compared. Both ion species were treated using GROMOS8721 or parameters matching solution activities using

KirkwoodBuff theory (KBFF).26,27 Alternatively, potassium was also described using parameters by Dang et al.22 When Na+ and K+ were described with GROMOS or the parameters by Dang et al., Cl was treated using GROMOS. When Na+ and K+ were described with KBFF, Cl was also treated using KBFF.27 The interaction between GROMOS ions and the water or the lipids, as well as the interaction between KBFF-Na+, KBFF-K+, or KBFF-Cl with the lipids, was described using geometrical combination rules for the C6 and C12 parameters of the Lennard-Jones interactions, corresponding to the combination rules of the Berger force field. The interaction between Dang-K+ and the water or the lipids was treated using the modified combination rules by Vacha, Berkowitz, and Jungwirth (VBJ) who also combined DangK+ with GROMOS-Cl, SPC water, and Berger lipids. This force field was used by VBJ to model a cell plasma membrane with an asymmetric multicomponent composition.11 We considered this force field in our work to validate the previous findings by VBJ. Each system was simulated for 100 ns at 300 K and 1 bar. First the simulations using the GROMOS force field for the ions were conducted. The final configurations of the corresponding systems were used for subsequent simulations with other ion force fields. The water molecules were kept rigid using SETTLE,33 and the lengths of covalent bonds in POPC were constrained using LINCS.34 This facilitated the usage of a 2 fs time step for the integration of the equations of motion. The weak coupling technique by Berendsen et al. was used to control the temperature or pressure with relaxation times of 0.1 and 1.0 ps, respectively.35 A pressure of 1 bar normal and lateral to the bilayer corresponding to a tension-free membrane were maintained by semiisotropic scaling of the simulation box. Lennard-Jones interactions were treated using a 1 nm cutoff. Full electrostatic interactions were evaluated using the particle mesh Ewald technique36 with a 1 nm cutoff for direct space, a 0.12 nm grid spacing, and an interpolation order of 4. Nonbonded neighbor lists were updated every 5 steps. All simulations were performed using the GROMACS package, version 3.37 Configurations were saved every 10 ps for further analysis. 3.2. Analysis. The relaxation of the system was monitored from the time evolution of the potential energy. Depending on the cation chloride concentration the relaxation required 120 ns. This period was omitted in the further analysis. This way of estimating the relaxation time of the system was appropriate, as indicated by the fact that the errors of all observables presented here are small and block averages did not show drifts. Number densities for subcomponents of the system as a function of the position normal to the bilayer, z, were obtained by dividing the box in the z direction into 200 slices and smoothening these profiles by convoluting the profiles with a Gaussian distribution with a fwhm of about 1 nm. The membrane thickness was determined from the distance between the maxima of the nitrogen distributions of the two leaflets. The average number of water molecules or ions in contact with phosphate or choline groups, or ester oxygens, and the number of water molecules in contact with a given ion were determined by considering the corresponding heavy atoms with a distance below a cutoff value. The latter was chosen as the position of the first minimum of the corresponding radial distribution function. The numbers of bound cations were obtained by analyzing the minimum distance between individual cations and the POPC membrane using a cutoff value. The latter was chosen as the position of the first minimum of the corresponding minimum distance distributions. 10589



ARTICLE

3.2.1. Apparent Binding Constants. To determine the apparent binding constant Kapp for the interaction of NaCl or KCl with POPC from the simulations, eq 1 was modified by dividing both the numerator and the enumerator of the left side by the total lipid concentration clip and denoting the fraction of lipids binding an ion by α, i.e., α¼

cblip clip

¼

b Nlip

Nlip

¼

b Nion Nlip

ð4Þ

Here, Nlip denotes the number of lipids available for binding. As both leaflets are available for binding in our simulations, Nlip is the total number of lipids in the system. The symbol Nblip denotes the number of lipids bound to an ion and Nbion denotes the number of ions bound to a lipid. Formally, assuming a 1:1 stochiometry, these two numbers are the same. In our analysis, α was determined from Nlip and Nbion. Hence, the apparent binding constant was estimated from Kapp ¼

α cð1 αÞ

ð5Þ

Here, c cion was taken as the effective ion concentration in the bulk. In the experiments, the bulk concentration of the ions is approximately equal to the total ion concentration (which is approximately equal to the concentration of unbound ions), because the bulk is large compared to the interfacial region. This can be deduced from the lipid concentration used in the experiments, 2 mM, which corresponds to almost 28 000 water molecules per lipid.8 The simulations contain only about 40 water molecules per lipid. The bulk concentration of the ions is not directly seen in the simulations due to the small size of the simulation box. The effective bulk concentration is determined by modification of an approach by B€ockmann et al.9 In that earlier work, the ion distributions in the whole membrane/water system were fitted to a PoissonBoltzmann distribution using assumptions about the dielectric permittivity as a function of the distance from the membrane center. As the solvation free energy of the ions was not considered in this approach, the ion densities in the lipid tail region were unphysically high. To avoid this artifact and assumptions on the dielectric permittivity profile at the water/membrane interface, we only considered the ion concentrations in the water, i.e., where the water density was equal to the bulk density. In this case, the ion distribution can be safely assumed to be solely governed by the electrostatic potential without effects from gradients in chemical potential. In this region, the concentration of cations c+(z) and anions c(z) at the distance z from the bilayer at given bulk concentration c0 are considered to obey the equation c( ðzÞ ¼ c0 exp½ q( ΦðzÞ=kB T

ð6Þ

where q( = (e is the charge of the ion with e denoting the elementary charge. Furthermore, Φ(z) is the electrostatic potential at distance z, kB is the Boltzmann constant, and T is the absolute temperature. Equation 6 yields pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi c0 ¼ cþ ðzÞ c ðzÞ ð7Þ Hence, c0 was determined by applying eq 7 to each xy slice in the water region and averaging over all the slices. The corresponding standard error was obtained from averages by dividing the water region into two halves.

Figure 1. Sodium and potassium adsorption at the POPC bilayer for different ion concentrations and different force fields. Shown are (a) a representative lipid configuration in stick representation, as well as number densities of (b) lipids and chlorides, (c) sodium, and (d) potassium. In (a), colors distinguish between hydrocarbons or carbons (cyan), oxygens (red), nitrogen (blue), and phosphate (tan). All density profiles are axis-symmetric at zero. For sodium and potassium the upper curve of one color represents the system with 30 cations and the lower curve of one color the system with 10 cations. The densities of the lipids are scaled for better comparison with the ion densities.

Errors for α, numbers of next neighbors and contacts, areas per lipid, and membrane thickness give the standard error from block averages obtained by dividing the trajectory into five segments. Errors for Kapp were computed using error propagation. 3.2.2. Intrinsic Binding Constants. To determine intrinsic binding constants, the surface concentration cs of cations in the different systems was determined. For most simulations, the cation concentration as a function of the distance from the membrane center showed a maximum at the headgrouptail interface of the bilayers corresponding to the layer of adsorbed ions. Between this maximum and the water region the cation concentration showed a minimum, and cs was taken as the concentration at this minimum. The standard error for cs,MD was obtained from the values for the two watermembrane interfaces. Hence, the corresponding intrinsic binding constant was determined from Kint ¼

α cs, MD ð1 αÞ

ð8Þ

To validate the values for Kint from the simulations, experimental intrinsic binding constants for the surface concentrations in the simulations cs with corresponding standard errors were determined for the bulk concentrations cexp,1, cexp,2, ..., cexp,n, with cexp,1 10590



ARTICLE

Table 2. Number of Choline Groups (N), Phosphate Groups (P), or Ester Oxygens (O) of the POPC Lipid in Contact with a Given Cationa no. of cations

ion-O

ion-P

ion-N

Na GROMOS 10

2.19 ( 0.25

0.0837 ( 0.0094

0.435 ( 0.047

30

1.82 ( 0.12

0.0329 ( 0.0053

0.477 ( 0.023

10

3.14 ( 0.14

0.150 ( 0.039

0.252 ( 0.042

30

2.170 ( 0.038

0.119 ( 0.013

0.313 ( 0.011

KBFF

K GROMOS 10

0.164 ( 0.020

0.157 ( 0.019

0.271 ( 0.011

30

0.085 ( 0.012

0.1015 ( 0.0068

0.2627 ( 0.0046

10

1.14 ( 0.14

1.39 ( 0.12

0.807 ( 0.071

30

1.019 ( 0.090

1.17 ( 0.11

0.959 ( 0.045

KBFF

a In each system 128 POPC lipids are surrounded by 50605120 water molecules depending on the ion concentration.

< cexp,2 < 3 3 3 < cexp,n, considered in the experiments. If there was a k for which cs,exp,k e cs,MD e cs,exp,k+1, the intrinsic binding constants for cs,exp,k and cs,exp,k+1 were calculated via eq 8. Here, cs, MD was replaced by cs,exp,k and cs,exp,k+1, and α was substituted by αk and αk+1. The latter were calculated from αi ¼

Kapp ci 1 þ Kapp ci

ð9Þ

with i = k and i = k + 1, correspondingly, and Kapp denoting the apparent binding constant from experiment. Hence, intrinsic binding constants Kint,1 and Kint,2 were obtained. The standard error for Kint,i, se(Kint,i), i = 1,2, was evaluated via error propagation. The experimental intrinsic binding constant for the surface concentration cs,MD was thus considered to be in the range [Kint,1 se(Kint,1), Kint,2 + se(Kint,2)].

4. RESULTS AND DISCUSSION Molecular dynamics simulations have been performed to study the interaction of potassium and sodium with POPC bilayers comparing different force fields for the ions. This section consists of two parts. In part one, the microscopic details of the interaction of NaCl and KCl with POPC for KBFF are compared with the properties from simulations with the GROMOS force field. In the second part, the membrane affinity of Na+ and K+ for GROMOS, KBFF, and Dang are validated against data from isothermal calorimetry and electrophoresis experiments. The systems considered comprised a bilayer containing 128 lipids, and 50005120 water molecules depending on the number of ions, with 10 or 30 ion pairs or without ions. 4.1. Microscopic Details of Ion Adsorption. Figure 1 shows the number density profiles of the lipids and the ions. GROMOSNa+, KBFF-Na+, and KBFF-K+ are strongly adsorbed in the headgroup region whereas Cl forms a diffuse layer in the water close to the membrane. The number densities of the adsorbed cations increase in the order KBFF-K+, GROMOS-Na+, and

Figure 2. Effect of ions on membrane structure. Shown are (a) the membrane thickness, (b) the area per lipid, and (c) the volume per lipid for a hydrated POPC bilayer in the presence of NaCl or KCl as a function of the effective ion concentration.

KBFF-Na+. In contrast, GROMOS-K+ is depleted from the headgroup region and the corresponding Cl distribution is more uniform. In general, the surface densities tend to be higher for Na+ than for K+. This means that the positive surface charge of the membrane from the adsorption of cations is higher for Na+ than for K+. Experimentally, the surface charge can be probed by means of electrophoresis experiments. These experiments indicate that POPC vesicles in the absence of alkali metal chlorides carry a weak negative charge. Addition of NaCl or KCl renders the surface charge less negative, the effect being larger for NaCl than for KCl, indicating that at given alkali metal chloride concentration the surface density of Na+ is higher than for K+, in agreement with our results. At the higher concentration the cation density profiles exhibit a minimum between the bulk water and the headgroup region. This indicates that the adsorbed cations repel unbound cations in the water region, an effect that is not significant at the lower ion concentration. To study the binding sites of the ions, the average number of ester oxygens, as well as phosphate and choline groups, in contact with a given Na+ or K+, averaged over all ions in the system, were 10591



ARTICLE

Table 3. PoissonBoltzmann Derived Effective Bulk Concentration (conc) and Number of Bound Cations per POPC Lipid from MD Simulations no. of cations

conc [mol/L]

no. of bound cations per lipid

Table 4. Apparent Binding Constant, Kapp, for Na+ and K+ from the Number of Bound Lipids at the Lowest Ion Concentration in Table 3 According to eq 5, as Well as Corresponding Experimental Values8

Na-GROMOS 10

0.0453 ( 0.0058

0.0579 ( 0.0062

30

0.207 ( 0.010

0.1500 ( 0.0024

2.28

experiment

1.25 ( 0.05

GROMOS

0.083 ( 0.020

0.0711 ( 0.0030

KBFF

0.1593 ( 0.0024

1.32 ( 0.29

Dang

0.63 ( 0.12

experiment

1.17 ( 0.13

30

0.190 ( 0.020

10

0.1244 ( 0.0076

0.0102 ( 0.0024

30

0.389 ( 0.016

0.0211 ( 0.0070

10

0.049 ( 0.011

0.0602 ( 0.0016

30

0.203 ( 0.019

0.152 ( 0.046

K-GROMOS

K-KBFF

K-Dang 10

0.0768 ( 0.0088

0.0359 ( 0.0094

30

0.275 ( 0.015

0.051 ( 0.028

determined for the two ion concentrations. The results are shown in Table 2. At the lowest ion concentration, GROMOS-Na+ is coordinated by about two and KBFF-Na+ by three ester oxygens on average. The average coordination number for GROMOS-Na+ corresponds to the most abundant binding ionlipid stochiometry of 1:2 found previously, though 1:1, 1:3, and 1:4 complexes were also observed.11 Both ions barely contact the phosphate groups. GROMOS-Na+ is in contact with choline half of the time and KBFF-Na+ with choline groups one-quarter of the time. GROMOS-K+ barely contacts ester oxygens or phosphate groups and forms more contacts with choline. This correlates with the density profiles where a gradual decay of the K+ concentrations is observed when the membrane is approached from the water phase. KBFF-K+ binds on average one ester oxygen and one phosphate group. KBFF-K+ also forms twice as many contacts with choline than Na+. The binding mode of KBFF-K+ and KBFF- as well as GROMOS-Na+ in terms of coordination numbers observed here correlates well with the position of the corresponding peak of the electron density profile, which resides less deeply in the headgroup region for K+ than for Na+. The adsorption of cations leads to a decrease in the area per lipid, as shown in Figure 2a. Although the effect of GROMOS-K+ is very small and the same for 10 and 30 ion pairs, a significant effect, increasing with ion concentration, is observed for the other systems. The magnitude of the effect increases in the order KBFF-K+, GROMOS-Na+, and KBFF-Na+. The magnitude of the effect thus correlates with the surface concentrations of the ions shown in Figure 1c,d. Comparison with Figure 2b,c shows that the reduction in area per lipid correlates with an increase in membrane thickness such that the volume per lipid remains largely constant. 4.2. Apparent Binding Constants. To validate the membrane affinity of the ions for different force fields against experimental data, the number of bound ions per lipid and the effective bulk concentration of the ions was determined to estimate apparent binding constants for the ions. The systems considered were those containing 10 or 30 ion pairs, from the simulations where sodium

and potassium were described using GROMOS or KBFF. In addition, simulations with 10 and 30 KCl ion pairs were conducted using ion parameters by Dang et al. The effective concentrations for the ions were evaluated by assuming that the ion distribution in the water is related to the electrostatic potential according to eq 6 such that the local bulk concentration is the geometric mean of the anion and cation concentration as in eq 7, averaged over the water layer. The effective bulk concentrations for NaCl and KCl for 10 or 30 ion pairs and the corresponding ratios α (number of bound ions)/ (number of lipids) are shown in Table 3. The effective concentrations for NaCl treated with GROMOS for 10 and 30 ion pairs were 0.0453((0.0058) and 0.207((0.010) M, respectively. B€ockmann et al. used the same simulation systems, force field, and time scales, but a slightly different method to determine the effective bulk concentration of the ions. Namely, they fitted the ion distribution from the solution of the PoissonBoltzmann equation for given charge distribution from the lipid and an assumed dielectric permittivity profile as a function of the distance from the membrane surface. In this approach, the distribution of the ions was assumed to be determined by the electrostatic potential only, even in the membrane region, neglecting the solvation free energy, which led to unrealistically high ion concentrations in the tail region of the membrane. We believe that our approach is more appropriate; nevertheless, it is useful to compare the outcomes of the two methods. The effective bulk concentration for the system containing 10 and 30 ion pairs obtained by B€ockmann et al. was 0.050 ( 0.030 and 0.220 ( 0.030 mol/L,9 respectively. Hence, B€ockmann’s estimates are associated with a higher statistical uncertainty. However, within the error, both methods yield the same effective ion concentration for these particular systems. As argued above, the apparent binding constants for NaCl or KCl interacting with POPC determined experimentally are effective binding constants for NaCl or KCl concentrations of at most cmax = 0.1 M. To be comparable to the experiments, the apparent binding constants from the simulations should be evaluated from systems with effective ion concentrations below cmax. As seen from Table 3, the effective ion concentration was below 0.1 M for most systems containing 10 ion pairs. An exception was the system for which KCl was treated using GROMOS; here, the effective ion concentration was 0.1244((0.0076) M and, hence, only slightly above the range to which the experimental binding constants referred. Hence, we used the systems containing 10 ion pairs to calculate the corresponding apparent binding constants, Kapp. The binding constants calculated according to eq 5 are shown in Table 4. 10592



ARTICLE

Table 5. ζ Potentials for POPC Vesicles in the Presence of NaCl and KCl at Concentrations c from Experiments8

Table 6. Intrinsic Binding Constants Kint for Na+ and K+ Interacting with POPC from Simulations and Experiments8a

ion

c (mM)

ζ (mV)

ion

Na+

0.050

1.75 ( 0.30

K+

0.100 0.035

1.14 ( 0.36 5.33 ( 0.40

0.100

3.22 ( 0.27

The binding constant for Na+ when treated with the GROMOS force field is Kapp,Na,GROMOS = 1.36((0.22) M1, which within the error agrees with the experimental value Kapp,Na,exp = 1.25((0.05) M1. The apparent binding constant for KBFF-Na+ is at least 2.28 M1 and, hence, exceeds the experimental value. The binding constants for Na+ are therefore ranked as follows: Kapp, Na, GROMOS ≈ Kapp, Na, exp < Kapp, Na, KBFF The binding constant for GROMOS-K+ is Kapp,K,GROMOS = 0.083((0.020) M1 and, hence, is more than an order of magnitude smaller than the experimental value, Kapp,K,exp = 1.17((0.13) M. The binding constant for Dang-K+ is Kapp,K,Dang = 0.63((0.12) M1 and, thus, is about half of the experimental value. The binding constant for KBFF-K+ is Kapp,K,KBFF = 1.32((0.29) M1 and, therefore, agrees with the experimental value within the error. Thus, the binding constants for K+ are ranked according to Kapp, K, GROMOS , Kapp, K, Dang < Kapp, K, KBFF ≈ Kapp, K, exp Hence, GROMOS-Na+ and KBFF-K+ yield the same apparent binding constants as the ITC experiments. That means that for given bulk concentration these force fields reproduce the number of cations per lipid indicated from the experimental data. In contrast to the ITC experiments that probe the ratio of the number of bound ions and the number of lipids, the electrophoresis experiments are sensitive to the surface charge, i.e., the number of ions per surface area. Although the binding constants for GROMOS-Na+ and KBFF-K+ are the same within the error (and reproducing the experimental values), the number densities for bound ions are higher for GROMOS-Na+ than for KBFF-K+, as shown in Figure 1c,d and discussed above. This is because, at the same ion concentration, Na+ induces a larger reduction in the area per lipid than K+. 4.3. Intrinsic Binding Constants. As the POPC bilayers in the presence of NaCl or KCl were positively charged in the simulations but negatively charged in the experiments, agreement of the apparent binding constants between the simulations and the experiments may not mean agreement of the intrinsic binding constants. The intrinsic binding constants for the simulations and the experiments were evaluated as described in the Methods. Comparison with experiments with our approach required us to use systems for which the surface concentrations of the ions were in the range of finite surface concentrations also considered in experiments. This was the case for 10 ion pairs with GROMOSK+ and for 30 ion pairs with the other ion force fields. Experimental surface concentrations were estimated from the ζ potentials given in Table 5. The surface concentration was taken as the cation concentration at the minimum between the peak from the membrane adsorbed cations and the distribution of cations in the water also seen in Figure 1. GROMOS-K+ did not show a maximum in the headgroup region, as seen in Figure 1d and discussed above; here,

N

Na+

GROMOS

30

0.0675 ( 0.0018

2.614 ( 0.036

K+

experiment GROMOS

10

0.0530.106 0.0433 ( 0.0053

1.1041.268 0.238 ( 0.05

KBFF

30

0.0683 ( 0.0050

2.62 ( 0.67

Dang

30

0.100 ( 0.015

1.30 ( 0.21

0.0420.114

0.9571.253

experiment

cs (M)

Kint (M1)

method

a

The intrinsic binding constants are effective values for the given range of surface concentrations, cs. For the simulations, the number of ion pairs N in the corresponding simulation systems is given.

the surface concentration was chosen as the concentration at the kink between the density profile for K+ in the membrane and that in solution, visible in Figure 1d. KBFF-Na+ was not considered in this analysis, as this force field was observed to overestimate the apparent binding constant of Na+ and, hence, should also overestimate the intrinsic binding constant of Na+. The results are shown in Table 6. The intrinsic binding constant for GROMOS-Na+ for a surface concentration of cs = 0.0675 ((0.0018) M was 2.614((0.036) M1, exceeding the experimental range for Kint (1.1041.268 M1) by a factor of 2.5. GROMOS-K+ showed a surface concentration of 0.0433((0.0053) M, yielding an intrinsic binding constant of 0.238((0.05) M being lower than the experimental range for Kint (0.9571.202 M1) by a factor of 4. The intrinsic binding constant for KBFF-K+ for cs = 0.0683((0.0050) M was 2.62((0.67) M1, exceeding the experimental range for Kint (0.9571.253 M1) by a factor of 2.5. DangK+ yielded an intrinsic binding constant of 1.30((0.21) M1 for cs = 0.100((0.015) M, in agreement with the experimental range for Kint considering the statistical error.

5. DISCUSSION We have compared apparent and intrinsic binding constants for the interaction of Na+ and K+ with POPC bilayers from MD simulations using different ion force fields with values from isothermal calorimetry. GROMOS-Na+ yields the same apparent binding constant as the experiments, whereas KBFF-Na+ overestimates Kapp by a factor of 2. GROMOS-K+ yields an apparent binding constant that is 1 order of magnitude lower than the experimental value. The apparent binding constant for Dang-K+ is a factor of 2 lower than the experimental value, whereas KBFFK+ reproduces the experimental binding constant. Agreement in terms of apparent binding constant means that for a given alkali metal chloride concentration the ratio of the number of bound ions and the number of lipids in the simulations is the same as in the experiments. Although the apparent binding constants for GROMOS-Na+ and KBFF-K+ are the same within the error, the surface density of membrane-bound GROMOS-Na+ at a given alkali metal chloride concentration is higher than that of membrane-bound KBFF-K+, because GROMOS-Na+ induces a stronger reduction in the area per lipid than KBFF-K+. This correlates well with the fact that for given alkali metal chloride concentration NaCl renders the electrophoretic mobilities of POPC more strongly positive than KCl.8 Although GROMOS-Na+ and KBFF-K+ reproduce the apparent binding constants, they overestimate the intrinsic binding constants by a factor of 3. GROMOS-K+ underestimates Kint by a factor of 4. Among the force fields considered, only 10593



ARTICLE

Dang-K+ reproduces the intrinsic binding constant from experiment. However, as mentioned above, its apparent binding constant is a factor of 2 lower than the experimental value. On the other hand, even though GROMOS-Na+ and KBFF-K+ yield the same surface coverage of POPC by ions as in experiments at neutral pH, the underlying mechanism is different. Strictly spoken, MD simulations of POPC bilayers, although meant to model conditions at neutral pH, may in fact mimic conditions at pH 4. Consistent modeling of the electrostatics of POPC bilayers at pH 7 would mean to include the negative surface charge indicated from electrophoresis experiments whose origin, though, is not understood. The difference in binding affinities between the ions might be explained in terms of differences in the ionic radii. Data from ITC and electrophoresis experiments show that Li+ binds more strongly to POPC than Na+.8 At the same time, a dehydrated Li+ ion is smaller than a dehydrated Na+ ion. This suggests that the membrane affinity of a monovalent atomic cation increases with decreasing size of the bare cation. One measure for the size of an ion is its van der Waals radius σ, which corresponds to the distance where the Lennard-Jones potential between like particles is zero. This parameter is given in Table 1. The ranking of the σ parameter for a given ion type, σNa, KBFF < σ Na, GROMOS

ð10Þ

and σK, KBFF < σK, Dang , σK, GROMOS

ð11Þ

is indeed opposite to the ranking of the corresponding apparent or intrinsic binding constants; i.e., smaller ions bind more strongly. In particular, GROMOS-K+, which shows a much lower surface affinity than the other ion models exhibits a van der Waals radius twice as large as that of KBFF-K+ and Dang-K+. However, the rule that smaller ions bind more strongly does not strictly apply when we compare force fields for Na+ and K+. Although σ is smaller for KBFF-K+ (0.303 nm) compared to GROMOS-Na+ (0.258 nm), both models show the same binding affinity. The strong binding affinity of KBFF-K+ despite the larger σ might be due to the relatively large depth of the Lennard-Jones potential, ε = 2.185 kJ/mol. This ε parameter is larger than for other K+ force fields; however, it is close to the value ε = 1.5 kJ/mol for parameter set 5a for K+ in SPC/E water parametrized on the basis of hydration free energies by Horinek et al.38 KBFF-K+ and KBFF-Na+ have been parametrized to match solution activities of KCl and NaCl. Solution activities, as also hydration free energies which provide another possible target for force field parametrization, are only accessible in a combined cationanion form. Hence, errors in the Na+ or the K+ parameters may be compensated by errors in the Cl parameters. On the other hand, the affinity of Na+ or K+ for POPC will be very sensitive to the cation parameters but largely insensitive to the Cl parameters. This could be the reason why the parameters matching solution activities are not accurate in terms of intrinsic binding constants for ionPOPC interactions. On the other hand, the membrane affinities of the ions will also depend on the parameters for the lipid headgroups. To not disturb other properties, the membrane affinities of ions for a given force field may be improved by breaking the combination rules for LJ interactions between the cations and the carbonyl or phosphate group oxygens of the lipids.

6. CONCLUSION Our work shows that the intrinsic binding constants for interactions of Na+ and K+ with POPC bilayers in conjunction with the Berger force field for the lipids and the SPC water model are too high for GROMOS-Na+ and KBFF-K+ but too low for GROMOS-K+ compared to experiment. Good agreement with experiment in terms of intrinsic binding constants is obtained for Dang-K+. The strong difference in the membrane affinities of Na+ and K+ with POPC bilayers seen by Vacha et al.11 could arise from an overestimation of the intrinsic binding affinity of Na+ by GROMOS. Although they overestimate the intrinsic binding constants, GROMOS-Na+ and KBFF-K+ reproduce apparent binding constants for POPC at neutral conditions in experiments. In agreement with results from isothermal calorimetry, for these force fields, Na+ and K+ show the same affinity in terms of the number of bound ions per number of lipids. Despite this, the surface density from membrane-bound cations is larger for Na+ than for K+, because Na+ causes a stronger reduction in the area per lipid than K+, resolving an apparent controversy between the results from ITC and electrophoresis experiments.8 The agreement of the apparent binding constants despite the disagreement of the intrinsic binding constants for GROMOSNa+ and KBFF-K+ is due to the difference in the surface charge of POPC between simulation and experiment. Whereas POPC bilayers are negatively charged in experiment at neutral conditions, even in the presence of NaCl or KCl, they are typically neutral (in the absence of ions) or positively charged (in the presence of NaCl or KCl) in MD simulations. The latter, though meant to model neutral pH, may mimic conditions at pH 4.29 To model POPC bilayers under physiological conditions reliably, future simulation studies should consider the negative surface charge of POPC indicated from the electrophoresis experiments (with the challenge that its origin is not unclear). Hence, a force field reproducing intrinsic binding constants for Na+ should be identified or developed. Such a force field may imply breaking the mixing rules for the interactions with the carbonyl oxygens of the lipids. On the other hand, more experiments on POPC should be conducted at pH 4 where the surface charge of pure POPC is zero, facilitating comparison with simulation studies. Similar arguments concerning the difference in surface charge between simulations and experiments and the need to validate simulations in terms of the membrane affinity of ions refer to other zwitterionic lipids (PC lipids with other tails as well as to phosphatidylethanolamine lipids).39,40 A framework for the validation of the membrane affinity of ions in simulations is provided by the present work. ’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ REFERENCES (1) McLaughlin, S. Annu. Rev. Biophys. Biophys. Chem. 1989, 18, 113–136. (2) Gurtovenko, A. A.; Vattulainen, I. J. Phys. Chem. B 2008, 112, 1953. (3) Cordomi, A.; Edholm, O.; Perez, J. J. J. Chem. Theory Comput. 2009, 5, 2125–2134. (4) Vrbka, L.; Vondrasek, J.; Jagoda-Cwiklik, B.; Vacha, R.; Jungwirth, P. Proc. Natl. Acad. Sci. U. S. A. 2006, 103, 15440–15444. 10594



ARTICLE

(5) Hess, B.; van der Vegt, N. F. A. Proc. Natl. Acad. Sci. U. S. A. 2009, 106, 13296–13300. (6) Eisenberg, M.; Gresalfi, T.; Riccio, T.; McLaughlin, S. Biochemistry 1979, 18, 5213–5223. (7) Tatulian, S. Eur. J. Biochem. 1987, 170, 413–420. (8) Klasczyk, B.; Knecht, V.; Lipowsky, R.; Dimova, R. Langmuir 2010, 26, 18951–18958. (9) B€ockmann, R. A.; Hac, A.; Heimburg, T.; Grubm€uller, H. Biophys. J. 2003, 85, 1647–1655. (10) Pandit, S. A.; Bostick, D.; Berkowitz, M. L. Biophys. J. 2003, 84, 3743–3750. (11) Vacha, R.; Berkowitz, M. L.; Jungwirth, P. Biophys. J. 2009, 96, 4493–4501. (12) Vacha, R.; Siu, S. W. I.; Petrov, M.; Boeckmann, R. A.; BaruchaKraszewska, J.; Jurkiewicz, P.; Hof, M.; Berkowitz, M. L.; Jungwirth, P. qJ. Phys. Chem. A 2009, 113, 7235–7243. (13) Vacha, R.; Jurkiewicz, P.; Petrov, M.; Berkowitz, M. L.; B€ockmann, R. A.; Barucha-Kraszewska, J.; Hof, M.; Jungwirth, P. J. Phys. Chem. B 2010, 114, 9504–9509. (14) Cordomi, A.; Edholm, O.; Perez, J. J. J. Phys. Chem. B 2008, 112, 1397–1408. (15) Leontidis, E.; Aroti, A. J. Phys. Chem. B 2009, 113, 1460–1467. (16) Garcia-Celma, J. J.; Hatahet, L.; Kunz, W.; Fendler, K. Langmuir 2007, 23, 10074–10080. (17) Lee, S.-J.; Song, Y.; Baker, N. A. Biophys. J. 2008, 94, 3565–3576. (18) Berger, O.; Edholm, O.; J€ahnig, F. Biophys. J. 1997, 72, 2002–2013. (19) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; D. Reidel Publishing Co.: Dordrecht, The Netherlands, 1981; pp 331342. (20) Jorgensen, W.; Chandrasekhar, J.; Madura, J.; Impey, R.; Klein, M. J. Chem. Phys. 1983, 79, 926–935. (21) van Gunsteren, W. F.; Berendsen, H. J. C. GROMOS; Biomos: Groningen, The Netherlands, 1987. (22) Dang, L. X.; Schenter, G. K.; Glezakou, V.-A.; Fulton, J. L. J. Phys. Chem. B 2006, 110, 23644–23654. (23) Aqvist, J. J. Phys. Chem. 1990, 94, 8021–8024. (24) Straatsma, T.; Berendsen, H. J. C. J. Chem. Phys. 1988, 89, 5876– 5886. (25) Reif, M. M.; Huenenberger, P. H. J. Chem. Phys. 2011, 134, 144104. (26) Klasczyk, B.; Knecht, V. J. Chem. Phys. 2010, 132, 024109. (27) Weerasinghe, S.; Smith, P. E. J. Chem. Phys. 2003, 119, 11342. (28) Hunter, R. J. Zeta Potential in Colloid Science; Academic Press: London, 1981. (29) Zhou, Y.; Raphael, R. M. Biophys. J. 2007, 92, 2451–2462. (30) Pincet, F.; Cribier, S.; Perez, E. Eur. Phys. J. B 1999, 11, 127–130. (31) Tieleman, D. P.; Sansom, M. S. P.; Berendsen, H. J. C. Biophys. J. 1999, 76, 40–49. (32) Hermans, J.; Berendsen, H. J. C.; van Gunsteren, W. F.; Postma, J. P. M. Biopolymers 1984, 23, 1513–1518. (33) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952–962. (34) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463–1472. (35) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (36) Darden, T.; York, D.; Pedersen, L. J. Chem. Phys. 1993, 98, 10089. (37) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701–1718. (38) Horinek, D.; Mamatkulov, S. I.; Netz, R. R. J. Chem. Phys. 2009, 130. (39) Shchelokovskyy, P.; Tristram-Nagle, S.; Dimova, R. New J. Phys. 2011, 13. (40) Roy, M.; Gallardo, M.; Estelrich, J. J. Colloid Interface Sci. 1998, 206, 512–517.

10595


Validating Affinities for IonâLipid Association from Simulation against

Recommend Documents

Validating Affinities for IonâLipid Association from Simulation against