Subscriber access provided by Maastricht University Library
Article
Exploring the Local Elastic Properties of Bilayer Membranes using Molecular Dynamics Simulations Gilles Paul Pieffet, Alonso Botero, Günther H. Peters, Manu Forero-Shelton, and Chad Leidy J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/jp504427a • Publication Date (Web): 17 Oct 2014 Downloaded from http://pubs.acs.org on October 25, 2014
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 29
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
Exploring the Local Elastic Properties of Bilayer Membranes using Molecular Dynamics Simulations Gilles Pieffet†,§, Alonso Botero†, Günther H. Peters‡, Manu Forero-Shelton† and Chad Leidy*,†
† Department of Physics, Universidad de los Andes, Carrera 1 No 18A 10, Bogotá, Colombia
‡ Department of Chemistry, Technical University of Denmark, 2800 Kgs. Lyngby, Denmark
§ Present address: Facultad Ciencias, Universidad Antonio Nariño. Bogotá, Colombia
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
Page 2 of 29
ABSTRACT: Membrane mechanical elastic properties regulate a variety of cellular processes involving local membrane deformation, such as ion channel function and vesicle fusion. In this work, we use molecular dynamics simulations to estimate the local elastic properties of a membrane. For this, we calculated the energy needed to extract a DOPE lipid molecule, modified with a linker chain, from a POPC bilayer membrane using the umbrella sampling technique. Although the extraction energy entails several contributions related not only to elastic deformation but also to solvation, careful analysis of the potential of mean force (PMF) allowed us to dissect the elastic contribution. With this information, we calculated an effective linear spring constant for the DOPC membrane of 44±4 kJ·nm-2·mol-1, in agreement with experimental estimates. The membrane deformation profile was determined independently during the stretching process at molecular detail, allowing us to fit this profile to a previously proposed continuum elastic model. Through this approach, we calculated an effective membrane spring constant of 42 kJ·nm-2·mol-1, which is in good agreement with the PMF calculation. Furthermore, the solvation energy we derived from the data is shown to match with the solvation energy estimated from critical micelle formation constants. This methodology can be used to determine how changes in lipid composition or the presence of membrane modifiers can affect the elastic properties of a membrane at a local level.
Keywords: lipid, deformation, solvation, continuum elastic model INTRODUCTION Membrane elastic properties are known to regulate the activity of integral membrane proteins by altering the energy cost of membrane deformation associated with protein conformational changes1–5. Ion channel dimerization and vesicle fusion depend strongly on these elastic properties1–3, which are influenced by changes in membrane composition through the introduction of lipid metabolites such as 2
ACS Paragon Plus Environment
Page 3 of 29
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
lysophospholipids6, cholesterol7 and polyunsaturated fatty acids8. For example, the presence of cholesterol has been shown to significantly increase the bending rigidity of the membrane for particular lipid species7. A key parameter is the effective spring constant, which is determined by the energy required to stretch the membrane in a direction normal to the surface. The effective spring constant has been used to characterize the energy contribution due to the insertion of a protein inclusion presenting a hydrophobic mismatch with the membrane9,10, for example, in studies of gramicidin dimerization3 and its modulation11. These elastic properties are normally measured experimentally using micropipette aspiration12–17, fluctuation analysis18,19 or can be estimated through x-ray diffraction20–22. However, these experimental approaches only provide average global elastic behavior, and cannot be used to investigate elastic behavior at a local level. Here, we estimate the effective spring constant at a local level for a membrane using molecular dynamics (MD) simulations and umbrella sampling23 (US) simulations. We chose a 1-palmitoyl-2-oleoylphosphatidylcholine (POPC) membrane where we included a 1,2-dioleoyl-phosphatidyl ethanolamine lipid molecule modified with a linker chain (mDOPE) as the inclusion being pulled. This is meant to ensure consistency with future force spectroscopy pulling experiments using an avidin-biotin coupled system A potential of mean force (PMF) calculation of a lipid extraction process has been previously reported on a different membrane system24. However, the authors focused on determining the energetics of desorption, and the elastic properties of the membrane were not considered. From our simulations, we were also able to estimate the energy cost of pulling the mDOPE lipid molecule from a POPC bilayer membrane and separate the energy contributions from membrane deformation and solvation. Alternatively, we calculated the deformation energy based on the continuum elastic model2,4 using the average deformation of the membrane at different distances of 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 29
pulled mDOPE, and compare the resulting effective membrane spring constant with that one calculated from the PMF profile. We expect that our methodology will provide a useful tool for estimating the energy contribution to the mechanical properties of integral membrane proteins resulting from local membrane deformation. Furthermore, the method could be used to study the effect of lipid composition on these elastic properties, and how the addition of membrane modifiers can influence these properties.
4
ACS Paragon Plus Environment
Page 5 of 29
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
METHODS Molecular dynamics and umbrella sampling simulations The simulated system consists of a lipid bilayer comprising 127 POPC lipids and one modified DOPE lipid (mDOPE). The modified DOPE corresponds to a DOPE molecule to which a small extension was added (Fig. 1A and 1C) that models the small linker connecting a DOPE lipid to a Biotin tag (Fig. 1B).
A
B
C
D Figure 1. Chemical structures of the lipids used in this study: (A) POPC, (B) biotinylated DOPE and (C) modified DOPE (mDOPE) designed with a linker. (D) Snapshots taken from simulations corresponding to three different stages of the extraction process.
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 29
The force field describing the lipid molecules (including the modified DOPE lipid and its extension) uses the Berger parameters25, while water molecules were represented by the SPC model26. Molecular Dynamics (MD) simulations were performed using the Gromacs software package version 4.5.427–31. Temperature was kept constant by weak coupling the system32 to a temperature bath at 300 K, with water and lipid coupled separately. Pressure coupling was performed using the Parrinello-Rahman barostat33 with a semi-isotropic scheme where lateral and perpendicular pressures are coupled independently34,35. The pressure was set to 1 bar. Short range electrostatic interactions below 1 nm were calculated at every time step while long range electrostatic interactions were calculated using the particle mesh Ewald method36. Simulations used a 1 nm cutoff for Lennard-Jones interactions. The bond lengths and the angles of the water molecules were constrained using the SETTLE algorithm37, while in the other molecules, bond lengths were constrained using the LINCS algorithm38. A time step of 2 fs was used in all simulations. The equilibrium MD simulations (equilibration of 200 ps and production run of 50 ns), performed to assess the membrane behavior, included 2460 water molecules. The umbrella sampling23 (US) simulations, performed to calculate the PMF, were done using a larger simulation box containing 7552 water molecules. The PMF corresponds to the free energy since it is the energy spent along the reaction coordinate. We switched to the US technique, in order to obtain statistics in regions of the conformational space with low probabilities but that are still important to the extraction process (such as when the lipid is partially outside of the membrane). The US method applies a biasing potential to restrict the sampling around these regions. This umbrella potential was applied to the atom C56 (Fig. 1C) as a harmonic constraint. The extraction process of mDOPE from inside the membrane to its complete exposure in the water phase was split into 54 umbrella windows between which the lipid was shifted by 0.058 nm on average. The starting structures of the umbrella windows were obtained from a steered MD simulation where mDOPE was pulled to increasing z-positions (i.e. along the bilayer normal) and into the solvent. The unbiased potential of mean force was then calculated using the weighted 6
ACS Paragon Plus Environment
Page 7 of 29
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
histogram analysis method (WHAM)39 as implemented in the Gromacs software package utility g_wham40. Unbiased mean forces in each umbrella window could also be directly calculated from the biased forces using the derivation of the potential of mean force from the umbrella integration method41,42. From these equations, when the reaction coordinate is equal to the (biased) mean value of the reaction coordinate in a given umbrella simulation, the unbiased mean force is equal to the biased mean force for that umbrella window. Calculation of the number of solvent exposed carbon atoms We find that by plotting the PMF as a function of the number of carbon atoms that become exposed to the solvent, it becomes easy to distinguish the purely elastic region of the membrane from a mixed region that involves deformation and lipid extraction. Through this methodology, it also becomes easy to compare our results with results derived from critical micelle concentration data. At a molecular level, it is not necessarily straightforward to define what it means for an atom to be solvent exposed when it is located at or near the solvent interface. For this calculation, we considered atoms to be exposed to the solvent when their z-positions were higher than the average z-position of the phosphate groups of the neighboring lipids located within a 1.2 nm radius from the center of mDOPE. With this procedure we define the number of solvated carbon atoms as the average number of CH2 alkyl groups per chain that are being exposed to the solvent. Continuum elastic model The continuum elastic model proposed by Nielsen et al.2,4 was applied to calculate independently the membrane force constant using the structural information of the lipids from the MD simulations. This model has originally been developed to predict the membrane deformation induced by the presence of a membrane protein with a hydrophobic length different from the hydrophobic length of the membrane 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 29
(hydrophobic mismatch). The membrane will tend to deform around the protein in order to reduce the water exposure of the hydrophobic areas of the protein. This concept can also be applied to our system, since the process of lipid extraction is analogous; the different windows of the umbrella sampling simulations can be seen as a series of protein inclusions of increasing hydrophobic lengths h0 that result in membrane deformation and increasing hydrophobic length mismatches u0 = u(r0) as shown in Fig. 2 (with h0 = u0 + d0/2). Here, the radius of the inclusion, r0, is equivalent to the radius of mDOPE, which varied between 4.4 and 4.5 Å in the US simulations.
Figure 2: Lipid pulling induced bilayer deformation. The membrane in the umbrella sampling windows has a radially dependent deformation of the monolayer containing the constrained lipid, where the deformation is denoted z = u(r) at the distance r from the center of the constrained lipid. The maximum monolayer thickness h0 is observed at the boundary position r0 of the constrained lipid, when the maximum membrane deformation u0 occurs. The boundary conditions of the continuum elastic model are also represented in the figure.
8
ACS Paragon Plus Environment
Page 9 of 29
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
We briefly review our adaptation of the continuum model of Nielsen et al.2 used to estimate the energetics of a deformed lipid bilayer. Assuming rotational symmetry, deformations are measured relative to the asymptotic membrane height by the function u(r), where r is the radial (horizontal) distance to the center point located between the two chains of the pulled lipid at the height of the phosphate positions of the first neighboring lipids, as shown in Fig. 2. In Nielsen et al.2, the change in the free energy due to a deformation of a lipid bilayer, assuming the same deformation for both layers, is modeled by the surface integral
∆Ga =
Ka 2 1 2 2 2 da u + K ( ∇ u ) + α ∇ u c 2 ∫ d 02
[1]
where the three terms in brackets account for the three contributions to the deformation energy, as defined in Nielsen and Andersen4: i) compression-expansion with deformation modulus Ka; ii) splaydistortion, with splay-distortion modulus Kc; and iii) surface tension, with interfacial tension α. The constant d0 is equal to the asymptotic separation between the two layers. For a single layer, the free energy is simply half that of the bilayer, with d0/2 interpreted as the thickness of the single layer. In equilibrium, the deformation satisfies the Euler-Lagrange equations derived from Eq.1, which may be cast as
∇ 4 u − γ∇ 2u + β u = 0 ,
[2]
where
γ=
α Kc
,β =
Ka , d02 K c
[3]
The boundary conditions are chosen such that at infinity, 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
lim u(r) = 0, lim u′(r) = 0 r→∞
r→∞
Page 10 of 29
[4]
and at r0,
u(r0 ) = u0 , u′(r0 ) = u0
[5]
as shown in Fig. 2. As explained in Nielsen et al.2, a general solution to the Euler-Lagrange equation is of the form
u(r) = ∑ Ai K 0 (ki r) [6] i
where K0(x) is the zeroth-order modified Bessel function of the first kind, and the ki are the four roots of the characteristic equation
k 4 − γ k 2 + β = 0. [7] In general, and specifically for the parameter values appropriate to our problem, the four roots of the characteristic equation are complex; from symmetry under parity and under complex conjugation of the *
characteristic equation, they can be written as: k+, k+*, -k+, and -k+ where Re(k+) >0, and Im(k+)>0. Of these, only the roots k+ and k+* have positive real parts, corresponding to solutions that decay at infinity. Thus, by fitting the conditions at r0, the solution to the deformation profile can be written as
u(r) = 2 Re [ A+ K 0 (k+r)]
[8]
where
A+ =
k− K1 (k− r0 )u0 + K 0 (k−r0 )s k− K1 (k− r0 )K 0 (k+r0 ) − k+ K 0 (k− r0 )K1 (k+ r0 )
[9]
10
ACS Paragon Plus Environment
Page 11 of 29
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
Substituting into Eq. 1 yields an expression for the deformation contribution to the free energy that depends on only the three geometric parameters s, r0, and u0:
∆Ga = −2π r0 K c Re A+ k+2 (sK 0 (k+r0 ) + k+u0 K1 (k+r0 )) − πα u0r0 s
[10]
To compare the results of the layer containing the pulled lipid in the US simulations and in the continuum model, we used the reference values given in Table 2 of Nielsen et al.2 for the constants Ka = 142.5 pN/nm, Kc = 28.5 pN·nm, and α = 3.0 pN/nm. We also used the fixed value d0 = 38.4 Å, which was determined by averaging the distance between the phosphate groups of the bilayer membrane in equilibrium. To match properly the boundary conditions of the continuum model, the value of r0 for each data set was taken to be the horizontal distance between the center of the pulled lipid (as previously defined) and its closest lipid neighbor in the layer, while u0 was taken to be the average height of the first neighbors of the pulled lipid instead of the pulled lipid. These criteria were chosen to match our system as close as possible to the protein inclusion model proposed by Nielsen et al., avoiding artifacts in the choice of u0 resulting from the decoupling of the pulled lipid from the membrane. The only parameter in the elastic model that cannot be reliably inferred from the simulated profiles due to the noise is the slope s at r0. Therefore for each value h0 of the reaction path the slope s was taken as the value that maximized the likelihood between the observed profile and the continuum elastic profile u(r, s) from Eq. 8 throughout the whole range of r. With the dependence of the slope s on h0 determined, the effective dependence of the energy on h0 can be obtained using Eq. 10.
11
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 12 of 29
RESULTS AND DISCUSSION After equilibration of the membrane, one of the lipids is pulled from the membrane until it is extracted. Using the umbrella sampling technique, the force that the lipid exerts on a harmonic potential is then obtained for a collection of points along the extraction trajectory chosen in order to obtain statistics in regions of the conformational space with low probabilities but that are still important to the extraction process. In this manner, it is possible to obtain the free energy profile or potential of mean force43, from which we estimate the effective spring constant of the membrane (see Methods section for more details).
Three regimes can be discerned during the pulling A modified lipid (mDOPE) was pulled from an equilibrated POPC membrane via a linker chain as described in the Methods section. Umbrella sampling simulations were then performed at different distances from the membrane to calculate the forces necessary to keep the lipid in place at different positions during pulling. The force was calculated using a harmonic potential that restrains the pulled atom, where we define Zpull as the distance between the C56 carbon (Fig. 1B) of mDOPE and the center of mass of the membrane (x = y = z = 0 Å). In Fig. 3A, the mean force is plotted as a function of the displacement distance (Zpull) for mDOPE. The linker adds about 0.5 nm to the average height of the phosphate group of the unmodified DOPE lipid, resulting in a distance of ~2.5 nm for the phosphate group of mDOPE with respect to the center of the membrane when in equilibrium at the minimum force value. The force-displacement curve, resulting from the umbrella sampling simulations shown in Fig. 3A, shows three distinct regimes. Initially, the force increases linearly up to a displacement distance of 4 nm. 12
ACS Paragon Plus Environment
Page 13 of 29
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
The linear force buildup in this regime suggests that this corresponds to the elastic deformation of the membrane. The force then reaches a plateau at 4.0 nm where it remains constant before dropping down to zero at 4.7 nm. The plateau is identified as the second regime. When the force returns to zero, which we identify as the third regime, the mDOPE snaps out of the membrane as shown in snapshot 3 of Fig. 1D. To further characterize these three distinct regimes, we plotted the average z-position of the lipids neighboring mDOPE (annular lipids) as a function of Zpull. Fig. 3B shows that the neighboring lipid position is coupled to Zpull up to 4.0 nm, after which the neighboring lipids remain at an approximately constant z-position (regime 2). At 4.7 nm, the POPC lipids surrounding mDOPE relax to their initial position (regime 3), which correlates to the force going back to zero. It should be noted that only one leaflet is shown in Fig. 3B because the other one is not deformed during the process. Implications for the continuum elastic model are discussed in the corresponding section. Fig. 3C shows the PMF profile of the extraction process. We see that the energy needed to completely extract mDOPE is ~106 kJ·mol-1, which is similar to the value calculated by Tieleman and Marrink24 at 80 kJ·mol-1 for a DPPC membrane, also calculated via the umbrella sampling method. When looking at the PMF profile, we find the same three regimes previously observed in Fig. 3A: i) a presumably quadratic relationship for the energy associated with the deformation of the membrane, ii) a linear relationship associated to the extraction regime, and finally iii) a constant value associated with the membrane relaxation regime after mDOPE is extracted from the membrane.
13
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 14 of 29
Figure 3. (A) Mean force applied to mDOPE to keep the lipid at a given distance along the reaction coordinate, Zpull. (B) Average height of the annular lipid phosphates within a 1.25 nm radius of mDOPE
14
ACS Paragon Plus Environment
Page 15 of 29
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
as a function of Zpull. (C) Potential of mean force of mDOPE as a function of Zpull when extracted from a POPC membrane. The solvation of the lipid marks the transition between the first and the second regime Fig. 4A and B show the solvent exposed surface area of the acyl chains (solid line) calculated using the double cubic lattice method44, and the number of solvated carbon atoms (see methods for details) as a function of Zpull. We can see that in the first regime, the exposed surface area of the acyl chain remains constant and low. The non-zero solvent exposure in this regime corresponds to interfacial water molecules that contact the top groups of the acyl chains. At 3.5 nm, the solvent exposure begins to increase with Zpull. This distance in Fig. 4A (3.5 nm) is smaller than the one observed for the mechanical decoupling of the lipid to its neighbors in Fig. 3A (4.0 nm) suggesting that water exposure takes place before mechanical decoupling. Finally, at 4.7 nm, there is a considerable increase in the exposed solvent area corresponding to the full exposure to the water phase of the acyl chain of mDOPE. This matches the point in Zpull where the force drops to zero. The dashed line in Fig. 4A shows the total surface area of the acyl chains. At 4.7 nm, it merges with the solvent exposed surface area curve confirming that the lipid has been extracted at this point. The downward trend of the surface area as a function of Zpull is concurrent with an increase in the order parameter of the chains (data not shown) as the restrained chains will have less surface area when they bundle. This can be interpreted by pointing out that the acyl chains of mDOPE are initially moving freely in a hydrophobic environment. These chains come together as they are pulled out. When mDOPE is completely outside the membrane, the solvent exposed area of the acyl chain drops in order to minimize the exposure of the hydrophobic chains to the solvent. A
B 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
Page 16 of 29
Figure 4. (A) Solvent accessible surface area and total surface area of the hydrophobic tail of mDOPE as a function of the reaction coordinate, Zpull. This area starts to increase after 3.5 nm (red arrow). (B) Number of CH2 alkyl groups of the lipid tail that are exposed to the solvent as a function of the reaction coordinate, Zpull. The values reported correspond to the average calculated from both chains.
The calculated free energy during lipid solvation is consistent with experimentally determined solution to micelle transfer free energy Fig. 5 shows the PMF profile from Fig. 3C as a function of the number of solvent exposed CH2 acyl groups of the lipid tail (labeled NbC, light blue circles). We also observe the three regimes previously described: i) from 0 to about ~30 kJ/mol, there are no exposed carbons, consistent with membrane deformation; ii) from ~30 to ~100 kJ/mol, there is an increase in energy with carbon number that should correspond to lipid solvation; and finally iii) a plateaus region at ~100 kJ/mol, where the lipid is fully solvated. In order to see if the calculated solvation energies are consistent with experimentally derived energies45, we calculated the transfer free energy of lipids or acyl chains from solution to micelle using the critical micelle concentration (CMC) with the following formula24: 16
ACS Paragon Plus Environment
Page 17 of 29
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 ∆G 0 = −RT ln = RT ln(CMC) . CMC The dark blue circles in Fig. 5 represent the transfer free energy for diacyl chains of varying lengths from 5 to 16 carbon atoms from experimentally determined CMC data45 and a linear fit with respect to the length of the acyl chain. The fit results in a free energy change of ∆Gotransf = 4.4 ± 0.2 kJ·mol-1 per extra CH2 alkyl group added to the acyl chain. Likewise, from the simulation data in Fig. 5, we determined that the change in free energy for carbons 5 to 13 is ∆Gosolv = 5.1 ± 0.5 kJ·mol-1 for each extra CH2 alkyl group solvated. There is good agreement between the solvation energy per alkyl group deduced from the CMC data and the simulation results for carbons 5 to 13 (i.e. Zpull values between 4.0 and 4.7 nm). However, it is clear from the figure that the slope of the curve between 30 and 60 kJ/mol, corresponding to the solvation of the first five carbon atoms, is different from that for atoms 5 to 13 (see mixed regime in Table 1). This corresponds to a Zpull value between 3.5 and 4.0 nm. Moreover, in Fig. 3A, it appears that the elastic regime ends at about 4 nm while in Fig. 4, solvation starts at 3.5nm. Thus, the region between 3.5 and 4.0 nm can then be interpreted as a mixed deformation/extraction behavior, corresponding to the transition between the regimes. Consequently, the change in energy seen in Fig. 5 at NbC = 0 corresponds to a pure deformation regime where no solvent exposure is measured. Also, the change in energy for NbC between 5 and 13 corresponds to a pure solvation regime, as it matches well with the CMC45 data and no further deformation is seen in Fig. 3B (regime 2).
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
Page 18 of 29
Figure 5. Experimental solvation free energy and extraction free energy. The light blue circles represent the calculated free energy for the extraction of mDOPE as a function of the number of CH2 alkyl groups (NbC) that are becoming exposed to the solvent. The dashed red line corresponds to the fit, with a correlation of r = 0.96. The dark blue circles represent the transfer free energy change for a series of diacyl PC chains of various lengths and calculated from experimentally determined CMCs45. The dashed blue line corresponds to the best fit for the transfer free energy points, with a correlation r > 0.99. Most of the energy spent in extracting the lipid is due to its solvation energy. The energy contributions to the pulling of the lipid of each of the three regimes — deformation, mixed deformation/solvation and solvation — are given in Table 1 and correspond to energies of 28 kJ·mol-1, 36 kJ·mol-1 and 42 kJ·mol-1, respectively. By extrapolating the solvation data, the expected contribution of the solvation energy to the mixed regime (the first 5 carbons) should be 26 ± 3 kJ·mol-1, so the contribution of the membrane deformation would be the remaining energy, 10 ± 3 kJ·mol-1. It follows that the total membrane deformation energy is 38 ± 3 kJ·mol-1, while the total solvation energy is about 68 ± 3 kJ·mol-1.
18
ACS Paragon Plus Environment
Page 19 of 29
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
Table 1: Energy associated to the deformation and extraction regime of the membrane.
Regime
Deformation
Mixed
Solvation
Relaxation
Zpull (nm)
2.50 - 3.50
3.50 - 4.00
4.00 - 4.74
> 4.74
NbC
0
1-5
5 - 13
> 13
Free energy (kJ·mol-1)
28
36
42
0
Extrapolated energy contributions (Pure deformation and solvation)
38 ± 3
-
68 ± 3
-
Membrane deformation properties from umbrella sampling agree with those derived from a continuum elastic model A linear regression of the points between 2.5 nm and 3.5 nm in the force-distance curve obtained from umbrella sampling (Fig. 3A) results in an apparent membrane spring constant of 44 ± 4 kJ·nm-2·mol1
. Similarly, fitting the regime of the PMF corresponding to the membrane deformation in Fig. 3C to a
quadratic expression also leads to an apparent membrane spring constant of 40 kJ·nm-2·mol-1. These results are discussed below together with the boundary conditions of the model as the experimental results depend on them. An alternate method for calculating the elastic spring constant is to fit the deformation profile of the membrane during the extraction process by using the continuum elastic model proposed by Nielsen et al.2, and originally developed to quantify the hydrophobic mismatch of proteins. Here, mDOPE at increasing Zpull position can be seen as a series of protein inclusions of increasing hydrophobic lengths 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
Page 20 of 29
that result in membrane deformation and increasing hydrophobic length mismatches (see Methods section). In Fig. 6A and 6B, we plot the average z-position of the phosphates of all the POPC lipids from each leaflet as a function of their horizontal distance to the phosphate of mDOPE. In Fig. 6A, mDOPE is in its equilibrium position (Zpull = 2.49 nm) in the absence of an applied force, and no deformation is observed. A noticeable deformation of the neighboring POPC lipids can be perceived in Fig. 6B when mDOPE is close to being extracted out of the membrane (Zpull = 4.74 nm). At this point, the membrane has already reached a limit in its level of deformation (since Zpull > 4.0 nm, Fig. 2A), and the energetic cost of further deforming the membrane is now higher than the energy cost of acyl chain exposure to the water. The deformation energy can be calculated by fitting the positions of the neighboring POPC lipids using the continuum elastic model proposed by Nielsen et al.2,4. Fitting the deformation energy scatter-plot to a quadratic expression yields an apparent spring constant for the membrane. The deformation energy and the fit are shown in Fig. 6D. The apparent membrane spring constant calculated with the continuum elastic model using both experimental data and the structural data extracted from the simulation is ~84 kJ·nm-2·mol-1, which is when both leaflets are assumed to be undergoing a symmetric deformation. The deformation energy for a single leaflet should therefore be ~42 kJ·nm-2·mol-1. This last value is comparable to the value of 44 kJ·nm-2·mol-1 found with the umbrella sampling method, which only involves the deformation of one leaflet as Fig. 6B shows that the lower leaflet is not deformed during the pulling. The agreement between the two methods indicates that the US technique is a quick and adequate method for calculating the effective spring constant of the membrane by MD simulations.
20
ACS Paragon Plus Environment
Page 21 of 29
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
A
B
C
D
Figure 6: Top: Average height of the phosphate groups as a function of their distance from the phosphate of the pulled mDOPE lipid. (A) At the beginning of the pulling at Zpull = 2.49 nm and (B) just before the lipid leaves the membrane at Zpull = 4.74 nm. Bottom: Continuum elastic model. (C) Slope s and (D) energy associated with the membrane deformation as a function of the reaction coordinate h0 and obtained by applying the continuum elastic model to the US simulations of the deformation regime. It is noteworthy to mention that the continuum elastic model depends on experimentally reported stretching and bending moduli to estimate the effective spring constant (see Methods section), while the US technique relies only on the force field.
21
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 22 of 29
Additionally, the original derivation of the continuum elastic model requires that the slope s between the inclusion and the first neighbors (Fig. 2) follow one of two different criteria: i) a no-slip criterion between neighboring lipids, denominated the constrained boundary condition, where s is equal to zero or ii) a minimum energy criterion, denominated the free boundary condition, where s is chosen to minimize the deformation energy. In our adaptation of the model, however, we actually fit this parameter to the simulation data. Fig. 6C shows a series of fits of the parameter s at various h0 values corresponding to different Zpull positions of mDOPE. We find that s varies linearly with the height h0 of the hypothetical inclusion, which is consistent with the free boundary condition, and indicates considerable slippage between mDOPE and its first neighbors. In our case, the free boundary condition (i.e. when s is assumed to be non zero) is expected to prevail since lipid packing defects at the interface between the membrane and the lipid being pulled would be reduced by the hydrocarbons filling the voids created at the interface, leading to slippage. The values of the apparent spring constant calculated from the US simulations (44 ± 4 kJ·nm-2·mol-1) and the continuum elastic model (~42 kJ·nm-2·mol-1) are comparable to the experimental estimate3 calculated in the case of free boundary conditions (between 23 ± 6 and 31 ± 2 kJ·nm-2·mol-1). The higher value in the simulations is consistent with the fact that they were performed with a solvent more polar than the experiments and therefore, the hydrophobic coupling is expected to be stronger.
CONCLUSION We pulled an mDOPE lipid from a POPC bilayer, and we calculated the free energy profile of the lipid extraction from the US simulations. The energy cost for completely extracting the lipid was found to be
22
ACS Paragon Plus Environment
Page 23 of 29
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
~106 kJ·mol-1. We were able to distinguish three different regimes in the process, as well as a mixed regime during the transition between the first two. The first regime was characterized by a nonlinear increase in energy with respect to the pulling coordinate, in a manner consistent with membrane deformation. From the PMF, we were able to calculate the apparent spring constant of the membrane, which was found to be 44 kJ·nm-2·mol-1. We also calculated the spring constant using the positions of the lipids surrounding mDOPE by means of the continuum elastic model, and we found that the two were similar and consistent with previous experimental estimates3, with the latter having a value of ~42 kJ·nm-2·mol-1. Moreover, the free parameter s describing the slope between the inclusion and first neighbors was found to be nonzero, indicating slippage. The second regime appeared after a short transition and was characterized by a linear increase in energy with respect to the pulling coordinate and was consistent with the solvation of the lipid´s tail. When comparing the extraction free energy with experimentally derived solvation free energy from critical micelle concentration studies45, we found that the energy required to solvate an extra CH2 group were similar giving 4.4 ± 0.2 kJ·mol-1 for the experimental data and 5.1 ± 0.5 kJ·mol-1 for the simulation data. The third regime, which was characterized by a constant energy and the relaxation of the membrane, is consistent with the snapping of the lipid out of the membrane and its solvation. The methodology presented provides a computational tool for determining membrane elastic properties at a local level (within 2-3 lipid shells) and is applicable for studying local changes in membrane elastic properties as a function of composition and in the presence of membrane modifiers.
23
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 24 of 29
* Corresponding author. E-mail:
[email protected] Supporting Information Available: Supplemental Methods, additional figure, equations and references. This material is available free of charge via the Internet at http://pubs.acs.org. ACKNOWLEDGMENTS Simulations were performed at the Danish Center for Scientific Computing at the University of Southern Denmark. Dr. G.Pieffet was supported by a postdoctoral position financed by the Department of Physics at the Universidad de los Andes. REFERENCES (1)
Huang, H. W. Deformation Free Energy of Bilayer Membrane and Its Effect on Gramicidin Channel Lifetime. Biophys. J. 1986, 50, 1061–1070.
(2)
Nielsen, C.; Goulian, M.; Andersen, O. S. Energetics of Inclusion-Induced Bilayer Deformations. Biophys. J. 1998, 74, 1966–1983.
(3)
Lundbæk, J. A.; Andersen, O. S. Spring Constants for Channel-Induced Lipid Bilayer Deformations Estimates Using Gramicidin Channels. Biophys. J. 1999, 76, 889–895.
(4)
Nielsen, C.; Andersen, O. S. Inclusion-Induced Bilayer Deformations: Effects of Monolayer Equilibrium Curvature. Biophys. J. 2000, 79, 2583–2604.
(5)
Wiggins, P.; Phillips, R. Analytic Models for Mechanotransduction: Gating a Mechanosensitive Channel. Proc. Natl. Acad. Sci. 2004, 101, 4071–4076.
(6)
Lundbæk, J. A.; Andersen, O. S. Lysophospholipids Modulate Channel Function by Altering the Mechanical Properties of Lipid Bilayers. J. Gen. Physiol. 1994, 104, 645–673.
(7)
Lundbæk, J. A.; Birn, P.; Girshman, J.; Hansen, A. J.; Andersen, O. S. Membrane Stiffness and Channel Function†. Biochemistry (Mosc.) 1996, 35, 3825–3830.
(8)
Bruno, M. J.; Koeppe, R. E.; Andersen, O. S. Docosahexaenoic Acid Alters Bilayer Elastic Properties. Proc. Natl. Acad. Sci. 2007, 104, 9638–9643. 24
ACS Paragon Plus Environment
Page 25 of 29
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
(9)
Parton, D. L.; Klingelhoefer, J. W.; Sansom, M. S. P. Aggregation of Model Membrane Proteins, Modulated by Hydrophobic Mismatch, Membrane Curvature, and Protein Class. Biophys. J. 2011, 101, 691–699.
(10)
Jensen, M. Ø.; Mouritsen, O. G. Lipids Do Influence Protein Function—the Hydrophobic Matching Hypothesis Revisited. Biochim. Biophys. Acta BBA - Biomembr. 2004, 1666, 205–226.
(11)
Kelkar, D. A.; Chattopadhyay, A. Modulation of Gramicidin Channel Conformation and Organization by Hydrophobic Mismatch in Saturated Phosphatidylcholine Bilayers. Biochim. Biophys. Acta BBA - Biomembr. 2007, 1768, 1103–1113.
(12)
Evans, E. A. Minimum Energy Analysis of Membrane Deformation Applied to Pipet Aspiration and Surface Adhesion of Red Blood Cells. Biophys. J. 1980, 30, 265–284.
(13)
Evans, E. A. Bending Elastic Modulus of Red Blood Cell Membrane Derived from Buckling Instability in Micropipet Aspiration Tests. Biophys. J. 1983, 43, 27–30.
(14)
Evans, E.; Rawicz, W. Entropy-Driven Tension and Bending Elasticity in Condensed-Fluid Membranes. Phys. Rev. Lett. 1990, 64, 2094–2097.
(15)
Henriksen, J. R.; Ipsen, J. H. Measurement of Membrane Elasticity by Micro-Pipette Aspiration. Eur. Phys. J. E Soft Matter 2004, 14, 149–167.
(16)
Ly, H. V.; Longo, M. L. The Influence of Short-Chain Alcohols on Interfacial Tension, Mechanical Properties, Area/Molecule, and Permeability of Fluid Lipid Bilayers. Biophys. J. 2004, 87, 1013– 1033.
(17)
Longo, M. L.; Ly, H. V. Micropipet Aspiration for Measuring Elastic Properties of Lipid Bilayers. Methods Mol. Biol. Clifton NJ 2007, 400, 421–437.
(18)
Faucon, J. F.; Mitov, M. D.; Méléard, P.; Bivas, I.; Bothorel, P. Bending Elasticity and Thermal Fluctuations of Lipid Membranes. Theoretical and Experimental Requirements. J. Phys. 1989, 50, 2389–2414.
(19)
Henriksen, J.; Rowat, A. C.; Ipsen, J. H. Vesicle Fluctuation Analysis of the Effects of Sterols on Membrane Bending Rigidity. Eur. Biophys. J. 2004, 33, 732–741.
(20)
Nagle, J. F.; Tristram-Nagle, S. Structure of Lipid Bilayers. Biochim. Biophys. Acta 2000, 1469, 159– 195.
(21)
Tristram-Nagle, S.; Nagle, J. F. Lipid Bilayers: Thermodynamics, Structure, Fluctuations, and Interactions. Chem. Phys. Lipids 2004, 127, 3–14.
25
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 26 of 29
(22)
Raghunathan, M.; Zubovski, Y.; Venable, R. M.; Pastor, R. W.; Nagle, J. F.; Tristram-Nagle, S. Structure and Elasticity of Lipid Membranes with Genistein and Daidzein Bioflavinoids Using XRay Scattering and MD Simulations. J. Phys. Chem. B 2012, 116, 3918–3927.
(23)
Torrie, G. M.; Valleau, J. P. Monte Carlo Free Energy Estimates Using Non-Boltzmann Sampling: Application to the Sub-Critical Lennard-Jones Fluid. Chem. Phys. Lett. 1974, 28, 578–581.
(24)
Tieleman, D. P.; Marrink, S.-J. Lipids Out of Equilibrium: Energetics of Desorption and Pore Mediated Flip-Flop. J. Am. Chem. Soc. 2006, 128, 12462–12467.
(25)
Berger, O.; Edholm, O.; Jahnig, F. Molecular Dynamics Simulations of a Fluid Bilayer of Dipalmitoylphosphatidylcholine at Full Hydration, Constant Pressure, and Constant Temperature. Biophys. J. 1997, 72, 2002–2013.
(26)
Berendsen, H. J. C.; Postma, J. P. M.; Van Gunsteren, W. F.; Hermans, J. Interaction Models for Water in Relation to Protein Hydration. Intermolecular Forces 1981, 11, 331–342.
(27)
Berendsen, H. J. C.; van der Spoel, D.; van Drunen, R. GROMACS: A Message-Passing Parallel Molecular Dynamics Implementation. Comput. Phys. Commun. 1995, 91, 43–56.
(28)
Lindahl, E.; Hess, B.; van der Spoel, D. GROMACS 3.0: A Package for Molecular Simulation and Trajectory Analysis. J Mol Mod 2001, 7, 306–317.
(29)
Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, Flexible, and Free. J. Comput. Chem. 2005, 26, 1701–1718.
(30)
Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447.
(31)
Pronk, S.; Páll, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; Spoel, D. van der; et al. GROMACS 4.5: A High-Throughput and Highly Parallel Open Source Molecular Simulation Toolkit. Bioinformatics 2013, 29, 845–854.
(32)
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.
(33)
Parrinello, M.; Rahman, A. Polymorphic Transitions in Single Crystals: A New Molecular Dynamics Method. J. Appl. Phys. 1981, 52, 7182–7190.
(34)
Zhang, Y.; Feller, S. E.; Brooks, B. R.; Pastor, R. W. Computer Simulation of Liquid/liquid Interfaces. I. Theory and Application to Octane/water. J. Chem. Phys. 1995, 103, 10252–10266.
(35)
Anézo, C.; de Vries, A. H.; Höltje, H.-D.; Tieleman, D. P.; Marrink, S.-J. Methodological Issues in Lipid Bilayer Simulations. J. Phys. Chem. B 2003, 107, 9424–9433.
26
ACS Paragon Plus Environment
Page 27 of 29
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
(36)
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.
(37)
Miyamoto, S.; Kollman, P. A. Settle: An Analytical Version of the SHAKE and RATTLE Algorithm for Rigid Water Models. J. Comput. Chem. 1992, 13, 952–962.
(38)
Hess, B.; Bekker, H.; Berendsen, H. J.; Fraaije, J. G. LINCS: A Linear Constraint Solver for Molecular Simulations. J. Comput. Chem. 1997, 18, 1463–1472.
(39)
Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. THE Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011–1021.
(40)
Hub, J. S.; de Groot, B. L.; van der Spoel, D. g_wham—A Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation Estimates. J. Chem. Theory Comput. 2010, 6, 3713–3720.
(41)
Kästner, J.; Thiel, W. Bridging the Gap between Thermodynamic Integration and Umbrella Sampling Provides a Novel Analysis Method: “Umbrella Integration.” J. Chem. Phys. 2005, 123, 144104.
(42)
Kästner, J. Umbrella Sampling. Wiley Interdiscip. Rev. Comput. Mol. Sci. 2011, 1, 932–942.
(43)
Kirkwood, J. G. Statistical Mechanics of Fluid Mixtures. J. Chem. Phys. 1935, 3, 300–313.
(44)
Eisenhaber, F.; Lijnzaad, P.; Argos, P.; Sander, C.; Scharf, M. The Double Cubic Lattice Method: Efficient Approaches to Numerical Integration of Surface Area and Volume and to Dot Surface Contouring of Molecular Assemblies. J. Comput. Chem. 1995, 16, 273–284.
(45)
Avanti Polar Lipids. https://www.avantilipids.com/index.php?option=com_content&view=article&id=1703&Itemid=4 22 (accessed Oct 14, 2014).
27
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 28 of 29
28
ACS Paragon Plus Environment
Page 29 of 29
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
TOC graphic:
29
ACS Paragon Plus Environment