Materials Science and Engineering C 36 (2014) 49–56
Contents lists available at ScienceDirect
Materials Science and Engineering C journal homepage: www.elsevier.com/locate/msec
Molecular dynamics simulation of interlayer water embedded in phospholipid bilayer Won Bae Han a, Suk Jun Kim b, Hyeun Hwan An a, Hee-Soo Kim a, Yongdeok Kim a, Chong Seung Yoon a,⁎ a b
Division of Materials Science and Engineering, Hanyang University, Seoul 133-791, Republic of Korea Samsung Electronics, Yongin-si, Gyeonggi-do 446-712, Republic of Korea
a r t i c l e
i n f o
Article history: Received 2 August 2013 Received in revised form 14 October 2013 Accepted 22 November 2013 Available online 4 December 2013 Keywords: Lipid bilayer Interlayer water Molecular dynamics simulation
a b s t r a c t 1,2-dioleoyl-sn-glycero-3-phosphocholine lipid bilayer with a thin layer of water molecules inserted in the hydrophobic region was simulated at 300 K to observe the pore structure formation during escape of the water molecules from the hydrophobic region. The transformation of the water slab into a cylindrical droplet in the hydrophobic region, which locally deformed the lipid monolayer, was prerequisite to the pore formation. If the thickness of the interlayer water was increased beyond a critical value, the local deformation was suppressed as such deformation would rupture the lipid bilayer. Hence, it was demonstrated that the pore structure formation or local permeability of the lipid membrane is closely related to the rigidity of the lipid membrane. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Lipid bilayer, which is formed spontaneously due to the amphipathic tendency of the lipid molecules, constitutes an integral component of the cell membrane. The lipid bilayer separates the cell interior from the extracellular environment while the semi-fluidic nature of the lipid bilayer also allows the selected ions and organic molecules to be transported across the plasma membrane via transmembrane protein [1]. Penetration of chemical species across a lipid bilayer without the aid of transmembrane protein is an important issue, especially in direct delivery of exogenous molecules into cells [2,3]. Electric field (electroporation) [4,5], ultrasound (sonoporation) [6,2], and laser irradiation [7] have been introduced to cause localized formation of pores within the lipid barrier which temporarily increases the local permeability of the cell membrane in the targeted area. In sonoporation in which ultrasonic wave is used to temporarily increase the membrane permeability, the ultrasound and ultrasound contrast agent in combination create bubble cavitation [8,9]. The bubble formation near the membrane is expected to generate shock waves that cause localized disorder within the lipid bilayer. The localized disorder, in turn, can temporarily increase the permeability of the lipid bilayer [10]. It is also conjectured that the following implosion of the cavitation bubbles may form localized jet streams acting as micro needles which penetrate the cell membrane [11]. However, the actual pore formation or temporary increase in the local permeability in the cell membrane is difficult to access experimentally due to the short duration, molecular dynamic (MD) simulation has been mainly employed to study the instantaneous pore formation [12,13]. K. Koshiyama et al. ⁎ Corresponding author. E-mail address: [email protected] (C.S. Yoon). 0928-4931/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.msec.2013.11.033
have suggested that the shock waves generated by repetitive ultrasound irradiation could lead to penetration of water molecules into the hydrophobic region of the lipid bilayer [14,13]. In a recent article, the same group also demonstrated that the water molecules embedded into the hydrophobic region can trigger a short-lived (nanosecond scale) formation of a water-filled pore structure in a model bilayer as the water molecules escape from the hydrophobic region [15]. The simulation showed that the pore formation required disruption and severe deformation of the lipid bilayer due to the segregation of the water molecules embedded among the hydrophobic tails. Sonoporation in which ultrasonic wave is used to temporarily increase the membrane permeability, the ultrasound shock wave can carry the water molecules inside of a model lipid bilayer on the picosecond time scale through MD simulation [15]. As the presence of water molecules in the hydrophobic region is energetically unstable, the water molecules will spontaneously condense into a droplet or into a continuous slab and eventually escape from the lipid bilayer. The escape of the water molecules will likely depend on the rigidity and the geometry of the water coalescence. In this work, increasingly thicker film of water molecules in a predetermined geometry is inserted into the interlayer space between the hydrophobic tails of a model lipid bilayer to explicitly observe the segregation of the water molecules and subsequent deformation of the lipid bilayer during the escape of the water molecules. The stiffness of the DOPC bilayer was also altered to see its effect on the water escape from the DOPC bilayer. For the simulation of the model bilayer, 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) was used instead of the previously studied 1,2-dipalmitoyl-sn-glycero-3phosphocholine (DPPC) because it was experimentally observed that the permeability of the phospholipid membrane can be greatly affected by its phase state (i.e. gel phase vs. liquid crystalline) [16]. DPPC with saturated aliphatic tails transforms from a gel phase to a liquid
50
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
Fig. 1. Simulated lipid bilayer model at 300 K consisting of 128 DOPC molecules (64 on top and bottom sides) with the water layer inserted between the two hydrophobic tails of the bilayer.
Fig. 3. Snap shots from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer having the 0.7-nm-thick interlayer water slab in the membrane (a), (b) after 2 ns; (c), (d) after 6 ns; and (e), (f) after 8 ns ((b),(d),and (f) show the water molecules only).
Fig. 2. Snap shots from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer having the 0.5-nm-thick interlayer water slab in the membrane (a), (b) after 2 ns; (c), (d) after 6 ns; and (e), (f) after 8 ns ((b),(d), and (f) show the water molecules only).
Fig. 4. (a), (b) Snap shots from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer having the 1.0-nm-thick interlayer water slab in the membrane after 10 ns ((b) shows the water molecules only).
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
crystalline state at a comparatively high temperature (~40 °C in excess water) [17], thus making the DPPC bilayer relatively rigid at room temperature; hence, MD simulation for the DPPC bilayer is usually performed above room temperature which imparted extra kinetic energy to the water molecules. DOPC, on the other hand, having cis double bond in the acyl chain transforms from the rigid gel phase to the liquid crystalline phase at −20 °C [17] so that the DOPC bilayer is in a liquidcrystalline state at room temperature and remains flexible. Hence, using DOPC enabled the simulation to be performed at 300 K while providing sufficient flexibility to the model bilayer. 2. Simulation methods Atomistic MD simulation of hydrated DOPC bilayer was performed using GROMACS 4.5.3 simulation package [18] and GROMOS96 53a6 force field modified for DOPC [19]. All initial configurations were built from an equilibrated hydrated DOPC lipid bilayer consisting of 128 DOPC molecules (64 in each lipid monolayer) and 3600 water molecules [19]. To insert the water molecules in the hydrophobic region,
51
the top and bottom DOPC monolayer was separated by 0.5, 0.7, and 1.0 nm. The resulting interfacial gap was filled with water molecules. The simulation box size was 6.5104 nm in the x- and y-directions and 9.2937 nm, 9.493 nm, and 9.7937 in the z-direction for the 0.5-, 0.7-, and 1.0-nm-thick water interlayers, respectively. After the water interlayer insertion, the system is relaxed by minimizing the total energy using the steepest descent method followed by a short 10 ps equilibration in the NVT ensemble at 300 K. The NPT simulation was performed with a time step of 2 fs at 300 K by weakly coupling lipids and water separately to a heat bath at 300 K [20] with a coupling time constant of 0.1 ps. Pressure was maintained at 1 bar with the semi-isotropic scaling using the weakly coupled Berendsen barostat [20] and a coupling constant of 1.0 ps. The simple point charge model was used for water. Within the cutoff distance of 1.2 nm, the electrostatic interaction was treated as pairwise summation over all neighboring particles using the Coulomb law and outside the cutoff distance the particle-mesh Ewald method was used for computational efficiency [21] with periodic boundary conditions in all three directions. The maximum Fourier spacing was 0.16 nm and the 4th order fitting function was used. A cutoff of
Fig. 5. Spatial (x-axis, y-axis and z-axis) distribution of interlayer water molecules for the cases of (a) the 0.5-nm-thick water slab; (b) the 0.7-nm-thick water slab; and (c) the 1.0-nmthick water slab at different time intervals.
52
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
the interlayer) and 1.0 nm (1380 water molecules in the interlayer) to study the effect of the volumes of the water molecules embedded in the hydrophobic region on the membrane stability. Fig. 2 shows snap shots from the NPT simulation of the DOPC bilayer having the 0.5-nm-thick interlayer water slab in the membrane. The water molecules were immediately rearranged into a single cylinder extending through the lipid bilayer after 2 ns as can be seen from Fig. 2(a). The cylindrical structure can be clearly seen in Fig. 2(b) which visualizes only the water molecules. The diameter of the embedded cylinder formed by the water molecules was approximately 2 nm which agrees with the value if assumed that the entire embedded water molecules were used in forming the cylinder. After 6 ns (see Fig. 2(c) and (d)), a pore structure was formed in the bottom lipid monolayer near the region where the lipid bilayer was most severely distorted due to the inclusion of the water column. Through the pore, the water molecules spontaneously diffused out of the hydrophobic region. Eventually after 8 ns (see Fig. 2(e) and (f)), the majority of the water molecules escaped from the hydrophobic region and the lipid bilayer returned to its original flat configuration. A similar sequence was observed in the NPT
Fig. 5 (continued).
1.2 nm was also used for the van der Waals interactions and bond constraints were handled using the LINCS algorithm with the highest order in the expansion of the constraint coupling matrix and the single correction iteration for rotational lengthening [22]. VMD package was used for both model building and visualization of the simulated DOPC bilayer systems [23]. Production run time was 10 ns for the 0.5-nm-thick water layer and 14 ns for the 0.7- and 1.0 nm-thick water layer cases. Snapshots of the simulation run were extracted from the trajectory file every 1 ns for comparison. 3. Simulation results Fig. 1 illustrates the lipid bilayer model at 300 K consisting of 128 DOPC molecules (64 on top and bottom sides) with the water layer inserted between the two hydrophobic tails of the bilayer. The interlayer water was introduced by raising the top lipid layer from the bottom lipid layer by 0.5 nm and solvating the entire system with water. The 0.5-nm-thick interlayer contained 725 water molecules that were initially distributed uniformly in a thin slab. Thickness of this water interlayer was increased further to 0.7 nm (1004 water molecules in
Fig. 6. DOPC bilayer thickness for (a) the 0.5-nm-thick water slab; (b) the 0.7-nm-thick water slab; and (c) the 1.0-nm-thick water slab at different time intervals.
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
53
coalescence of multiple cylindrical water columns instead of one large column or formation of spherical droplets did not occur in the lipid bilayer. Instead, the flat slab configuration remained unchanged throughout the simulation as indicated by the z-axis water density profile. Fig. 5 suggests that a water layer trapped in the hydrophobic region of the lipid bilayer will escape first by coalescing into a cylindrical droplet which, in turn, locally deforms the lipid layer. The local deformation causing the lipid layer to bend severely created pore(s) and allowed the water molecules to escape through the pore structure. To better illustrate the local deformation of the lipid layer, the bilayer thickness profile along the x-axis was calculated at 2 ns intervals by averaging the difference of the z-positions of the P atoms in the lower and upper lipid layer and is shown in Fig. 6. As expected, the lipid bilayer thickness greatly increased where the water molecules coalesced for the 0.5- and 0.7-nm-thick water interlayer cases whereas in the case of the 1.0-nm thick water interlayer, the bilayer thickness profile remained unchanged. For both the 0.5- and 0.7-nm-thick water interlayer lipid
Fig. 7. Time-averaged deuterium order parameter for carbon atoms in (a) the Sn-1; and (b) Sn-2 oleoyl chains of DOPC molecules in the lipid layer.
simulation of the DOPC bilayer with the water slab thickness increased to 0.7 nm as shown in Fig. 3. The interlayer water molecules spontaneously formed a cylindrical column and the subsequent pore formation allowed the water molecules to diffuse from hydrophobic region. In contrast, as can be seen from Fig. 4, when the water interlayer thickness was increased to 1.0 nm, the water slab remained flat even after 10 ns and no noticeable distortion of the water slab was observed. To better examine the spatial distribution of interlayer water molecules for the above 3 cases (0.5, 0.7, and 1.0 nm), time-evolution of the local density of the interlayer water molecules (calculated by first creating an index file for the interlayer water molecules and using the g_density command in the GROMACS subroutines) along the directions parallel (x- and y-axis) and normal (z-axis) to the lipid bilayer is plotted in Fig. 5. A large peak in the interlayer water density in the x-axis and a relatively flat density profile (notwithstanding local fluctuations) in the y-axis are consistent with the formation of a cylindrical water column as shown in Fig. 2. The z-axis density profile clearly demonstrates the escape of the interlayer water molecules since the density profile became increasingly flat and the central peak approached zero as the simulation proceeded. With the 0.7-nm-thick water slab, the cylindrical column apparently formed in the y-direction with a larger diameter due to the larger amount of the embedded water molecules as shown in Fig. 5(b). The water molecules were slower to escape from the interior of the bilayer compared to the case with the 0.5-nm-thick water interlayer as indicated by the z-axis profiles in Fig. 5(a) and (b). Fig. 5(c) confirms that the flat water interlayer was stable when the water interlayer thickness was 1.0 nm. Both x- and y-axis interlayer water density profiles were flat with relatively low local fluctuations, suggesting that
Fig. 8. Density profiles of the DOPC head groups in the z-direction for the cases of (a) the 0.5-nm-thick water slab; (b) the 0.7-nm-thick water slab; and (c) the 1.0-nm-thick water slab.
54
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
Fig. 9. Normalized average system potential energies for the DOPC bilayers containing 0.5-, 0.7-and 1.0-nm-thick water slab.
Fig. 10. (a) Snap shot after 10 ns from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer containing 0.5-nm-thick water slab using the stiffened dihedral twisting forces for the lipid hydrocarbon chains (see text), (b) same snap shot as (a) without the DOPC molecules showing the water molecules only.
layers, the bilayer thickness increased approximately by 1 nm (0.5-nmthick water layer) and 2 nm (0.7-nm-thick water layer) from the initial respective unrelaxed configurations during the water escape. If each lipid layer is deformed by more than ~1 nm in the z-direction by curving,then pores are generated in the lipid layer. From the data, the critical lipid layer thickness before the pore generation due to the local deformation appears to be 1–2 nm. To quantitatively assess the pore formation in the lipid layer, the carbon-deuterium order parameter, SCD was calculated from SCD = 1/2 b3cos2θ − 1 N where θ is the angle between the acyl chain segment and the membrane normal [24,25]. SCD was calculated for the each carbon atoms in the acyl chains of the lipid layer through which the water molecules diffused out. Fig. 7 plots the average order parameter for each carbon atom in the Sn-1 and Sn-2 oleoyl chains, starting from the ester carbonyl carbon near the lipid head [26,27]. The order parameter which is time-averaged over 10 ns indicates that SCD for the 1.0-nm-thick interlayer case which remained stable is consistently higher than those of the 0.5- and 0.7-nm-thick interlayer cases in which pores were created in the respective lipid layer. Stability of the lipid bilayer can be also inferred from the density profiles of the DOPC head groups in the z-direction calculated by first creating an index file for the head group atoms and using the g_density command in the GROMACS subroutines. As the water molecules diffuse across the bottom DOPC monolayer for the cases of the 0.5- and 0.7-nm-thick interlayer water slabs, the density profile of the head group for the bottom DOPC monolayer, shown in Fig. 8(a) and (b), became progressively broadened due to the pore formation whereas the profiles for the top layer exhibited a sharp peak as the top layer straightened from the escape of the water molecules. As expected, the density profile of the DOPC groups in the case of the 1.0-nm-thick interlayer water slab in Fig. 8(c) was symmetric and maintained the two sharp peaks corresponding to the top and bottom DOPC monolayers. 4. Discussion A water layer trapped in the hydrophobic region of the lipid bilayer will escape first by coalescing into a cylindrical droplet which, in turn, severely deforms the lipid layer. Pore(s) is(are) formed in the deformed lipid layer where the local lipid density is lowest due to the curvature of the deformation. As seen from the simulation, the interlayer water
Table 1 Parameters for the Ryckaert–Bellemans potential used to model the lipid hydrocarbon chains.
Original Berger lipid Modified
C0
C1
C2
C3
C4
C5
9.2789 9.784
12.156 12.156
−13.120 −13.120
−3.097 −3.097
26.240 26.240
−31.495 −32.0
Fig. 11. Spatial (x-axis, y-axis and z-axis) distribution of interlayer water molecules for the DOPC bilayer containing 0.5-nm-thick water slab using the stiffened dihedral twisting forces for the lipid hydrocarbon chains at different time intervals.
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
molecules gradually escape through the pore structure. The spontaneous re-arrangement of the water molecules is thus prerequisite for the pore formation through which the water molecules can diffuse out of the lipid bilayer. This spontaneous transformation of the interlayer water molecules into the cylindrical droplet is energetically favored due to the repulsive force between the hydrophobic tail and the water molecules; however, the re-arrangement of the water molecules inside the bilayer leads to a local deformation of the lipid layer which increases the potential energy of the lipid bilayer. As can be seen from the curves of the system potential energies for the three DOPC bilayers plotted in Fig. 9, the system potential energy curves for the bilayers with 0.5and 0.7-nm-thick interlayer water decay slowly as the deformed DOPC bilayer returns to its initial flat state whereas the potential energy for the bilayers with 1.0-nm-thick interlayer water remains more or less unchanged. The formation of a cylindrical droplet observed during the MD simulation is reasonable since a spherical droplet will further increase the local deformation. As the water layer thickness is increased beyond the critical value, the energy lowered by the coalescence of the water molecules no longer exceeds the energy gained from the local deformation of the lipid bilayer and the water layer remains trapped in the hydrophobic region. In fact, if the water molecules in the 1.0-nm-thick layer were coalesced into a cylindrical bilayer, the diameter of the cylinder would be 2.8 nm which will bring the head group of the lipid monolayer in contact with the hydrophobic tail of the adjacent lipids; hence, such coalescence of water molecules within the hydrophobic region would likely rupture the lipid bilayer. Therefore, it appears that the re-arrangement of the interlayer water molecules and the subsequent escape rate of water molecules will likely depend on the rigidity of the lipid bilayer. To test the proposal that the escape of the interlayer water molecules would depend on the rigidity of the lipid bilayer, the DOPC bilayer was stiffened by changing the force parameters in the Ryckaert– Bellemans potential function which treats the dihedral twisting forces for the lipid hydrocarbon chains [28]. The original and new parameters for the Ryckaert–Bellemans potential are listed in Table 1 and the plot of the modified potential function is provided in the electronic supplemental material. Using the stiffer potential function, the simulation was carried out for 10 ns on the lipid bilayer containing the 0.5-nmthick water interlayer, following the identical procedure as the previous runs. The water molecules in the interlayer quickly coalesced into a cylindrical shape, but were unable to escape from the DOPC bilayer and the water molecules remained trapped within the lipid bilayer as can be seen from the snap of the simulation after 10 ns shown in Fig. 10. The time-evolution of the local density of the interlayer water molecules along the directions parallel (x- and y-axis) and normal (z-axis) to the
Fig. 12. (a) Snap shot after 10 ns from the NPT simulation of the DOPC bilayer containing 0.5-nm-thick water slab with the DOPC molecules kept at 200 K and water molecules at 300 K, (b) same snap shot as (a) without the DOPC molecules showing the water molecules only.
55
lipid bilayer plotted in Fig. 11 confirms that the interlayer indeed remained trapped inside the lipid bilayer. The lipid molecules extended in the vertical direction in order to accommodate the coalescence of the water molecules as the intermolecular potential interaction was left unchanged. However, the increased stiffness of the hydrocarbon chains prevented the formation of a pore in the lipid layer. Furthermore, to see the effect of the slowing down of the lipid layer molecules relative to the water molecules, a simulation run carried out with the DOPC molecules rescaled at 200 K while maintaining the water molecules at 300 K. As expected, the “cold” DOPC molecules considerably retarded the coalescence of the interlayer water molecules as the lipid layer remained rigid. A simulation run lasting for 10 ns showed that the water molecules remained trapped in the DOPC bilayer as can be seen from Fig. 12. Again the interlayer water density shown in Fig. 13 demonstrates the relatively slow coalescence of the water molecules. The simulation results in Figs. 10–13 evince that the rigidity of the lipid bilayer does have a large effect on the movement of the water
Fig. 13. Spatial (x-axis, y-axis and z-axis) distribution of interlayer water molecules for the DOPC bilayer containing 0.5-nm-thick water slab with the DOPC molecules kept at 200 K and water molecules at 300 K.
56
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
molecules embedded in the lipid bilayer and the permeability of the lipid membrane. Although a thick flat layer of water molecules may remain stable microscopically as the MD simulation demonstrated, the finite size of the simulation cell limits the long range fluctuation or undulation [29,30]. The long range fluctuation in the bilayer geometry will likely rupture the bilayer and each lipid monolayer would fold back and form micelles or vesicles. A large simulation cell for an extended period would reproduce the transformation of the lipid layer into other geometrical shapes by inserting the water molecules into the hydrophobic region which is equivalent to separating the two lipid monolayers. 5. Conclusion It was demonstrated through MD simulation that a pore structure in a model lipid bilayer was formed when a layer of water was introduced into the hydrophobic regions. Re-arrangement of the water molecules into a cylindrical droplet which locally deformed the lipid monolayer preceded the pore formation. The pore structure formation or local permeability of the lipid membrane is thus closely related to the rigidity of the lipid membrane. Acknowledgment This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2012000815). Appendix A. Supplementary data The dynamics of the DOPC bilayer having the 0.5, 0.7 and 1.0-nmthick interlayer water slab in the membrane and the plot of the modified Ryckaert–Bellemans potential function can be found online. Supplementary data to this article can be found online at http://dx.doi. org/10.1016/j.msec.2013.11.033.
References [1] C.W. Pratt, K. Cornely, Essential Biochemistry, John Wiley & Sons, New Jersey, 2004. [2] T. Kodama, Y. Tomita, K. Koshiyama, M.J.K. Blomley, Ultrasound Med. Biol. 32 (2006) 905–914. [3] S. Mehier-Humbert, R.H. Guy, Adv. Drug Deliv. Rev. 57 (2005) 733–753. [4] D.P. Tieleman, BMC Biochem. 5 (2004) 10. [5] J. Gehl, Acta Physiol. Scand. 177 (2003) 437–447. [6] S. Bao, B.D. Thrall, D.L. Miller, Ultrasound Med. Biol. 23 (1997) 953–959. [7] W. Tao, J. Wilkinson, E.J. Stanbridge, M.W. Berns, Proc. Natl. Acad. Sci. U. S. A. 84 (1987) 4180–4184. [8] Y. Yamakoshi, T. Miwa, Jpn. J. Appl. Phys. 50 (2011) 07HF01. [9] T. Kodama, N. Tomita, S. Horie, N. Sax, H. Iwasaki, R. Suzuki, K. Maruyama, S. Mori, F. Manabu, J. Electron Microsc. 59 (2010) 187–196. [10] T. Kodama, K. Koshiyama, Y. Tomita, M. Suzuki, T. Yano, S. Fujikawa, Am. Inst. Phys. Conf. Proc. 829 (2006) 34–38. [11] M.J. Haynes, Ultrasound Med. Biol. 26 (2000) 57–62. [12] K. Koshiyama, S. Wada, J. Biomech. 44 (2011) 2053–2058. [13] K. Koshiyama, T. Kodama, T. Yano, S. Fujikawa, Biochim. Biophys. Acta Biomembr. 1778 (2008) 1423–1428. [14] K. Koshiyama, T. Kodama, T. Yano, S. Fujikawa, Biophys. J. 91 (2006) 2198–2205. [15] K. Koshiyama, T. Yano, T. Kodama, Phys. Rev. Lett. 105 (2010) 018105. [16] H.H. An, J.H. Kim, J.H. Lee, D.H. Kwon, H.S. Kim, Y.H. Kim, C.S. Yoon, Chem. Commun. 46 (2010) 9238–9240. [17] A. Terheiden, C. Mayer, K. Moh, B. Stahlmecke, S. Stappert, M. Acet, B. Rellinghaus, Appl. Phys. Lett. 84 (2004) 3891–3893. [18] B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, J. Chem. Theory Comput. 4 (2008) 435–447. [19] D. Poger, A.E. Mark, J. Chem. Theory Comput. 6 (2010) 325–336. [20] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684–3690. [21] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577–8593. [22] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, J. Comput. Chem. 18 (1997) 1463–1472. [23] W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 14 (1996) 33–38. [24] M. Patra, M. Karttunen, M.T. Hyvonen, E. Falck, P. Lindqvist, I. Vattulainen, Biophys. J. 84 (2003) 3636–3645. [25] S. Esteban-Martin, J. Salgado, Biophys. J. 92 (2007) 903–912. [26] A.H. de Vries, A.E. Mark, S.J. Marrink, J. Phys. Chem. B 108 (2004) 2454–2463. [27] S.A. Pandit, S. Vasudevan, S.W. Chiu, R.J. Mashl, E. Jakobsson, H.L. Scott, Biophys. J. 87 (2004) 1092–1100. [28] J.-P. Ryckaert, A. Bellemans, Chem. Phys. Lett. 30 (1975) 123–125. [29] I. Chandrasekhar, D. Bakowies, A. Glättli, P. Hünenberger, C. Pereira, W.F. van Gunsteren, Mol. Simul. 31 (2005) 543–548. [30] N. Ganesan, B.A. Bauer, T.R. Lucas, S. Patel, M. Taufer, J. Comput. Chem. 15 (2011) 2958–2973.
Contents lists available at ScienceDirect
Materials Science and Engineering C journal homepage: www.elsevier.com/locate/msec
Molecular dynamics simulation of interlayer water embedded in phospholipid bilayer Won Bae Han a, Suk Jun Kim b, Hyeun Hwan An a, Hee-Soo Kim a, Yongdeok Kim a, Chong Seung Yoon a,⁎ a b
Division of Materials Science and Engineering, Hanyang University, Seoul 133-791, Republic of Korea Samsung Electronics, Yongin-si, Gyeonggi-do 446-712, Republic of Korea
a r t i c l e
i n f o
Article history: Received 2 August 2013 Received in revised form 14 October 2013 Accepted 22 November 2013 Available online 4 December 2013 Keywords: Lipid bilayer Interlayer water Molecular dynamics simulation
a b s t r a c t 1,2-dioleoyl-sn-glycero-3-phosphocholine lipid bilayer with a thin layer of water molecules inserted in the hydrophobic region was simulated at 300 K to observe the pore structure formation during escape of the water molecules from the hydrophobic region. The transformation of the water slab into a cylindrical droplet in the hydrophobic region, which locally deformed the lipid monolayer, was prerequisite to the pore formation. If the thickness of the interlayer water was increased beyond a critical value, the local deformation was suppressed as such deformation would rupture the lipid bilayer. Hence, it was demonstrated that the pore structure formation or local permeability of the lipid membrane is closely related to the rigidity of the lipid membrane. © 2013 Elsevier B.V. All rights reserved.
1. Introduction Lipid bilayer, which is formed spontaneously due to the amphipathic tendency of the lipid molecules, constitutes an integral component of the cell membrane. The lipid bilayer separates the cell interior from the extracellular environment while the semi-fluidic nature of the lipid bilayer also allows the selected ions and organic molecules to be transported across the plasma membrane via transmembrane protein [1]. Penetration of chemical species across a lipid bilayer without the aid of transmembrane protein is an important issue, especially in direct delivery of exogenous molecules into cells [2,3]. Electric field (electroporation) [4,5], ultrasound (sonoporation) [6,2], and laser irradiation [7] have been introduced to cause localized formation of pores within the lipid barrier which temporarily increases the local permeability of the cell membrane in the targeted area. In sonoporation in which ultrasonic wave is used to temporarily increase the membrane permeability, the ultrasound and ultrasound contrast agent in combination create bubble cavitation [8,9]. The bubble formation near the membrane is expected to generate shock waves that cause localized disorder within the lipid bilayer. The localized disorder, in turn, can temporarily increase the permeability of the lipid bilayer [10]. It is also conjectured that the following implosion of the cavitation bubbles may form localized jet streams acting as micro needles which penetrate the cell membrane [11]. However, the actual pore formation or temporary increase in the local permeability in the cell membrane is difficult to access experimentally due to the short duration, molecular dynamic (MD) simulation has been mainly employed to study the instantaneous pore formation [12,13]. K. Koshiyama et al. ⁎ Corresponding author. E-mail address: [email protected] (C.S. Yoon). 0928-4931/$ – see front matter © 2013 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.msec.2013.11.033
have suggested that the shock waves generated by repetitive ultrasound irradiation could lead to penetration of water molecules into the hydrophobic region of the lipid bilayer [14,13]. In a recent article, the same group also demonstrated that the water molecules embedded into the hydrophobic region can trigger a short-lived (nanosecond scale) formation of a water-filled pore structure in a model bilayer as the water molecules escape from the hydrophobic region [15]. The simulation showed that the pore formation required disruption and severe deformation of the lipid bilayer due to the segregation of the water molecules embedded among the hydrophobic tails. Sonoporation in which ultrasonic wave is used to temporarily increase the membrane permeability, the ultrasound shock wave can carry the water molecules inside of a model lipid bilayer on the picosecond time scale through MD simulation [15]. As the presence of water molecules in the hydrophobic region is energetically unstable, the water molecules will spontaneously condense into a droplet or into a continuous slab and eventually escape from the lipid bilayer. The escape of the water molecules will likely depend on the rigidity and the geometry of the water coalescence. In this work, increasingly thicker film of water molecules in a predetermined geometry is inserted into the interlayer space between the hydrophobic tails of a model lipid bilayer to explicitly observe the segregation of the water molecules and subsequent deformation of the lipid bilayer during the escape of the water molecules. The stiffness of the DOPC bilayer was also altered to see its effect on the water escape from the DOPC bilayer. For the simulation of the model bilayer, 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC) was used instead of the previously studied 1,2-dipalmitoyl-sn-glycero-3phosphocholine (DPPC) because it was experimentally observed that the permeability of the phospholipid membrane can be greatly affected by its phase state (i.e. gel phase vs. liquid crystalline) [16]. DPPC with saturated aliphatic tails transforms from a gel phase to a liquid
50
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
Fig. 1. Simulated lipid bilayer model at 300 K consisting of 128 DOPC molecules (64 on top and bottom sides) with the water layer inserted between the two hydrophobic tails of the bilayer.
Fig. 3. Snap shots from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer having the 0.7-nm-thick interlayer water slab in the membrane (a), (b) after 2 ns; (c), (d) after 6 ns; and (e), (f) after 8 ns ((b),(d),and (f) show the water molecules only).
Fig. 2. Snap shots from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer having the 0.5-nm-thick interlayer water slab in the membrane (a), (b) after 2 ns; (c), (d) after 6 ns; and (e), (f) after 8 ns ((b),(d), and (f) show the water molecules only).
Fig. 4. (a), (b) Snap shots from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer having the 1.0-nm-thick interlayer water slab in the membrane after 10 ns ((b) shows the water molecules only).
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
crystalline state at a comparatively high temperature (~40 °C in excess water) [17], thus making the DPPC bilayer relatively rigid at room temperature; hence, MD simulation for the DPPC bilayer is usually performed above room temperature which imparted extra kinetic energy to the water molecules. DOPC, on the other hand, having cis double bond in the acyl chain transforms from the rigid gel phase to the liquid crystalline phase at −20 °C [17] so that the DOPC bilayer is in a liquidcrystalline state at room temperature and remains flexible. Hence, using DOPC enabled the simulation to be performed at 300 K while providing sufficient flexibility to the model bilayer. 2. Simulation methods Atomistic MD simulation of hydrated DOPC bilayer was performed using GROMACS 4.5.3 simulation package [18] and GROMOS96 53a6 force field modified for DOPC [19]. All initial configurations were built from an equilibrated hydrated DOPC lipid bilayer consisting of 128 DOPC molecules (64 in each lipid monolayer) and 3600 water molecules [19]. To insert the water molecules in the hydrophobic region,
51
the top and bottom DOPC monolayer was separated by 0.5, 0.7, and 1.0 nm. The resulting interfacial gap was filled with water molecules. The simulation box size was 6.5104 nm in the x- and y-directions and 9.2937 nm, 9.493 nm, and 9.7937 in the z-direction for the 0.5-, 0.7-, and 1.0-nm-thick water interlayers, respectively. After the water interlayer insertion, the system is relaxed by minimizing the total energy using the steepest descent method followed by a short 10 ps equilibration in the NVT ensemble at 300 K. The NPT simulation was performed with a time step of 2 fs at 300 K by weakly coupling lipids and water separately to a heat bath at 300 K [20] with a coupling time constant of 0.1 ps. Pressure was maintained at 1 bar with the semi-isotropic scaling using the weakly coupled Berendsen barostat [20] and a coupling constant of 1.0 ps. The simple point charge model was used for water. Within the cutoff distance of 1.2 nm, the electrostatic interaction was treated as pairwise summation over all neighboring particles using the Coulomb law and outside the cutoff distance the particle-mesh Ewald method was used for computational efficiency [21] with periodic boundary conditions in all three directions. The maximum Fourier spacing was 0.16 nm and the 4th order fitting function was used. A cutoff of
Fig. 5. Spatial (x-axis, y-axis and z-axis) distribution of interlayer water molecules for the cases of (a) the 0.5-nm-thick water slab; (b) the 0.7-nm-thick water slab; and (c) the 1.0-nmthick water slab at different time intervals.
52
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
the interlayer) and 1.0 nm (1380 water molecules in the interlayer) to study the effect of the volumes of the water molecules embedded in the hydrophobic region on the membrane stability. Fig. 2 shows snap shots from the NPT simulation of the DOPC bilayer having the 0.5-nm-thick interlayer water slab in the membrane. The water molecules were immediately rearranged into a single cylinder extending through the lipid bilayer after 2 ns as can be seen from Fig. 2(a). The cylindrical structure can be clearly seen in Fig. 2(b) which visualizes only the water molecules. The diameter of the embedded cylinder formed by the water molecules was approximately 2 nm which agrees with the value if assumed that the entire embedded water molecules were used in forming the cylinder. After 6 ns (see Fig. 2(c) and (d)), a pore structure was formed in the bottom lipid monolayer near the region where the lipid bilayer was most severely distorted due to the inclusion of the water column. Through the pore, the water molecules spontaneously diffused out of the hydrophobic region. Eventually after 8 ns (see Fig. 2(e) and (f)), the majority of the water molecules escaped from the hydrophobic region and the lipid bilayer returned to its original flat configuration. A similar sequence was observed in the NPT
Fig. 5 (continued).
1.2 nm was also used for the van der Waals interactions and bond constraints were handled using the LINCS algorithm with the highest order in the expansion of the constraint coupling matrix and the single correction iteration for rotational lengthening [22]. VMD package was used for both model building and visualization of the simulated DOPC bilayer systems [23]. Production run time was 10 ns for the 0.5-nm-thick water layer and 14 ns for the 0.7- and 1.0 nm-thick water layer cases. Snapshots of the simulation run were extracted from the trajectory file every 1 ns for comparison. 3. Simulation results Fig. 1 illustrates the lipid bilayer model at 300 K consisting of 128 DOPC molecules (64 on top and bottom sides) with the water layer inserted between the two hydrophobic tails of the bilayer. The interlayer water was introduced by raising the top lipid layer from the bottom lipid layer by 0.5 nm and solvating the entire system with water. The 0.5-nm-thick interlayer contained 725 water molecules that were initially distributed uniformly in a thin slab. Thickness of this water interlayer was increased further to 0.7 nm (1004 water molecules in
Fig. 6. DOPC bilayer thickness for (a) the 0.5-nm-thick water slab; (b) the 0.7-nm-thick water slab; and (c) the 1.0-nm-thick water slab at different time intervals.
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
53
coalescence of multiple cylindrical water columns instead of one large column or formation of spherical droplets did not occur in the lipid bilayer. Instead, the flat slab configuration remained unchanged throughout the simulation as indicated by the z-axis water density profile. Fig. 5 suggests that a water layer trapped in the hydrophobic region of the lipid bilayer will escape first by coalescing into a cylindrical droplet which, in turn, locally deforms the lipid layer. The local deformation causing the lipid layer to bend severely created pore(s) and allowed the water molecules to escape through the pore structure. To better illustrate the local deformation of the lipid layer, the bilayer thickness profile along the x-axis was calculated at 2 ns intervals by averaging the difference of the z-positions of the P atoms in the lower and upper lipid layer and is shown in Fig. 6. As expected, the lipid bilayer thickness greatly increased where the water molecules coalesced for the 0.5- and 0.7-nm-thick water interlayer cases whereas in the case of the 1.0-nm thick water interlayer, the bilayer thickness profile remained unchanged. For both the 0.5- and 0.7-nm-thick water interlayer lipid
Fig. 7. Time-averaged deuterium order parameter for carbon atoms in (a) the Sn-1; and (b) Sn-2 oleoyl chains of DOPC molecules in the lipid layer.
simulation of the DOPC bilayer with the water slab thickness increased to 0.7 nm as shown in Fig. 3. The interlayer water molecules spontaneously formed a cylindrical column and the subsequent pore formation allowed the water molecules to diffuse from hydrophobic region. In contrast, as can be seen from Fig. 4, when the water interlayer thickness was increased to 1.0 nm, the water slab remained flat even after 10 ns and no noticeable distortion of the water slab was observed. To better examine the spatial distribution of interlayer water molecules for the above 3 cases (0.5, 0.7, and 1.0 nm), time-evolution of the local density of the interlayer water molecules (calculated by first creating an index file for the interlayer water molecules and using the g_density command in the GROMACS subroutines) along the directions parallel (x- and y-axis) and normal (z-axis) to the lipid bilayer is plotted in Fig. 5. A large peak in the interlayer water density in the x-axis and a relatively flat density profile (notwithstanding local fluctuations) in the y-axis are consistent with the formation of a cylindrical water column as shown in Fig. 2. The z-axis density profile clearly demonstrates the escape of the interlayer water molecules since the density profile became increasingly flat and the central peak approached zero as the simulation proceeded. With the 0.7-nm-thick water slab, the cylindrical column apparently formed in the y-direction with a larger diameter due to the larger amount of the embedded water molecules as shown in Fig. 5(b). The water molecules were slower to escape from the interior of the bilayer compared to the case with the 0.5-nm-thick water interlayer as indicated by the z-axis profiles in Fig. 5(a) and (b). Fig. 5(c) confirms that the flat water interlayer was stable when the water interlayer thickness was 1.0 nm. Both x- and y-axis interlayer water density profiles were flat with relatively low local fluctuations, suggesting that
Fig. 8. Density profiles of the DOPC head groups in the z-direction for the cases of (a) the 0.5-nm-thick water slab; (b) the 0.7-nm-thick water slab; and (c) the 1.0-nm-thick water slab.
54
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
Fig. 9. Normalized average system potential energies for the DOPC bilayers containing 0.5-, 0.7-and 1.0-nm-thick water slab.
Fig. 10. (a) Snap shot after 10 ns from the NPT simulation (at 300 K and 1 bar) of the DOPC bilayer containing 0.5-nm-thick water slab using the stiffened dihedral twisting forces for the lipid hydrocarbon chains (see text), (b) same snap shot as (a) without the DOPC molecules showing the water molecules only.
layers, the bilayer thickness increased approximately by 1 nm (0.5-nmthick water layer) and 2 nm (0.7-nm-thick water layer) from the initial respective unrelaxed configurations during the water escape. If each lipid layer is deformed by more than ~1 nm in the z-direction by curving,then pores are generated in the lipid layer. From the data, the critical lipid layer thickness before the pore generation due to the local deformation appears to be 1–2 nm. To quantitatively assess the pore formation in the lipid layer, the carbon-deuterium order parameter, SCD was calculated from SCD = 1/2 b3cos2θ − 1 N where θ is the angle between the acyl chain segment and the membrane normal [24,25]. SCD was calculated for the each carbon atoms in the acyl chains of the lipid layer through which the water molecules diffused out. Fig. 7 plots the average order parameter for each carbon atom in the Sn-1 and Sn-2 oleoyl chains, starting from the ester carbonyl carbon near the lipid head [26,27]. The order parameter which is time-averaged over 10 ns indicates that SCD for the 1.0-nm-thick interlayer case which remained stable is consistently higher than those of the 0.5- and 0.7-nm-thick interlayer cases in which pores were created in the respective lipid layer. Stability of the lipid bilayer can be also inferred from the density profiles of the DOPC head groups in the z-direction calculated by first creating an index file for the head group atoms and using the g_density command in the GROMACS subroutines. As the water molecules diffuse across the bottom DOPC monolayer for the cases of the 0.5- and 0.7-nm-thick interlayer water slabs, the density profile of the head group for the bottom DOPC monolayer, shown in Fig. 8(a) and (b), became progressively broadened due to the pore formation whereas the profiles for the top layer exhibited a sharp peak as the top layer straightened from the escape of the water molecules. As expected, the density profile of the DOPC groups in the case of the 1.0-nm-thick interlayer water slab in Fig. 8(c) was symmetric and maintained the two sharp peaks corresponding to the top and bottom DOPC monolayers. 4. Discussion A water layer trapped in the hydrophobic region of the lipid bilayer will escape first by coalescing into a cylindrical droplet which, in turn, severely deforms the lipid layer. Pore(s) is(are) formed in the deformed lipid layer where the local lipid density is lowest due to the curvature of the deformation. As seen from the simulation, the interlayer water
Table 1 Parameters for the Ryckaert–Bellemans potential used to model the lipid hydrocarbon chains.
Original Berger lipid Modified
C0
C1
C2
C3
C4
C5
9.2789 9.784
12.156 12.156
−13.120 −13.120
−3.097 −3.097
26.240 26.240
−31.495 −32.0
Fig. 11. Spatial (x-axis, y-axis and z-axis) distribution of interlayer water molecules for the DOPC bilayer containing 0.5-nm-thick water slab using the stiffened dihedral twisting forces for the lipid hydrocarbon chains at different time intervals.
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
molecules gradually escape through the pore structure. The spontaneous re-arrangement of the water molecules is thus prerequisite for the pore formation through which the water molecules can diffuse out of the lipid bilayer. This spontaneous transformation of the interlayer water molecules into the cylindrical droplet is energetically favored due to the repulsive force between the hydrophobic tail and the water molecules; however, the re-arrangement of the water molecules inside the bilayer leads to a local deformation of the lipid layer which increases the potential energy of the lipid bilayer. As can be seen from the curves of the system potential energies for the three DOPC bilayers plotted in Fig. 9, the system potential energy curves for the bilayers with 0.5and 0.7-nm-thick interlayer water decay slowly as the deformed DOPC bilayer returns to its initial flat state whereas the potential energy for the bilayers with 1.0-nm-thick interlayer water remains more or less unchanged. The formation of a cylindrical droplet observed during the MD simulation is reasonable since a spherical droplet will further increase the local deformation. As the water layer thickness is increased beyond the critical value, the energy lowered by the coalescence of the water molecules no longer exceeds the energy gained from the local deformation of the lipid bilayer and the water layer remains trapped in the hydrophobic region. In fact, if the water molecules in the 1.0-nm-thick layer were coalesced into a cylindrical bilayer, the diameter of the cylinder would be 2.8 nm which will bring the head group of the lipid monolayer in contact with the hydrophobic tail of the adjacent lipids; hence, such coalescence of water molecules within the hydrophobic region would likely rupture the lipid bilayer. Therefore, it appears that the re-arrangement of the interlayer water molecules and the subsequent escape rate of water molecules will likely depend on the rigidity of the lipid bilayer. To test the proposal that the escape of the interlayer water molecules would depend on the rigidity of the lipid bilayer, the DOPC bilayer was stiffened by changing the force parameters in the Ryckaert– Bellemans potential function which treats the dihedral twisting forces for the lipid hydrocarbon chains [28]. The original and new parameters for the Ryckaert–Bellemans potential are listed in Table 1 and the plot of the modified potential function is provided in the electronic supplemental material. Using the stiffer potential function, the simulation was carried out for 10 ns on the lipid bilayer containing the 0.5-nmthick water interlayer, following the identical procedure as the previous runs. The water molecules in the interlayer quickly coalesced into a cylindrical shape, but were unable to escape from the DOPC bilayer and the water molecules remained trapped within the lipid bilayer as can be seen from the snap of the simulation after 10 ns shown in Fig. 10. The time-evolution of the local density of the interlayer water molecules along the directions parallel (x- and y-axis) and normal (z-axis) to the
Fig. 12. (a) Snap shot after 10 ns from the NPT simulation of the DOPC bilayer containing 0.5-nm-thick water slab with the DOPC molecules kept at 200 K and water molecules at 300 K, (b) same snap shot as (a) without the DOPC molecules showing the water molecules only.
55
lipid bilayer plotted in Fig. 11 confirms that the interlayer indeed remained trapped inside the lipid bilayer. The lipid molecules extended in the vertical direction in order to accommodate the coalescence of the water molecules as the intermolecular potential interaction was left unchanged. However, the increased stiffness of the hydrocarbon chains prevented the formation of a pore in the lipid layer. Furthermore, to see the effect of the slowing down of the lipid layer molecules relative to the water molecules, a simulation run carried out with the DOPC molecules rescaled at 200 K while maintaining the water molecules at 300 K. As expected, the “cold” DOPC molecules considerably retarded the coalescence of the interlayer water molecules as the lipid layer remained rigid. A simulation run lasting for 10 ns showed that the water molecules remained trapped in the DOPC bilayer as can be seen from Fig. 12. Again the interlayer water density shown in Fig. 13 demonstrates the relatively slow coalescence of the water molecules. The simulation results in Figs. 10–13 evince that the rigidity of the lipid bilayer does have a large effect on the movement of the water
Fig. 13. Spatial (x-axis, y-axis and z-axis) distribution of interlayer water molecules for the DOPC bilayer containing 0.5-nm-thick water slab with the DOPC molecules kept at 200 K and water molecules at 300 K.
56
W.B. Han et al. / Materials Science and Engineering C 36 (2014) 49–56
molecules embedded in the lipid bilayer and the permeability of the lipid membrane. Although a thick flat layer of water molecules may remain stable microscopically as the MD simulation demonstrated, the finite size of the simulation cell limits the long range fluctuation or undulation [29,30]. The long range fluctuation in the bilayer geometry will likely rupture the bilayer and each lipid monolayer would fold back and form micelles or vesicles. A large simulation cell for an extended period would reproduce the transformation of the lipid layer into other geometrical shapes by inserting the water molecules into the hydrophobic region which is equivalent to separating the two lipid monolayers. 5. Conclusion It was demonstrated through MD simulation that a pore structure in a model lipid bilayer was formed when a layer of water was introduced into the hydrophobic regions. Re-arrangement of the water molecules into a cylindrical droplet which locally deformed the lipid monolayer preceded the pore formation. The pore structure formation or local permeability of the lipid membrane is thus closely related to the rigidity of the lipid membrane. Acknowledgment This work was supported by the National Research Foundation of Korea (NRF) grant funded by the Korean government (MEST) (No. 2012000815). Appendix A. Supplementary data The dynamics of the DOPC bilayer having the 0.5, 0.7 and 1.0-nmthick interlayer water slab in the membrane and the plot of the modified Ryckaert–Bellemans potential function can be found online. Supplementary data to this article can be found online at http://dx.doi. org/10.1016/j.msec.2013.11.033.
References [1] C.W. Pratt, K. Cornely, Essential Biochemistry, John Wiley & Sons, New Jersey, 2004. [2] T. Kodama, Y. Tomita, K. Koshiyama, M.J.K. Blomley, Ultrasound Med. Biol. 32 (2006) 905–914. [3] S. Mehier-Humbert, R.H. Guy, Adv. Drug Deliv. Rev. 57 (2005) 733–753. [4] D.P. Tieleman, BMC Biochem. 5 (2004) 10. [5] J. Gehl, Acta Physiol. Scand. 177 (2003) 437–447. [6] S. Bao, B.D. Thrall, D.L. Miller, Ultrasound Med. Biol. 23 (1997) 953–959. [7] W. Tao, J. Wilkinson, E.J. Stanbridge, M.W. Berns, Proc. Natl. Acad. Sci. U. S. A. 84 (1987) 4180–4184. [8] Y. Yamakoshi, T. Miwa, Jpn. J. Appl. Phys. 50 (2011) 07HF01. [9] T. Kodama, N. Tomita, S. Horie, N. Sax, H. Iwasaki, R. Suzuki, K. Maruyama, S. Mori, F. Manabu, J. Electron Microsc. 59 (2010) 187–196. [10] T. Kodama, K. Koshiyama, Y. Tomita, M. Suzuki, T. Yano, S. Fujikawa, Am. Inst. Phys. Conf. Proc. 829 (2006) 34–38. [11] M.J. Haynes, Ultrasound Med. Biol. 26 (2000) 57–62. [12] K. Koshiyama, S. Wada, J. Biomech. 44 (2011) 2053–2058. [13] K. Koshiyama, T. Kodama, T. Yano, S. Fujikawa, Biochim. Biophys. Acta Biomembr. 1778 (2008) 1423–1428. [14] K. Koshiyama, T. Kodama, T. Yano, S. Fujikawa, Biophys. J. 91 (2006) 2198–2205. [15] K. Koshiyama, T. Yano, T. Kodama, Phys. Rev. Lett. 105 (2010) 018105. [16] H.H. An, J.H. Kim, J.H. Lee, D.H. Kwon, H.S. Kim, Y.H. Kim, C.S. Yoon, Chem. Commun. 46 (2010) 9238–9240. [17] A. Terheiden, C. Mayer, K. Moh, B. Stahlmecke, S. Stappert, M. Acet, B. Rellinghaus, Appl. Phys. Lett. 84 (2004) 3891–3893. [18] B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, J. Chem. Theory Comput. 4 (2008) 435–447. [19] D. Poger, A.E. Mark, J. Chem. Theory Comput. 6 (2010) 325–336. [20] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, J. Chem. Phys. 81 (1984) 3684–3690. [21] U. Essmann, L. Perera, M.L. Berkowitz, T. Darden, H. Lee, L.G. Pedersen, J. Chem. Phys. 103 (1995) 8577–8593. [22] B. Hess, H. Bekker, H.J.C. Berendsen, J.G.E.M. Fraaije, J. Comput. Chem. 18 (1997) 1463–1472. [23] W. Humphrey, A. Dalke, K. Schulten, J. Mol. Graph. 14 (1996) 33–38. [24] M. Patra, M. Karttunen, M.T. Hyvonen, E. Falck, P. Lindqvist, I. Vattulainen, Biophys. J. 84 (2003) 3636–3645. [25] S. Esteban-Martin, J. Salgado, Biophys. J. 92 (2007) 903–912. [26] A.H. de Vries, A.E. Mark, S.J. Marrink, J. Phys. Chem. B 108 (2004) 2454–2463. [27] S.A. Pandit, S. Vasudevan, S.W. Chiu, R.J. Mashl, E. Jakobsson, H.L. Scott, Biophys. J. 87 (2004) 1092–1100. [28] J.-P. Ryckaert, A. Bellemans, Chem. Phys. Lett. 30 (1975) 123–125. [29] I. Chandrasekhar, D. Bakowies, A. Glättli, P. Hünenberger, C. Pereira, W.F. van Gunsteren, Mol. Simul. 31 (2005) 543–548. [30] N. Ganesan, B.A. Bauer, T.R. Lucas, S. Patel, M. Taufer, J. Comput. Chem. 15 (2011) 2958–2973.
