| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||

* Division of Mechanical and Space Engineering, Hokkaido University, Sapporo, Japan; and
Biomedical Engineering Research Organization, Tohoku University, Sendai, Japan
Correspondence: Address reprint requests to Shigeo Fujikawa, E-mail: fujikawa{at}eng.hokudai.ac.jp.
| ABSTRACT |
|---|
|
|
|---|
| INTRODUCTION |
|---|
|
|
|---|
The biological cell membranes are, as well known, composed of phospholipids, proteins, and carbohydrates, and models describing its structure continue to be revised and refined, beginning from the fluid mosaic model (8
) to the lipid rafts microdomain models (9
,10
). Although the cell membranes in general have such complex structures, the fundamental and common element is the phospholipid bilayer. The method of MD simulations for cell membranes has been developed for the lipid bilayer accordingly. There are a number of important contributions about the numerical methodology, such as the force field (11
), the charges of polar molecules (12
), the effect of choice of ensembles (13
), the treatment of the electrostatic interactions (14
), and so on. On the basis of these contributions, works on the molecular transport across the bilayer are being published (15
18
). The simulations of pore formation in lipid bilayer have also been performed recently (19
22
).
However, we emphasize that much effort has been devoted to the research of an almost static configuration or slow translocation of lipid molecules. Dynamical processes of cell membranes under the action of shock waves have not been investigated with MD simulations so far, to our knowledge. The shock wave is a high pressure wave with a steep wave front that propagates at a supersonic speed, and it passes the cell membrane within a very short time of the order of picoseconds. Therefore, understanding a high-speed phenomenon induced by the interaction of shock waves with cell membranes should be indispensable to the advance of the gene/drug delivery technique with shock waves. We address this challenging problem in this work using unsteady nonequilibrium MD simulations.
| METHODS |
|---|
|
|
|---|
Lipid bilayer/water system
The simulations here have first been performed with a small bilayer/water system consisting of 32 dipalmitoylphosphatidylcholine (DPPC) lipids and 4800 water molecules in a rectangular simulation box. Thereafter, we have validated the results by conducting further simulations with a larger system of 128 lipids and 19,200 water molecules. This section explains the simulation method of the equilibrium bilayer/water system, which was utilized for an initial condition of shock wave simulation described in subsequent sections.
The refined united atom force field (23
) with AMBER99 force field (24
) was applied to the DPPC lipids. The atomic charges calculated by Chiu et al. (12
) were used for DPPC, and they are illustrated in Fig. 1 with the molecular structure and AMBER atom types. The simulation procedure was substantially the same as an already established one (23
). For water, the simple point charge model (25
) was used. Lennard-Jones interactions were cut off at 1.0 nm, and long-range electrostatic interactions were handled by means of the particle mesh Ewald method (26
) with Ewald coefficient 0.275.
|
T = 0.5 ps, T = 323 K,
p = 1.0 ps, and p = 101.3 kPa. The equations of motion of molecules were integrated employing the leap frog method with time step
= 2.0 fs, and the three-dimensional periodic boundary conditions were applied. After an
8-ns run of equilibration simulation with the SHAKE algorithm (28
Here, we remark that our bilayer/water system included a large water layer of thickness
14 nm, whereas the thickness of the bilayer was
4 nm. This is because this simulation of shock waves required a large number of water molecules, as explained below. Such a large water layer has also been used by Marrink et al. in an MD simulation of lipid adhesion (29
).
Because we are interested in the dynamical process of structural change in bilayer resulting from shock wave irradiation, it may be better to remove the constraint of molecular bond lengths based on the SHAKE algorithm. We therefore continued the equilibrium simulation of the bilayer/water system without using SHAKE and confirmed that the equilibrium configuration of lipid and water molecules remained plausible at least up to 12 ns without SHAKE, as compared with not only earlier numerical results (11
,13
,14
,23
) but also experimental ones (30
). For example, the average interfacial areas per lipid molecule were 67 Å2 in the larger system and 64 Å2 in the smaller one, and they are included in the range of experimental values from 56 Å2 to 72 Å2 reported in a review article (30
). The obtained equilibrium bilayer/water system was utilized for the preparation of the initial configuration for the shock wave simulations.
Modeling of shock waves
In the previous experiment (6
), it was shown that the transport of fluorophores across the membranes of leukemia cells is governed by the shock impulse rather than the maximum pressure. The shock impulse per unit area I is defined as
![]() | (1) |
From the definition of impulse, the shock impulse I can be regarded as the increment of momentum of water divided by an area A on which the pressure p(t) is exerted. That is,
![]() | (2) |
The shock waves in the experiment were pulses with duration times of at least 180 ns, corresponding to the pulse width of
270 µm in water. It is impossible or prohibitively expensive to reproduce such shock waves in the simulation box with a periodic boundary condition. In this study, we numerically expressed the shock impulse by the momentum M(t+) of water molecules. More precisely, at the beginning of shock wave simulation, we added the momentum M(t+) = I x A to water molecules in a volume A x LZ adjacent to the bilayer, where the area A is the cross-sectional area normal to the z direction of the simulation box, e.g., A = 7.18 x 5.99 nm = 43.01 nm2 in the larger system, and LZ is the length of the volume in the z direction (see Fig. 2). The choice of LZ is arbitrary, and we put LZ = 4 nm, which is almost equal to the initial thickness of bilayers, because this simulation was focused on the behavior of bilayers with the excess momentum M(t+) added by shock wave. The momentum change of water molecules at the beginning of shock simulation was numerically implemented by the addition of an average velocity V to the thermal velocity of water molecules in the equilibrated bilayer/water system, as shown in Fig. 2. The average velocity V is given by
![]() | (3) |
|
= 2.0 fs to 0.78 fs according to the size of V to avoid the excess approach of molecules with large velocities. The AMBER 7 set of programs (31
Owing to the periodic boundary conditions, the simulations should be terminated at the time when the effect of the shock impulse reached the boundary at the opposite side of the simulation box. To simulate the essential part of structural changes of the bilayer before the termination of the simulation, sufficiently thick water regions were required in both sides of the bilayer in the simulation box. We therefore used the bilayer/water system with a large water layer of thickness
14 nm.
| ANALYSIS |
|---|
|
|
|---|
In the collapse and rebound stages of bilayers, the ordering in the tails (acyl chains) of lipid molecules was considerably different from that in the equilibrium state. The order parameter SCD (11
,14
,32
) is widely known to be a useful measure for the ordering of lipid molecules in equilibrium states, and usually SCD is evaluated in a long-time average. In the analysis of collapsing and rebounding bilayers, the application of long-time averages should be bypassed, and therefore we used an averaged instantaneous order parameter
CD defined as
![]() | (4) |
i is the angle between the axis of the ith molecular axis and the bilayer normal (the z axis) and NC (= 28) is the number of carbons in both sn-1 and sn-2 chains (see Fig. 1). The angle
is were evaluated from the instantaneous configurations of lipid molecules. By definition, in the equilibrium states,
CD is equal to the average of the usual SCD over all carbons in acyl chains, the value of which is estimated as
0.16 from an earlier numerical result (11
The fluidity of lipid molecules within the bilayer plane (the xy plane) has often been studied by means of the diffusion coefficient in the equilibrium MD simulations (33
). However, the diffusion coefficient was not an appropriate measure of the fluidity in a high-speed phenomenon in the unsteady and nonequilibrium state. To quantify the lateral fluidity of lipids in the bilayer plane in unsteady states, we introduced an accumulated lateral displacement ALD defined as
![]() | (5) |
t (= 8
) is a time interval for producing a configuration frame, k is the frame number, and n = t/
t. Note that ALD is independent of the choice of the origin of the coordinate system, and it increases with increase in the movement of lipid molecules in the bilayer. We also examined the growth rate of ALD, defined as dALD/dt, which should be a constant in an equilibrium state.
The transport of water molecules into the bilayer was the most important phenomenon in the results of these simulations. It has been suggested experimentally (34
,35
) and confirmed numerically (11
,13
) that very few water molecules exist beyond the carbonyl group in the sn-1 chain (see Fig. 1) in equilibrium states. Accordingly, when a water molecule entered a hydrophobic region between the carbonyl groups in sn-1 chains in upper and lower sides of bilayers in the simulation, we decided that the water molecule was transported into the bilayer.
| RESULTS AND DISCUSSION |
|---|
|
|
|---|
350 fs from the instant of shock impulse impact, and then it rebounded slowly. The numerical result in the smaller system completely agrees with that in the larger system.
|
|
CD, because the chain length becomes small if the disorder of chain bend angles increases. In Fig. 5, we present the temporal changes of
CDs of the upper and lower layers separately for the case of I = 50 mPa·s. Before the beginning of collapse, both
CDs were equal to 0.16 as expected from knowledge of the equilibrium state (11
CD of the upper layer and then that of the lower layer began to decrease, and the time lag was
0.1 ps, which corresponded to the time that the shock impulse passed the bilayer. The two
CDs were rapidly reduced until reaching their minimum values, and the following rebound process was relatively slow. From the comparison with Fig. 4, it is immediately realized that the change of bilayer thickness was due to the change of chain bend angles.
|
|
0.1 nm/ps (see Fig. 6 a). However, despite the nonzero growth rate, very few water molecules enter the hydrophobic region in the equilibrium state, as mentioned before. Therefore, the continuance of high fluidity, i.e., the large growth rate for some period, seems to play an essential role for the penetration of water molecules because the growth rate can be regarded as a measure of the lateral fluidity of lipid molecules.
Water penetration into the hydrophobic region
The application of shock impulse to bilayer induced the penetration of water molecules in the hydrophobic region of the bilayer. In the following, we demonstrate how the water molecules penetrated into the bilayer and discuss this important phenomenon.
Fig. 7 is a series of snapshots for the same result as that shown in Fig. 3 (in a slight close-up). In the figure, acyl chains are eliminated and water molecules are indicated by blue spots. The penetration of water molecules into the hydrophobic region of the bilayer is shown clearly. Note that, as can be seen from Fig. 7 a, even in the equilibrium state (before the arrival of the shock impulse), several water molecules exist near the lipid headgroups shown in yellow. However, they are outside the hydrophobic region between the carbonyl groups in sn-1 chains in upper and lower sides of the bilayer, and they cannot enter the hydrophobic region without the action of the shock impulse in the timescale of MD simulations.
|
|
2.5 molecules/ps·nm2 in the case of I = 50 mPa·s and it exceeded 7.5 molecules/ps·nm2 in the case of I = 100 mPa·s. We remark that the trend observed in Fig. 8 b is similar to that of the growth rate of ALD shown in Fig. 6 b. That is, as the shock impulse increases, the lateral fluidity of lipid molecules increased and so did the number of penetrated water molecules.
As shown in Fig. 8 b, the number of penetrated water molecules increased with increasing shock impulse. This qualitatively agrees with the previous experiment (6
), although there were a few differences: in the experiment, calcein (622 Da) and fluorescein isothiocyanate dextran (71.6 kDa) were delivered into the cells by the shock impulse of 100 Pa·s, whereas in the simulation, water (18 Da) penetrated into the bilayer by the impulse of at most 100 mPa·s. This may be a reflection of the fact that the transport of large molecules across cell membranes requires a large amount of impulse.
The mass density profile of water molecules in the case of I = 50 mPa·s is illustrated in Fig. 9, where Fig. 9 a shows the initial profile and the shading indicates the hydrophobic region. From the figure, it is clear that water molecules are carried into the hydrophobic region (see Fig. 7) and the shrinkage and recovery of the hydrophobic region are induced by high density water. An empty region behind the high density water (12 nm < z < 15 nm) is caused by the velocity difference in the initial condition.
|
| SUMMARY |
|---|
|
|
|---|
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
Submitted on January 18, 2006; accepted for publication June 13, 2006.
| REFERENCES |
|---|
|
|
|---|
2. Doukas, A. G., D. J. McAuliffe, and T. J. Flotte. 1993. Biological effects of laser-induced shock waves: structural and functional cell damage in vitro. Ultrasound Med. Biol. 19:137146.[CrossRef][Medline]
3. Gambihler, S., M. Delius, and J. W. Ellwart. 1994. Permeabilization of the plasma membrane of L1210 mouse leukemia cells using lithotripter shock waves. J. Membr. Biol. 141:267275.[Medline]
4. Delius, A., and G. Adams. 1999. Shock wave permeabilization with ribosome inactivating proteins: a new approach to tumor therapy. Cancer Res. 59:52275232.
5. Rosenthal, I., J. Z. Sostaric, and P. Riesz. 2004. Sonodynamic therapya review of the synergistic effects of drugs and ultrasound. Ultrason. Sonochem. 11:349363.[Medline]
6. Kodama, T., M. R. Hamblin, and A. G. Doukas. 2000. Cytoplasmic molecular delivery with shock waves: importance of impulse. Biophys. J. 79:18211832.
7. Kodama, T., A. G. Doukas, and M. R. Hamblin. 2002. Shock wave-mediated molecular delivery into cells. Biochim. Biophys. Acta. 1542:186194.[Medline]
8. Singer, S. J., and G. L. Nicholson. 1972. The fluid mosaic model of the structure of cell membranes. Science. 175:720731.
9. Simons, K., and E. Ikonen. 1997. Functional rafts in cell membranes. Nature. 387:569572.[CrossRef][Medline]
10. Edidin, M. 2003. The state of lipid rafts: from model membranes to cells. Annu. Rev. Biophys. Biomol. Struct. 32:257283.[CrossRef][Medline]
11. Berger, O., O. Edholm, and F. Jahnig. 1997. Molecular dynamics simulations of a fluid bilayer of dipalmitoylphosphatidylcholine at full hydration, constant pressure, and constant temperature. Biophys. J. 72:20022013.
12. Chiu, S. W., M. Clark, V. Balaji, S. Subramaniam, H. L. Scott, and E. Jakobsson. 1995. Incorporation of surface tension into molecular dynamics simulation of an interface: a fluid phase lipid bilayer membrane. Biophys. J. 69:12301245.
13. Shinoda, W., N. Namiki, and S. Okazaki. 1997. Molecular dynamics study of a lipid bilayer: convergence, structure, and long-time dynamics. J. Chem. Phys. 106:57315743.[CrossRef]
14. Patra, M., M. Karttunen, M. T. Hyvonen, E. Falck, P. Lindqvist, and I. Vattulainen. 2003. Molecular dynamics simulations of lipid bilayers: major artifacts due to truncating electrostatic interactions. Biophys. J. 84:36363645.
15. Tieleman, D. P., and J. Bentz. 2002. Molecular dynamics simulation of the evolution of hydrophobic defects in one monolayer of a phosphatidylcholine bilayer: relevance for membrane fusion mechanisms. Biophys. J. 83:15011510.
16. Zhu, F., E. Tajkhorshid, and K. Schulten. 2002. Pressure-induced water transport in membrane channels studied by molecular dynamics. Biophys. J. 83:154160.
17. Ulander, J., and A. D. J. Haymet. 2003. Permeation across hydrated DPPC lipid bilayers: simulation of the titrable amphiphilic drug valproic acid. Biophys. J. 85:34753484.
18. Bemporad, D., C. Luttmann, and J. W. Essex. 2004. Computer simulation of small molecule permeation across a lipid bilayer: dependence on bilayer properties and solute volume, size, and cross-sectional area. Biophys. J. 87:113.
19. Tieleman, D. P. 2004. The molecular basis of electroporation. BMC Biochem. 5:10:112.[CrossRef][Medline]
20. Leontiadou, H., A. E. Mark, and S.-J. Marrink. 2004. Molecular dynamics simulations of hydrophilic pores in lipid bilayers. Biophys. J. 86:21562164.
21. Tarek, M. 2005. Membrane electroporation: a molecular dynamics simulation. Biophys. J. 88:40454053.
22. Gurtovenko, A. A., and I. Vattulainen. 2005. Pore formation coupled to ion transport through lipid membranes as induced by transmembrane ionic charge imbalance: atomistic molecular dynamics study. J. Am. Chem. Soc. 127:1757017571.[CrossRef][Medline]
23. Smondyrev, A. M., and M. L. Berkowitz. 1999. United atom force field for phospholipid membranes: constant pressure molecular dynamics simulation of dipalmitoylphosphatidicholine/water system. J. Comput. Chem. 20:531545.[CrossRef]
24. Wang, J. M., P. Cieplak, and P. A. Kollman. 2000. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 21:10491074.[CrossRef]
25. Berendsen, H. J. C., J. P. M. Postma, W. F. van Gunsteren, and J. Hermans. 1981. Interaction models for water in relation to protein hydration. In Intermolecular Forces. B. Pullman, editor. Reidel Publishing Company, Dordrecht, The Netherlands. 331342.
26. Essmann, U., L. Perera, M. L. Berkowitz, T. Darden, H. Lee, and L. G. Pedersen. 1995. A smooth particle mesh Ewald method. J. Chem. Phys. 103:85778593.[CrossRef]
27. Berendsen, H. J. C., J. P. M. Postma, W. F. van Gunsteren, A. DiNola, and J. R. Haak. 1984. Molecular dynamics with coupling to an external bath. J. Chem. Phys. 81:36843690.[CrossRef]
28. Ryckaert, J. P., G. Ciccotti, and H. J. C. Berendsen. 1977. Numerical integration of the Cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes. J. Comput. Phys. 23:327341.[CrossRef]
29. Marrink, S.-J., O. Berger, P. Tieleman, and F. Jähnig. 1998. Adhesion forces of lipids in a phospholipid membrane studied by molecular dynamics simulations. Biophys. J. 74:931943.
30. Nagle, J. F., and S. Tristram-Nagle. 2000. Structure of lipid bilayers. Biochim. Biophys. Acta. 1469:159195.[Medline]
31. Pearlman, D. A., D. A. Case, J. W. Caldwell, W. S. Ross, T. E. Cheatham, S. Debolt, D. Ferguson, G. Siebel, and P. Kollman. 1995. AMBER, a package of computer-programs for applying molecular mechanics, normal-mode analysis, molecular dynamics and free-energy calculations to simulate the structural and energetic properties of molecules. Comput. Phys. Commun. 91:141.[CrossRef]
32. Egberts, E., and H. J. C. Berendsen. 1988. Molecular dynamics simulation of a smectic liquid crystal with atomic detail. J. Chem. Phys. 89:37183732.[CrossRef]
33. Almeida, P. F. F., and W. L. C. Vaz. 1995. Lateral diffusion in membranes. In Handbook of Biological Physics, Vol. 1. Structure and Dynamics of Membrane. R. Lipowsky and E. Sackmann, editors. Elsevier, Amsterdam. 305357.
34. Zaccai, G., J. K. Blasie, and B. P. Schoenborn. 1975. Neutron diffraction studies on the location of water in lecithin bilayer model membranes. Proc. Natl. Acad. Sci. USA. 72:376380.
35. Büldt, G., H. U. Gally, and J. Seelig. 1979. Neutron diffraction studies on phosphatidylcholine model membranes. J. Mol. Biol. 134:673691.[CrossRef][Medline]
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |