10016
Langmuir 2006, 22, 10016-10024
Atomistic Simulations of Langmuir Monolayer Collapse Christian D. Lorenz* and Alex Travesset Department of Physics and Astronomy, Iowa State UniVersity and Ames Laboratory, Ames, Iowa 50011 ReceiVed June 28, 2006. In Final Form: September 9, 2006 Monolayers at the vapor/water interface collapse by exploring the third dimension at sufficient lateral compression, either by forming three-dimensional structures or by solubilization into the aqueous solution. In this paper, we provide an atomistic description of collapse from molecular dynamics (MD) simulations. More specifically, we investigate monolayers of arachidic acids spread on pure water and in an aqueous solution with Ca2+ ions in the subphase. In both cases, it is found that the collapsed systems generally lead to the formation of multilayer structures, which in the system with Ca2+ ions, proceeds by an intermediate regime where the monolayer exhibits significant roughness (of the order of 4 Å). If no roughness is present, the system forms collapsed structures into the aqueous solution. The computational cost of atomic MD limits our simulations to relatively small system sizes, fast compression rates, and temporal scales on the order of a nanosecond. We discuss the issues caused by these limitations and present a detailed discussion of how the collapse regime proceeds at long time scales. We conclude with a summary of the implications of our results for further theoretical and experimental studies.
I. Introduction The phase diagram of insoluble amphiphiles at the vapor/ water interface (Langmuir monolayers) consists of several phases1,2 whose structure and properties, at least for saturated fatty acids, are reasonably well understood.3 At large enough surface pressures (π), monolayers are often found in a crystalline state with a close-packing molecular area (A0). On further increasing the surface pressure, the monolayer undergoes breakage and/or folding in a process referred to as collapse. Monolayer collapse has been extensively investigated in recent years, as it provides an appropriate model system to investigate the mechanics of surface tension regulation in lungs.4-7 Other investigations are motivated by the notion that the underlying principles governing collapse are relevant for understanding a variety of problems that involve the response of amphiphilic molecules to stress, such as membrane fusion and fission in biological membranes.8 The collapse of self-assembled monolayers is also relevant when considering the tribological properties of these monolayers and their use as lubricants, which has been the focus of considerable experimental9-12 and simulation research.13-15 * To whom correspondence should be addressed. E-mail: cdlorenz@ iastate.edu. (1) Gaines, G. In Insoluble Monolayers at the Liquid-Gas Interface, Interscience, New York, 1966. (2) Lundqvist, M. Prog. Chem. Fats Lipids 1978 16, 101. (3) Kaganer, V.; Mo¨hwald, H.; Dutta, P. ReV. Mod. Phys. 1999, 71, 779. (4) Wang, Z. D.; Hall, S. B.; Notter, R. H. J. Lipid Res. 1995, 36, 1283. (5) Robertson, B.; Halliday, H. L. Biochim. Biophys. Acta 1998, 1408, 346. (6) Lipp, M. M.; Lee, K. Y. C.; Takamoto, D. Y.; Zasadzinski, J. A.; Waring, A. J. Phys. ReV. Lett. 1998, 81, 1650. (7) Warriner, H. E.; Ding, J.; Waring, A. J.; Zasadzinski, J. A. Biophys. J. 2002, 82, 835 and references therein. (8) Chernomordik, L. V.; Kozlov, M. M. Annu. ReV. Biochem. 2003, 72, 25. (9) Meyer, E.; Overney, R.; Luthi, R.; Brodbeck, D.; Howald, L.; Frommer, J.; Guntherodt, H. J.; Wolter, O.; Fujihira, M.; Takano, H.; Gotoh, Y. Thin Solid Films 1992, 220, 132. (10) Tsukruk, V. V.; Bliznyuk, V. N.; Hazel, J.; Visser, D.; Everson, M. P. Langmuir 1996, 12, 4840. (11) Major, R. C.; Kim, H. I.; Houston, J. E.; Zhu, X. Y. Tribol. Lett. 2003, 14, 237. (12) Maboudian, R.; Carraro, C. Annu. ReV. Phys. Chem. 2004, 55, 35. (13) Chandross, M.; Webb, E. B., III; Stevens, M. J.; Grest, G. S.; Garofalini, S. H. Phys. ReV. Lett. 2004, 93, 166103. (14) Chandross, M.; Lorenz, C. D.; Grest, G. S.; Stevens, M. J.; Webb, E. B., III. JOM 2005, 57, 55.
Models for monolayer collapse have been proposed on the basis of experimental observations or theoretical models. In ref 16, Ries proposed that the collapse of 2-hydroxytetracosanoic acid films proceeds by the formation of double-layer platelets (bonded at their headgroups), which are initially formed at weakened boundaries emanating from the monolayer toward the vacuum and grow in size until they break and stack on the monolayer to form an inhomogeneous trilayer-layer film. Subsequent imaging with light microscopy and topographic micrographs of transferred collapsed films of the same molecules provided evidence that collapse structures depend on environmental variables such as compression rate or temperature (which controls the phase of the system).17 Other studies of collapse with 10,12-pentacosaiyonic acid, however, suggested that the collapse proceeds by multilayer formation, with no evidence of folding or formation of platelets.18 More recently, experimental studies with tetracosanoic acid provided evidence of a surface roughening collapse and subsequent anisotropic cracking in the form of curved filaments.19 Theoretical descriptions usually assume that monolayer collapse is a long waVelength collective phenomena, which can be described by the use of continuum models. In ref 20, it was proposed that collapse should occur at a slightly negative surface tension, where the monolayer buckles out of the plane with a free energy cost limited by the bending rigidity. The previous model was extended to include spontaneous curvature,21 which should allow buckling more easily, and it was found that collapse occurs at marginally positive surface tensions. In ref 22, a collapse mechanism for mixtures of amphiphiles with different spontaneous curvature was developed. Continuum models for the collapse of solid monolayers have also been described.23 (15) Lorenz, C. D.; Chandross, M.; Grest, G. S.; Stevens, M. J.; Webb, E. B., III. Langmuir 2005, 21, 11744. (16) Ries, H. E., Jr. Nature 1979, 281, 287. (17) Ybert, C.; Lu, W.; Moller, G.; Knobler, C. J. Phys. Chem. B 2002, 106, 2004. (18) Gourier, C.; Knobler, C.; Daillant, J.; Chatenay, D. Langmuir 2002, 18, 9434. (19) Hatta, E. Langmuir 2004, 20, 4059. (20) Milner, S. T.;Joanny, J. F.; Pincus, P. Europhys. Lett. 1988, 9, 452. (21) Hu, J. G.; Granek, R. J. Phys. II (France) 1996, 6, 999. (22) Diamant, H.; Witten, T.; Ege, C.; Gopal, A.; Lee, K. Y. C. Phys. ReV. E 2001, 63, 61602. (23) Saint-Jalmes, A.; Gallet, F. Eur. Phys. J. B 1998, 2, 489.
10.1021/la061868r CCC: $33.50 © 2006 American Chemical Society Published on Web 10/21/2006
Atomistic Simulations of Langmuir Monolayer Collapse
The molecular structure of the amphiphiles plays a fundamental role in monolayer collapse, as evidenced by the different modes of collapse observed for different amphiphiles. Furthermore, the physical parameters at which collapse occurs are directly related to characteristics of the amphiphile at the molecular scale. In fatty acids, for example, the molecular area at collapse is identical to the close packing area of hydrocarbon chains A0 ≈ 19 Å2,2 and the structures formed strongly depend on the charge of the fatty acid headgroup, as demonstrated by the analysis of transferred collapsed fatty acid monolayers with divalent ions in the subphase to a solid support.19,24 These observations clearly point out that in many systems, such as fatty acids, collapse is, to a large extent, a short waVelength phenomena and is a direct consequence of the molecular structure of the amphiphile. In these cases, any theoretical description must incorporate a significant number of the details of the amphiphile at the molecular level. In this paper, we provide, for the first time, an all-atom molecular dynamics (MD) simulation of the collapse process in Langmuir monolayers. The particular system we consider is arachidic acid (AA) in the crystalline phase, both in contact with pure water and with a subphase containing divalent ions (Ca2+). Our theoretical study allows us to investigate the collapse regime with atomic resolution and discuss effects such as the role of water, ion binding, or hydrocarbon chain conformations, which are not usually accounted for within long-wavelength descriptions and, as discussed above, are expected to determine the pathway followed by the collapse process. As shown below, the computational cost of such a study is considerable, thus limiting our investigation to relatively small systems (∼100 AA molecules) and, more significantly, compression rates (∼1 m/s ) 10 Å/ns) that drive the system considerably out of equilibrium. In Section V, we provide a detailed discussion of how our results extrapolate to larger systems and slower compression rates, as well as the implications for experimental results. Previous MD simulations of Langmuir monolayer collapse26 (for phosphatidyl choline monolayers) have been performed with coarse-grained models, while our MD results are united atom simulations. While appropriate for systems where electrostatic charges do not play a critical role (as expected for dipolar lipids such as phosphatidyl choline), the description of charged amphiphiles is subtle within coarse-grained models, and united atoms are therefore more appropriate. Recent MD simulations27 have been focused on the collapse of amphiphilic polymers, which are not Langmuir monolayers, as they have significantly different properties from acyl chains. II. Computational Procedure The collapse of AA (CH3(CH2)17COOH) monolayers is studied with classical MD simulations. The general system studied here consists of the AA monolayer placed between two walls on top of water as shown in Figure 1. Specifically, two different systems have been studied: one in which the AA monolayer is on pure water (shown in Figure 1a) and the other where the AA monolayer is on a CaCl2 solution (shown in Figure 1c). For the system on CaCl2 solution, we make two simplifying assumptions: (1) all of the AA molecules have been deprotonated by the Ca2+ ions25 in solution, resulting in the molecules becoming ionic (CH3(CH2)17COO-) and (2) we only represent the Ca2+ ions in solution since they are the ions which will interact with the now negatively charged AA molecules. The walls in these simulations are modeled such that (24) Kundu, S.; Datta, A.; Hazra, S. Langmuir 2005, 21, 5894. (25) Travesset, A.; Vaknin, D. Europhys. Lett. 2006, 74,181. (26) Nielsen, S. O.; Lopez, C. F.; Moore, P. B.; Shelley, J. C.; Klein, M. L. J. Phys. Chem. B 2003 107, 13911. (27) Leith, D.; Morton-Blake, D. A. Mol. Sim. 2005, 31, 987.
Langmuir, Vol. 22, No. 24, 2006 10017
Figure 1. Snapshots of the initial configurations of the arachidic acid monolayer system consisting of the arachidic acid molecules at 20 Å2/molecule on (a-b) pure water (S4-S5) and (c-d) on Ca2+ solution (S1-S2). (a) and (c) show the initial configurations with all of the monolayer molecules and ions in the water, while (b) and (d) show the initial configurations with half of the monolayer molecules and ions are outside of the water. The water has been represented in these snapshots with dots so that the ions in solution can be seen more clearly. Colors: Arachidic acid, carbon (light blue), hydrogen (silver), and oxygen (red); Water, oxygen (cyan), hydrogen (cyan); Ca2+ (yellow); walls (green). they only interact with the AA molecules but not with the water or the ions in solution. This is the easiest way to simulate the real experimental situation, which are performed on a large trough of water with a small depth so that water or ions are able to flow around the walls without being compressed. The interaction of the walls with the AA molecules is defined with a repulsive Lennard-Jones potential with σ ) 1 (in Å) and ) 1 (in kcal/mol). We have found that the simulation results are not significantly affected by the values for the σ and of the walls (for example, surface pressures are the same). Each wall consists of 748 beads. The surface pressure is computed by summing the forces on each of the wall beads in a given wall and then dividing by its perimeter. Present MD studies are limited to time scales on the order of a nanosecond. These time scales are smaller than relaxation times in Langmuir monolayers, which in some cases reach seconds or even much longer.1 It should therefore be expected that our simulations will be significantly out of equilibrium and the collapse process may depend on the initial conditions. We therefore investigated two initial conditions which were then allowed to relax before starting the collapse process. The first condition consisted of COOH/COOgroups initially placed in the water layer, as shown in Figure 1a and c. The second initial condition is motivated by the hypothesis that the collapse proceeds by first breaking a small domain, which is later ejected during collapse. We therefore considered an initial configuration where half of the monolayer molecules were placed 4 Å above the water surface (Figure 1b). It should be expected that for distances much smaller than 4 Å the initial configuration would be equivalent to S1, while for much larger separations, the AA molecules would already be too far from the interface; thus, 4 Å provides a reasonable value. For the monolayer on the Ca2+ solution, half of the Ca2+ ions are placed above the water surface between the monolayer molecules and the water (Figure 1d). For each of the initial configurations, we studied the collapse for systems where AA molecules were placed at areas per molecule of 20 Å2 and 25 Å2. The dimensions of the system are 69.097 Å × 65.543 Å, which results in 196 AA molecules in the 20 Å2/molecule monolayers and 144 AA molecules in the 25 Å2/molecule monolayers. The water layer is 69.097 Å × 65.543 Å × 55.0 Å. Our system is periodic in the x and y dimensions and finite in the z dimension. To keep any water that may evaporate from the water layer within the box bounds,
10018 Langmuir, Vol. 22, No. 24, 2006
Lorenz and TraVesset
Table 1. System Descriptions system S1 S2 S3 S4 S5
description Ca2+
arachidic acid on solution (22 Å2/molecule). all monolayer molecules and Ca2+ ions are placed in the water layer. arachidic acid on Ca2+ solution (22 Å2/molecule). half the monolayer molecules and Ca2+ ions are placed above the water layer arachidic acid on Ca2+ solution (30 Å2/molecule). all monolayer molecules and Ca2+ ions are placed in the water layer. arachidic acid on pure water (22 Å2/molecule). all monolayer molecules are placed in the water layer. arachidic acid on pure water (22 Å2/molecule). half the monolayer molecules are placed above the water layer.
walls28
we used the standard technique of placing repulsive along the box boundaries in the z dimension (5 Å below the bottom of the water layer and 5 Å above the top of the walls). The repulsive walls were modeled with a Lennard-Jones 9-3 potential
[152 (σr ) - (σr ) ]
E)
9
3
Figure 2. Plot of the number density (no./Å3) and cos(θH2O) as a function of the z dimension for system S4.
(1)
for r < rc, where σ ) 1.0 Å, ) 1.0 kcal/mol and rc ) 0.86 Å. In the systems where the Ca2+ ions are present, the number of ions in the systems is equal to one-half the number of AA molecules, which ensures charge neutrality of the entire system. Each system consists of ∼33 000 atoms. Each of these systems has been summarized in Table 1 and has been given a system number by which they will be referred throughout the remainder of the paper. We use the Charmm 27 force field parameters to model the interactions of the AA molecules, the water molecules, and the Ca2+ ions.29,30 The interaction of the water molecules is described using the TIP3P model.31 Each system was initially equilibrated at constant area and constant temperature for 1 ns. Then, the systems were compressed by moving the walls in the x direction at constant velocities, V/2 and -V/2, which corresponds to a compression rate of V ) 2 m/s. How this compression rate relates to the experimental case is discussed extensively in the Conclusion section. Even at this fast compression rate, the simulation work is considerable, and a single simulation compression takes 10 days on 16 3.2 GHz Intel Xeon processor nodes. Simulations discussed in this paper were performed with the LAMMPS molecular dynamics code.28 In all simulations, the temperature is fixed at 298 K and controlled with a Langevin thermostat with damping applied with a time constant of 0.01 fs-1. The shake algorithm32 was used to constrain the bond lengths of the water molecule, and a 1 fs time step was used. The van der Waals interaction is cut off at 10 Å. A slab version of the particle-particle particle-mesh algorithm33 is used to compute the long-range Coulombic interactions. This version models an isolated slab that is infinitely periodic in the x and the y dimensions and finite in the z dimension.
III. Results The results of our simulations are reported as follows. We discuss the state of the monolayer before collapse in Section (28) Plimpton, S. J. J. Comp. Phys. 1995, 117, 1. (29) Follope, N.; MacKerell, A. D., Jr. J. Comp. Chem. 2000, 21, 86. (30) MacKerell, A. D., Jr.; Bashford, D.; Bellott, M.; Dunbrack, R. L., Jr.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; JosephMcCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos, C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E., III; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiorkiewicz-Kuczera, J.; Yin, D.; Karplus, M. J. Phys. Chem. B 1998, 102, 3586. (31) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (32) Ryckaert, J.-P. ; Ciccotti, G.; Berendsen, H. J. C. J. Comp. Phys. 1977, 23, 327. (33) Crozier, P. S.; Rowley, R. L.; Henderson, D. J. Chem. Phys. 2001, 114, 7513.
Figure 3. Plot of the number density (no./Å3) and cos(θH2O) as a function of the z dimension for system S1.
III.A. The collapse results are reported in Sections III.C and III.B for AA/pure water and the AA/Ca2+ solution systems, respectively. A. Constant Area Simulations. Systems S1 and S4 are studied holding the walls at constant separation. This allows us to investigate the properties of the system before collapse. These simulations were run for 2 ns while data was accumulated and averaged over the last nanosecond. The structure of water was characterized by the number density and the angle (θH2O), defined by the vector connecting the oxygen and the center of the two hydrogens in a given water molecule relative with the vector making up the positive z axis. The water molecule is aligned with its hydrogens pointing in the positive z direction if the cosine is negative and the oxygen pointing in the positive z direction if the cosine is positive. Figure 2 shows a plot of the number density and cos(θH2O) for the AA/pure water system (system S4). In this case, the number density has been defined by the number of atoms per unit volume, not number of molecules per unit volume. Therefore, the expected number density of bulk water is ∼0.1 atoms/Å3, which is the same as ∼0.0333 molecules/ Å3. The plot shows that the number density of water remains the bulk value almost up to the interface. Furthermore, water is randomly aligned in the bulk region (resulting in the net cosine being 0). At the water/AA molecule interface, however, the water is generally aligned with the oxygens pointing in the positive z direction. Therefore, the oxygen in the water molecules are hydrogen bonded with the hydrogens in the -COOH headgroups of the AA molecules. Figure 3 shows a plot of the number density and cos(θH2O) for the AA/Ca2+ system (system S1). The number density of water remains the bulk value almost up to the interface, but the
Atomistic Simulations of Langmuir Monolayer Collapse
Langmuir, Vol. 22, No. 24, 2006 10019
Figure 4. Snapshot of the bound Ca2+ ions to the AA monolayer molecules. The Ca2+ ions and oxygens have been labeled with ‘Ca’ and ‘O’, respectively. Colors: Arachidic acid, carbon (light blue), hydrogen (silver), oxygen (red), and Ca2+ (yellow). Figure 7. Plot of the normalized surface pressure (π/πo) as a function of area per molecule during compression at a rate of V ) 2 m/s for the AA/Ca2+ systems (S1, S2, S3). The value of πo is 4426 mN/m.
Figure 5. Plot of the distribution of the number of AA molecules that are bound to a Ca2+ ion.
Figure 6. Plot of the distribuion of the number of oxygens in AA molecules that are bound to a Ca2+ ion.
orientation of the water molecules at the interface is significantly different from the pure water case. Here the water molecules are aligned so that their hydrogens are pointed toward the AA/water interface, as it should be expected, as the interface is negatively charged. This results in an electric field that partially orders water molecules well into the water layer (∼30 Å) as opposed to just at the interface, as is observed in system S1. We investigated the local structure surrounding the Ca2+ ions in these systems. To do so, we calculated the number of oxygens from both water molecules and from AA molecules in the first hydration shell of the Ca2+ ions. Also, a Ca2+ ion is defined as bound to an AA molecule if an oxygen from the AA headgroup is within the hydration shell of the ion. We find that each Ca2+ ion is bound to 2.5 AA molecules. Figure 4 shows an example of the typical structures of the AA molecule headgroups and the Ca2+ ions. As can be seen, there are two AA molecules bound to one Ca2+ ion and four AA molecules bound to the other Ca2+ ion shown. These are two typical binding patterns observed in this system. Figure 5 shows the distribution of the number of AA molecules bound to each Ca2+ ion, and Figure 6 shows the distribution of the number of oxygens within AA molecules bound to each Ca2+ ion. From these figures, we observe that the
Ca2+ ions bind to 1-6 oxygens from AA molecules and that each ion is bound to 1-4 AA molecules with almost equal probability and a smaller probability for binding to 5 AA molecules. It is found that the first hydration shell of the bound Ca2+ ions includes ∼7 oxygens, 3 from AA molecules (2 from one AA molecule and 1 from another) and 4 from water molecules. We point out that more detailed studies on the structure of interfacial water have been reported for phospholipids, as, for example, reviewed in ref 34, where it is found that interfacial water structures are highly nonuniversal, so we caution that the results obtained in this paper, which are for fatty acids at much smaller molecular area, should not be expected to agree with simulations in phospholipid systems. B. Collapse of AA/Ca2+ Solution Systems. Figure 7 shows a plot of the surface pressure as a function of the area per molecule for the AA/Ca2+ solution systems (systems S1, S2, and S3). The actual values of the surface pressure obtained are unphysical, as will be discussed below, so we present the surface pressure in the y axis in terms of relative units π/πo, where πo ) 4426 mN/m. The general trends of the curves for the different AA/Ca2+ solution systems are similar, showing a sharp increase of the surface pressure as a function of molecular area until an area per molecule of ∼17 Å2/molecule (see Figure 7), where a maximum surface pressure is observed, which is followed by a subsequent drop in the surface pressure. This peak at the maximum of the surface pressure defines the ‘point of failure’ (POF) of the monolayer. System S3 is started with an initial configuration identical to S1, except that at a larger molecular area. Clearly, the surface pressure in the range between 20 and 30 Å2 reflects that the system is out of equilibrium, as the surface pressures are too high and do not coincide with those of system S1, when comparisons become available. We point out, however, that in the collapse regime, defined by molecular areas
