Subscriber access provided by UPSTATE Medical University Health Sciences Library
Article
Quantifying the Relationship Between Curvature and Electric Potential in Lipid Bilayers Dennis Skjøth Bruhn, Michael Andersen Lomholt, and Himanshu Khandelia J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.6b03439 • Publication Date (Web): 10 May 2016 Downloaded from http://pubs.acs.org on May 13, 2016
Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.
The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.
Page 1 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Quantifying the Relationship Between Curvature and Electric Potential in Lipid Bilayers Dennis S. Bruhn, Michael A. Lomholt, and Himanshu Khandelia∗ MEMPHYS - Center for Biomembrane Physics, Department of Physics, Chemistry and Pharmacy, University of Southern Denmark, Campusvej 55, 5230 Odense M, Denmark E-mail:
[email protected] Phone: +45 65 50 35 10
1
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Abstract Cellular membranes mediate vital cellular processes by being subject to curvature and transmembrane electrical potentials. Here we build upon the existing theory for flexoelectricity in liquid crystals to quantify the coupling between lipid bilayer curvature and membrane potentials. Using molecular dynamics simulations, we show that head group dipole moments, the lateral pressure profile across the bilayer and spontaneous curvature all systematically change with increasing membrane potentials. In particular, there is a linear dependence between the bending moment (the product of bending rigidity and spontaneous curvature) and the applied membrane potentials. We show that biologically relevant membrane potentials can induce biologically relevant curvatures corresponding to radii of around 500 nm. The implications of flexoelectricity in lipid bilayers are thus likely to be of considerable consequence both in biology and in model lipid bilayer systems.
Introduction A biological membrane needs to form compartments in order to sustain life. Compartmentalization is achieved by curving the lipid bilayer that makes up the membrane. The composition of the membrane is tuned such that curving it is possible, but resistance towards this exists, so that resulting compartments do not collapse. Apart from this “basic” use of curvature some organelles consist of highly curved segments of membrane, like the endoplasmic reticulum or the Golgi, and such segments are also found in transport processes, like exo- and endocytosis. 1 Transmembrane electric potential differences, or membrane potentials, are essential to life and sustain neuronal excitability, signal transmission, transport, and influence the activity of membrane proteins. 2 The potential difference is ≈70 mV in resting cells, and is mainly generated by ion gradients maintained by a combination of passive ion diffusion across ion channels and active uphill ion transport by ion pumps. However, the curvature of the mem2
ACS Paragon Plus Environment
Page 2 of 22
Page 3 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
brane makes an important, but somewhat under appreciated contribution to the membrane potential. For a membrane under zero electric fields, bending the membrane alters the charge and electric dipole distributions mainly in the head group regions, which creates a potential difference between the leaflets (or conversely, an existing potential can induce a curvature in the membrane). Such effects are described by flexoelectricity. 3 The implications of flexoelectricity in membranes are significant. Cells contain highly curved dynamic membranes, which may induce large electrical potentials locally affecting a myriad range of membrane processes including the activity of membrane proteins. In model systems too, vesicles of different curvatures (sizes) will have different electrical potentials across them. The flexoelectric theory was first developed in the context of liquid crystals and refers to the coupling between electric polarization and mechanical strain in dielectric materials. Although a lot of theoretical and experimental work has been done in the context of liquid crystals, 3,4 nearly all of the flexoelectricity work on lipid bilayers has been pioneered by a single research group. 3,4 Here we carry out extensive Molecular Dynamics (MD) simulations of a simple lipid bilayer model to examine and measure the change in the spontaneous curvature, or more precisely the bending moment obtained by multiplying spontaneous curvature with bending rigidity, upon application of membrane potentials of different magnitudes. We find that the membrane potential induces biologically relevant magnitudes of spontaneous curvature. To our best knowledge, this is the first instance of quantitative investigation of flexoelectricity using simulations. We also find that the lateral pressure profiles changes systematically with increasing membrane potentials. Such changes will influence conformational transitions in, for example, mechanosensitive proteins which respond to changes in the stress of their lipid environments. 5–7 Examples of such proteins are the channels Large- and Small-conductance mechanosensitive channel (MscL/S), some of which are present in neurons 8 that are subject
3
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 4 of 22
to very high action potentials. The conformational landscape of MscL has been shown to be very sensitive to both the overall membrane tension and the lateral pressure profile. 9
Theory and Methods The spontaneous curvature, c0 , of a bilayer is defined as the curvature which minimizes the mean curvature term in the Helfrich free energy per unit area fc : 10,11 1 fc = κ(c1 + c2 − c0 )2 + κ ¯ c1 c 2 , 2
(1)
where κ is the bending modulus, c1 and c2 are the principal radii, and κ ¯ is the saddle-splay modulus. The composition and environment determines c0 . The free energy written on the above form is convenient for studying membranes that are free to adjust their area. If we partly multiply out the square we can rewrite the free energy in a form where each term is proportional to its own elastic coefficient: 1 ¯ c1 c2 . fc = κ(c1 + c2 )2 − β(c1 + c2 ) + σ + κ 2
(2)
Here we have introduced the flat membrane tension σ = κc20 /2 and the bending moment β = κc0 . This second form is more convenient for the situation in this paper, where we study a topologically constrained closed membrane. In particular it allows for the possibility that constraints at the boundaries can alter the tension σ independently from the values of κ and c0 . We use periodic boundary conditions in the simulations described below, which eliminates net curvature. Thus, spontaneous curvature cannot be observed visually, or calculated using, for example, polynomial fitting. 12 Instead we calculate the bending moment by first computing the lateral pressure profile (or simply pressure profile), which, traditionally, is
4
ACS Paragon Plus Environment
Page 5 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
defined as: π(z) = PL (z) − PN (z),
(3)
where PL (z) = (Pxx (z) + Pyy (z))/2 is the pressure in the plane of the bilayer as a function of the height, z, and PN (z) is the pressure normal to the bilayer, along the z-axis. The first moment of π(z) gives the bending moment of the membrane 11 :
β = κc0 =
Z
(z − z0 )π(z)dz,
(4)
where z0 is at the center of the bilayer. The same approach has been used before to calculate bending moment and spontaneous curvature of asymmetric lipid bilayers in dissipative particle dynamics simulations. 13 Note that the pressure tensor, Pij , is the negative stress tensor, Tij : Pij = −Tij .
Local Stress due to External Electric Field In our simulations we use a transmembrane electric field to generate membrane potentials. The interactions between this constant external electric field, E ext , and the particles are not included in the standard implementation 14,15 of the local stress tensor calculation (see Simulation Details). Contributions to the pressure profile from these interactions therefore have to be calculated. ˆ where zˆ is a unit vector perpendicular to the plane The applied field is: E ext = E ext z, of the bilayer. The total electric field in the system is E = E ext + E p , where E p is the field generated by the charged particles:
E p (r) =
1 X qα (r − rα ) , 4πǫ0 α |r − rα |3
(5)
where the sum is over all charged particles in the system, qα and rα are the charge and position of the α’th particle, respectively.
5
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 6 of 22
In the following we will consider the stress tensor definition used by Irving and Kirkwood 16 to calculate contributions due to interactions between the particles and the field. The electrostatic interaction tensor in the Irving-Kirkwood formulation, TijIK , satisfies ∂j TijIK = ρEip , where ∂j = ∂/∂rj , ρ(r) =
P
α qα δ(r
(6)
− rα ) and Eip is the i’th component of E p . This shows
that TijIK is not the “full” stress tensor Tij , because this has to account for the full electric force on the particles, i.e., it has to satisfy:
∂j Tij = ρEi ,
(7)
where Ei is the i’th component of the total electric field. Hence we can split Tij into Tij = TijIK + Tijint where the interaction tensor must satisfy: ∂j Tijint = ρEiext .
(8)
A possible choice for Tijint is: Tijint = E ext δiz δjz
X
qα θ(z − zα )δ(x − xα )δ(y − yα ),
(9)
α
since ∂z θ(z − zα ) = δ(z − zα ) makes sure that this choice satisfies eq. (8). Here θ(z) is the Heaviside step function and the sum includes all charged particles. Note that one could have chosen another Tijint as long as it satisfies eq. (8). However, this choice for Tijint can be interpreted as the Irving-Kirkwood stress tensor contribution resulting from the membrane system’s interaction with two collections of opposite charges placed far away along the z-axis above/below the membrane, such that the electric field is constant and equal to E ext locally around the membrane. The direct lines of force between the two collections and the charges
6
ACS Paragon Plus Environment
Page 7 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
in the membrane system are then parallel to the z-axis due to the large distance, which is implemented by the z independent delta-functions in eq. (9). Thus, this particular choice has a physical justification. The contributions to the normal and in-plane pressure are obtained by averaging Tijint over an area A in the x-y-plane. For the normal pressure we find: ∞
∞
int
PN
Z Z E ext X dy δ(x − xα )δ(y − yα ) dx qα θ(z − zα ) =− A α −∞
ext
=−
E 2A
X
(10)
−∞
qα sign(z − zα ),
(11)
α
under assumption of charge neutrality
P
α qα
= 0 (which must be fulfilled when using particle
mesh Ewald for the electrostatic interactions). For the in-plane pressure we simply get
PLint = 0
(12)
due to the Kroenecker deltas in the interaction tensor. We therefore obtain the following addition to the lateral pressure profile due to the system’s interactions with the external electric field:
π add (z) =
E ext X qα sign(z − zα ). 2A α
(13)
The importance of this contribution is illustrated in the Supporting Information. In the analysis the z-axis is divided into slabs (see Simulation Details). The midpoints of these slabs are then used as the discretized z-values in eq. (13).
7
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 8 of 22
Flexoelectricity For lipid membranes Petrov 4 coins two flexoelectric effects: the direct and the converse effect. The direct effect refers to curvature-induced polarization: Upon deforming/curving the membrane, the distribution of charge/dipoles rearranges causing polarization of the headgroup regions, thus creating a potential difference across the membrane. The converse effect refers to the opposite situation: a membrane potential, or equivalently a transmembrane electric field, causes the bilayer to spontaneously curve, or at least a bending moment to be induced, because of the electrostatic interactions of the field with the charges and dipoles in the head groups. In this paper we focus on the converse effect: we measure the bending moment induced in a bilayer as a function of the applied membrane potential. Petrov 4 defines a relation between the induced spontaneous curvature and the applied electric field. Since the size of the electric field varies a lot between the interior of the bilayer, the head group region, and the bulk water region due to the different dielectric constants of these regions, we propose that a better way of characterizing the effect is via the membrane potential difference, ∆Φ. We therefore propose the following relation between ∆Φ and the resulting bending moment (similar to the definition used by Bivas and Kozlov 17 , who investigated the effects of direct flexoelectricity on membrane elasticity):
β = λf ∆Φ,
(14)
where λf is a coupling constant which depends on the composition of the bilayer.
Simulation Details The coarse-grained (CG) MARTINI force field 18,19 is used to model the lipid bilayer (required files were obtained from the MARTINI web site 20 in October 2014). A CG force field is sufficient to demonstrate the existence of, and to describe the procedure for the estimation 8
ACS Paragon Plus Environment
Page 9 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
of membrane potential induced bending moment. Furthermore, the current analysis protocol (see below) makes atomistic resolution practically impossible. The simulated system consists of a symmetric bilayer with 256 palmitoyl-oleoyl-phosphatidylcholine (POPC) molecules in each leaflet, solvated with 15.800 MARTINI polarizable water particles 21 . Na+ and CL – ions are present to ensure a salt concentration of 150 mM. This system resembles the experimental ones in both composition and ionic concentration. All simulations are performed using GROMACS 4.6 22,23 . The system was equilibrated using the Berendsen coupling 24 for both temperature and pressure. Production runs are carried out in the NPT ensemble for 125 ns using a time step of 20 fs. The temperature is controlled using the Nosé-Hoover thermostat 25 , with a reference value of 323 K, and the pressure is coupled semi-isotropically using the Parrinello-Rahman barostat 26 , with a reference value of 1 bar. The van der Waal’s forces are shifted to zero from 0.9 to 1.2 nm. Electrostatic interactions are calculated using the particle mesh Ewald (PME) method 27,28 with a real space cutoff of 1.4 nm and a global relative dielectric constant of ǫr = 2.5 as recommended by Yesylevskyy et al. 21 , when using polarizable water. To generate membrane potential differences, external electric fields of various strengths are applied normal to the plane of the bilayer (along the z-axis). The strengths are: 0 (no field), 0.0335, 0.067, 0.101 and 0.134 V nm−1 . These values are chosen to generate transmembrane potential differences of about 0, 0.5, 1.0, 1.5 and 2 V. The actual values of the potential difference however depend on the box dimensions which are slightly affected by the interactions between the electric field and the particles (see Results). The largest potentials used are one order of magnitude larger than standard biological values, to maximize the response of the system. The results will still be valid for smaller membrane potentials. 40 copies of each simulation were run with different starting velocity distributions. By monitoring the equilibration of the charge distribution we found that the final 75 ns of each simulation should be used for analysis. The pressure profiles are obtained using a custom version of GROMACS 4.5 which calcu-
9
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 10 of 22
lates the stress tensor locally as described by Vanegas et al. 15 . Their version is based on the one by Ollila et al. 14 . Here all forces are re-calculated and distributed locally in the system resulting in a 3D representation of the stress tensor. Both implementations use the stress tensor definition by Irving and Kirkwood 16 . For the calculations of the local stress tensors all parameters were kept the same as in the original simulations except for the handling of electrostatic interactions. As explained by Vanegas et al. 15 this implementation does not calculate the contributions from the reciprocal space when calculating electrostatic interactions, originally calculated using PME. Instead a plain cutoff is used, with typical values of 2.0 to 2.5 nm. However, in our systems, the applied electric field causes a more organized distribution of ions which need to interact across the lipid bilayer which is about 4 nm thick. Such long-range electrostatic forces will not be captured by the usual cut-off values, and we have used a cut-off of 6.0 nm. Each system dimensions must be at least 12.0 nm, and it although such sizes can be simulated on the all-atom level, the large number of simulations, and the very significant increase in the CPU-hours spent on the analysis render the all-atom simulations practically impossible to implement and analyse. The z-axis of the simulation box is divided into 75 slabs, yielding thicknesses around 0.2 nm, consistent with previous studies. 14
Results and Discussion To verify that the system reacts to the increasing applied electric field we investigate the resulting membrane potential, the distribution of ions, the dipole moment of the lipid head groups, and the lateral pressure profile. To determine the magnitude of the membrane potentials, we use Poisson’s equation: 29 ρ ∇2 Φ = − , ǫ
10
ACS Paragon Plus Environment
(15)
Page 11 of 22
where Φ(r) is the electric potential and ρ(r) is the charge density of the system. From this the electric potential can be found by integrating the charge density twice. The integrations are performed with respect to z since both Φ and ρ are independent of x and y. Doing this calculation for the charge distribution for each of the five systems averaged over all 40 simulations, yields the membrane potentials shown in Figure 1. 0 −0.5 Φ [V]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
−1
0.000 V nm−1 0.034 V nm−1 0.067 V nm−1
−1.5
0.101 V nm−1 0.134 V nm−1
−2
−5
0
5
z [nm] Figure 1: Membrane potentials, generated by applying transmembrane electric fields, for the five systems calculated using Poisson’s equation and the average charge distribution for each system. The membrane potential differences for the five systems are: 0.00, 0.48, 0.98, 1.44 and 1.76 V. These numbers are below the ones aimed for. This is because the interactions between the lipids and the applied electric field result in an increase of the area of the bilayer, thus resulting in a decrease of the height of the simulation box. Finding the average height of the box for each system and multiplying by the applied electric field strength gives potential differences within 5 % of the values found using Poisson’s equation. The MARTINI model of POPC used here has no other charges besides a P-N dipole in the head group, which points towards the water region. This simple model makes predicting the behavior of the head groups relatively simple when an electric field is applied. Dipoles in the upper leaflet are aligned with the electric field while dipoles in the lower are aligned oppositely. This causes additional stess in the lower leaflet because of the resulting torque on the dipoles. We can observe this change by measuring the normal component of the dipole moment, µz , separately for each leaflet, as a function of the applied membrane potential 11
ACS Paragon Plus Environment
The Journal of Physical Chemistry
(Figure 2). In the upper leaflet µz increases while it decreases in the lower, verifying that the systems react as expected with increasing applied electric fields. 9 µz /area [D/nm2 ]
−3 8 −4 7 −5 6 0 0.5 1 1.5 ∆Φ [V]
0 0.5 1 1.5 ∆Φ [V]
(a) Upper leaflet
−6
(b) Lower leaflet
Figure 2: Normal component of the head group dipole moment as a function of the membrane potential difference, ∆Φ, for the (a) upper leaflet, and (b) lower leaflet. Error bars represent standard deviations of 40 simulations for each value of ∆Φ. 10 ρ(z) [C/cm3 ]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Page 12 of 22
Na+ 0 0.00 V
Cl
0.48 V
−
0.98 V 1.44 V
−10 −6
1.74 V
−4
−2
0
2
4
6
z [nm] Figure 3: Charge densities for the Na+ and Cl – ions as a function of z for the five different membrane potentials. A symmetric ion distribution is observed for the system with no applied potential (Figure 3). As the electric field strength is increased, Na+ ions will gather at the lower leaflet while Cl – ions will gather at the upper. If the two types of ions interact differently with the two leaflets then this asymmetry could also contribute to the generation of a bending moment together with the dipole effects, 12
ACS Paragon Plus Environment
Page 13 of 22
described above. Apart from that possibility the contribution of the ions is mainly to screen the water regions from the electric field in the membrane, such that almost all of the potential drop will happen in the membrane. If the ionic strength is so weak that this is not the case, then the bending moment measured at fixed membrane potential will be reduced, since the electric field in the membrane will be weaker, see for instance Loubet et al. 30 200 100
π(z) [bar]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
0 −100 −200
0.00 0.48 0.98 1.44 1.76
−300 −400
−4
−2
0
2
V V V V V
4
z [nm] Figure 4: Average pressure profile for each value of the membrane potential difference showing increasing asymmetry with increasing potential difference. There is an increased asymmetry in the pressure profile as the membrane potential increases (Figure 4). The change is most evident in the upper leaflet, where the (negative) interfacial pressure increases while the (positive) repulsive head group term decreases. The altered pressure profile should influence mechanosensitive proteins, which react to the force they experience from their lipid surroundings. If the membrane potential alters the pressure profile this can be a mechanism used by proteins sensitive to this potential. Investigation of this possibility is however beyond the scope of the present work. Finally, to test whether an increase in the bending moment with the potential is measurable, as eq. (14) predicts, we use eq. (4) to calculate β for all simulations (Figure 5).
13
ACS Paragon Plus Environment
The Journal of Physical Chemistry
Mean values Fit
2 β [pJ/nm]
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
0 −2 −4 −6
0
0.5 1 ∆Φ [V]
1.5
Figure 5: The bending moment, β, as a function of the membrane potential difference, ∆Φ, together with a fit of eq. (14) constrained through the origin. Error bars represent uncertainties of the 40 simulations for each value of ∆Φ. Despite the large spread of the measurements that arises from very high pressure fluctuations, it is clear that the size of the average bending moment increases with increasing potential difference. For this system we expect a positive bending moment which means λf < 0 which is also the case. This fits expectations: overall the pressure should increase in the lower leaflet, because of the additional stress from the dipoles. In the upper leaflet, the pressure will decrease because of increased ordering of dipoles, and corresponding lower head group repulsion. Thus, with a lower pressure in the upper leaflet, and a higher pressure in the lower leaflet, a negative first moment is expected. The fit indicates that eq. (14) holds. Note that the fit is constrained through the origin because the spontaneous curvature of the reference system should be zero. The fluctuations can possibly be reduced by simulating larger systems, but stress tensor calculations become very expensive, because we use a cutoff of 6 nm for these calculations (see Methods section). Bilayers in the MARTINI force field tolerate much higher electrical potentials, 31 and thus our calculation of λf should be viewed as an estimate. Please note that even for all-atom force fields, the actual potentials required for eliciting response like electroporation do not match
14
ACS Paragon Plus Environment
Page 14 of 22
Page 15 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
experimental values. In experiments, pore formation occurs at 0.25 to 0.5 V. 32 In all-atom simulations, the minimum value is nearly 1.4 V 33 or as large as 3 V. 32 Thus, one cannot expect an accurate quantitative estimate of lambda from all-atom simulations either. Figure 5 can hence be used to predict the order of magnitude of the curvature induced by a given value of ∆Φ. From the fit |λf | is found to be 1.75 ± 0.38 pJ/Vm, which is comparable with values obtained experimentally for, e.g. phosphatidylcholine or phosphatidylserine bilayers, under similar ionic conditions. 4 The bending rigidity for a MARTINI POPC bilayer system, with the same salt concentration, has been calculated by Loubet et al. 34 to, roughly, κ = 10 × 10−20 J. This value does not depend on the membrane potential. 34 Inserting these values together with a transmembrane potential difference of ∆Φ = 100 mV one finds a spontaneous curvature between 1.37 × 106 and 2.13 × 106 m−1 corresponding to radii between 0.47 and 0.73 µm, which are of the same order as biologically relevant sizes.
Conclusions In conclusion, we use coarse grained MD simulations to show that the head group dipole moments, the lateral pressure profile and the spontaneous curvature of a POPC bilayer all change systematically when a transmembrane electrical potential difference is applied. Our results predict that typical membrane potentials can induce spontaneous curvatures of biologically relevant sizes and thus that the flexoelectric effect contributes in inducing curvature in biological membranes. The ramifications of the phenomena are easy to appreciate keeping in light the highly curved membranes of cells and the sensitivity of several key proteins and membrane processes to the transmembrane electric potential. More realistic estimates of the correlation between electric potential and spontaneous curvature will emerge from all-atom MD simulations, which take into account all the dipoles in the lipids, such as those present in the acyl tails and the glycerol head groups.
15
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
Acknowledgement D.S.B. and H.K. are funded by a Lundbeck Foundation Young Invesigator Grant. M.A.L. acknowledges funding from Danish Council for Independent Research — Natural Sciences (FNU), grant number 4002-00428B. The authors acknowledge DeIC National HPC Centre, University of Southern Denmark, for computing resources.
Supporting Information Available Data showing the equilibration of the system (area per lipid and charge density as a functiom of time). Normal pressure plots illustrating the importance of the additional contribution to the local stress tensor (eq. (13)).
References (1) McMahon, H. T.; Boucrot, E. Membrane Curvature at a Glance. J. Cell Sci. 2015, 128, 1065–1070. (2) Lodish, H.; Berk, A.; Matsudaira, P.; Kaiser, C. A.; Kriger, M.; Scott, M. P.; Zipursky, S. L.; Darnell, J. Molecular Cell Biology, 5th ed.; Freeman: New York, NY, U.S.A., 2004. (3) Ahmadpoor, F.; Sharma, P. Flexoelectricity in Two-Dimensional Crystalline and Biological Membranes. Nanoscale 2015, 7, 16555–16570. (4) Petrov, A. G. Flexoelectricity of Model and Living Membranes. Biochim. Biophys. Acta 2001, 1561, 1–25. (5) Morris, C. E. Mechanosensitive Ion Channels. J. Membrane Biol. 1990, 113, 93–107.
16
ACS Paragon Plus Environment
Page 16 of 22
Page 17 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
(6) Árnadóttir, J.; Chalfie, M. Eukaryotic Mechanosensitive Channels. Annu. Rev. Biophys. 2010, 39, 111–137. (7) Haswell, E.; Phillips, R.; Rees, D. Mechanosensitive Channels: What Can They Do and How Do They Do It? Structure 2011, 19, 1356–1369. (8) García-Añoveros, J.; Corey, D. P. Mechanosensation: Touch at the Molecular Level. Curr. Biol. 1996, 6, 541–543. (9) Ollila, O. H. S.; Louhivuori, M.; Marrink, S.; Vattulainen, I. Protein Shape Change Has a Major Effect on the Gating Energy of a Mechanosensitive Channel. Biophys. J. 2011, 100, 1651–1659. (10) Helfrich, W. Elastic Properties of Lipid Bilayers: Theory and Possible Experiments. Z. Naturforsch. C 1973, 28, 693–703. (11) Safran, S. A. Statistical Thermodynamics of Surfaces, Interfaces, and Membranes; Addison-Wesley: Reading, MA, USA, 1994. (12) Pezeshkian, W.; Chaban, V. V.; Johannes, L.; Shillcock, J.; Ipsen, J. H.; Khandelia, H. The Effects of Globotriaosylceramide Tail Saturation Level on Bilayer Phases. Soft Matter 2015, 11, 1352–1361. (13) Różycki, B.; Lipowsky, R. Spontaneous Curvature of Bilayer Membranes from Molecular Simulations: Asymmetric Lipid Densities and Asymmetric Adsorption. J. Chem. Phys. 2015, 142, DOI: 10.1063/1.4906149. (14) Ollila, O. H. S.; Risselada, H. J.; Louhivuori, M.; Lindahl, E.; Vattulainen, I.; Marrink, S. J. 3D Pressure Field in Lipid Membranes and Membrane-Protein Complexes. Phys. Rev. Lett. 2009, 102, DOI: 10.1103/PhysRevLett.102.078101. (15) Vanegas, J. M.; Torres-Sánchez, A.; Arroyo, M. Importance of Force Decomposition
17
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
for Local Stress Calculations in Biomembrane Molecular Simulations. J. Chem. Theory Comput. 2014, 10, 691–702. (16) Irving, J. H.; Kirkwood, J. G. The Statistical Mechanical Theory of Transport Processes. IV. The Equations of Hydrodynamics. J. Chem. Phys. 1950, 18, 817–829. (17) Bivas, I.; Kozlov, M. M. Effect of Flexoelectricity on the Curvature Elasticity of the Membrane. Liq. Cryst. 1990, 8, 813–817. (18) Marrink, S. J.; de Vries, A. H.; Mark, A. E. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 108, 750–760. (19) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812–7824. (20) MARTINI Website (accessed October 1, 2014): Downloads. http://www.cgmartini. nl/index.php/downloads. (21) Yesylevskyy, S. O.; Schäfer, L. V.; Sengupta, D.; Marrink, S. J. Polarizable Water Model for the Coarse-Grained MARTINI Force Field. PLoS Comput. Biol. 2010, 6, DOI: 10.1371/journal.pcbi.1000810. (22) Berendsen, H.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43– 56. (23) Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D. et al. GROMACS 4.5: a HighThroughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845–854.
18
ACS Paragon Plus Environment
Page 18 of 22
Page 19 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
(24) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684–3690. (25) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255–268. (26) Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190. (27) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N log(N ) Method for Ewald Sums in Large Systems. J. Phys. Chem. 1993, 98, 10089–10092. (28) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A Smooth Particle Mesh Ewald Method. J. Chem. Phys 1995, 103, 8577–8593. (29) Griffiths, D. J. Introdution to Electrodynamics, 4th ed.; Pearson: Upper Saddle River, NJ, U.S.A., 1999. (30) Loubet, B.; Hansen, P. L.; Lomholt, M. A. Electromechanics of a Membrane with Spatially Distributed Fixed Charges: Flexoelectricity and Elastic Parameters. Phys. Rev. E 2013, 88, DOI: 10.1103/PhysRevE.88.062715. (31) Polak, A.; Velikonja, A.; Kramar, P.; Tarek, M.; Miklavcic, D. Electroporation Threshold of POPC Lipid Bilayers with Incorporated Polyoxyethylene Glycol (C12 E8 ). J. Phys. Chem. B 2015, 119, 192–200. (32) Tieleman, P. D. The Molecular Basis of Electroporation. BMC Biochem 2004, 5, 10. (33) Böckmann, R. A.; de Groot, B. L.; Kakorin, S.; Neumann, E.; Grubmüller, H. Kinetics, Statistics, and Energetics of Lipid Membrane Electroporation Studied by Molecular Dynamics Simulations. Biophys. J. 2008, 95, 1837–1850.
19
ACS Paragon Plus Environment
The Journal of Physical Chemistry
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
(34) Loubet, B.; Lomholt, M. A.; Khandelia, H. Tension Moderation and Fluctuation Spectrum in Simulated Lipid Membranes Under an Applied Electric Potential. J. Chem. Phys 2013, 139, DOI: 10.1063/1.4826462.
20
ACS Paragon Plus Environment
Page 20 of 22
Page 21 of 22
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60
The Journal of Physical Chemistry
Graphical TOC Entry Apply electric eld
E
21
ACS Paragon Plus Environment
Apply electric eld
The Journal of Physical Page Chemistry 22 of 22
1 2ACS Paragon Plus Environment E 3 4 5 +