ARTICLE pubs.acs.org/Langmuir
Using Molecular Dynamics to Study Liquid Phase Behavior: Simulations of the Ternary Sodium Laurate/Sodium Oleate/Water System Dylan T. King, Dallas B. Warren, Colin W. Pouton, and David K. Chalmers* Medicinal Chemistry and Drug Action, Monash Institute of Pharmaceutical Sciences, Monash University, Parkville, Victoria, Australia
bS Supporting Information ABSTRACT: The prediction of surfactant phase behavior has applications in a wide range of areas. An accurate modeling of liquid phase behavior can aid our understanding of colloidal process or be used to design phases that respond in a defined way to their environment. In this work, we use molecular dynamics to model the phase behavior of the ternary sodium laurate/sodium oleate/water system and compare the simulation results to experimental data. Simulations were performed with the GROMOS 53A6 united-atom force field and cover the entire ternary phase diagram, producing micellar, hexagonal, and lamellar phases. The aggregate simulation time for the 33 simulations performed during this study is 4.4 μs. We find that the simulations were able to model the experimentally observed liquid phase behavior accurately, showing that the carboxylate and lipid parameters of the 53A6 force field give very good quality results for the in silico prediction of liquid system phase behavior.
1. INTRODUCTION Soaps, salts derived from long-chain carboxylic acids, are classical amphiphiles having a hydrophilic polar headgroup and a hydrophobic tail. Despite having a simple molecular structure, these compounds produce many different phases in water; at lower concentrations, soaps form micelles, intermediate concentrations produce a hexagonal phase, and high concentrations form lamellar structures. Such systems therefore represent a straightforward and useful model for investigating the prediction of liquid phase behavior. The sodium laurate/sodium oleate/ water system has been characterized experimentally by Mongondry et al.,1 and we have used this study as the basis for this work. Here we report the atomistic molecular dynamics (MD) simulations used to study mixtures of the sodium salts of long-chain carboxylic acids lauric acid (sodium dodecanoate) and oleic acid (sodium (Z)-octadec-9-enoate) in aqueous solution (Figure 1). Oleic and lauric acids and their derivatives (phospholipids and glycerides) are important biological compounds, and the study of these carboxylates is of direct relevance to a number of biological processes, particularly lipid digestion. Molecular dynamics studies of phase behavior have become routinely accessible to the research community by using parallel molecular dynamics software such as NAMD,2 GROMACS,3 GROMOS,4 or Amber.5 These packages, coupled with steady increases in computer power and the development of cluster computing, now allow the routine simulation of systems containing hundreds of thousands or even millions of particles. Studies can also be performed at different levels of detail by using an r 2011 American Chemical Society
all-atom force field that explicitly includes all atoms, such as united-atom force fields (nonpolar hydrogen atoms are incorporated into the parent carbon atom), coarse-grained force fields (with groups of multiple atoms represented by a single particle), or even more generalized molecular representations (for example, representing molecules as rigid geometric particles). To model phase phenomena, we need to be able to simulate molecular systems that are significantly larger than the noncovalent structures present in the phase. Smaller structures such as bilayers or micelles can be accessed using united- or all-atom simulations, whereas larger structures such as vesicles are usually modeled using coarse-grained approaches. Simplifying the model by reducing the number of particles represented (e.g., converting from all-atom to united-atom models) greatly increases the computational speed, enabling larger and longer simulations. However, the speedup may come at the cost of reduced accuracy of the simulation and the loss of direct correlation between the chemical structure of the components and their representation in the modeled system. Atomistic and coarse-grained molecular dynamics studies of liquid phase behavior are regularly reported. A large majority of the investigations cover only a single system or a small range of systems and predominately deal with only a single phase or a single phase transition and the characterization of those phases, Received: June 18, 2011 Revised: August 11, 2011 Published: August 12, 2011 11381
dx.doi.org/10.1021/la2022903 | Langmuir 2011, 27, 11381–11393
Langmuir
Figure 1. Stuctures of sodium oleate and sodium laurate showing atom numbering.
for example, the formation of lipid bilayers68 (also see the comprehensive review of lipid bilayer simulations by Marrink et al.9), the conversion between lipid bilayer gels and liquid crystal phases,1016 liquid crystal transitions (isotropicnematic17 and isotropiclamellar18), the formation of micelles1921 and reverse micelles,2224 conversion of micelles to bilayers,25 and transitions from reverse hexagonal to reverse cubic phases.26 Relatively few studies have attempted to conduct broad surveys of phase behavior by using larger numbers of MD simulations to predict more extensive phase diagrams. A limited number of reported examples include those of de Meyer et al., who have performed a study of dimyristoylphosphatidylcholine (DMPC) phospholipid bilayers using coarse-grained models,27 Liu et al., who used coarse-grained dynamics to study the effects of temperature and pressure changes on the shape transformations in vesicles,28 and McDonald and Hanna, who have calculated the pressure/temperature phase diagram of 40 -octyl-[1,10 -biphenyl]-4-carbonitrile (8 CB).29 The molecular dynamics simulation of liquid phase behavior is a general approach with very broad applicability. An accurate prediction of liquid phase behavior provides the potential for designed phases with specific geometrical or chemical properties that respond in a defined manner to changes in temperature, concentration, or other triggers. A notable example of designed phase behavior is the development of rationally designed cubic phases that have been developed for the crystallization of membrane-bound proteins.30 Our specific interest is the design of formulations for poorly water-soluble drugs.3139 Thus, modeling the phase behavior of ionic surfactants will ultimately allow the accurate modeling of the fate of drugs in the digestive tract.
2. METHODS Molecular dynamics simulations were performed using GROMACS3 version 4.0.7. Calculations were performed using the following computer systems: the Victorian Partnership for Advanced Computing (VPAC) Linux cluster (www.vpac.org) and the Victorian Life Sciences Computation Initiative (VLSCI) Linux cluster (www.vlsci.unimelb.edu.au). The VPAC LC consisted of 888 AMD 2.8 GHz dual CPUs connected using Myrinet and the VLSCI LC composed of 1088 Intel Nehalem computer cores running at 2.66 GHz and connected using InfiniBand. Parallel scaling of sample simulations (20 000 steps) was performed to determine the optimum conditions under which to operate the simulations to achieve an efficient use of CPU time. For the 10 nm3 simulations, the optimum utilization was determined to be 16 CPUs for both VPAC and VLSCI, and typical required CPU times were 24 CPU h/ns of simulation on VPAC and 13.7 on the VLSCI system.
ARTICLE
Sodium laurate and sodium oleate were represented using the GROMOSs 53a6 united-atom force field.40 (CH, CH2, and CH3 groups are represented as a single “atom”.) The cis double bond in sodium oleate was modeled using dihedral parameters developed by Barchar et al.41,42 Water was modeled using rigid SPC water43 and constrained using SETTLE.44 The remaining solute bonds were constrained by the LINCS algorithm.45 Isotropic periodic boundary conditions were employed. A cutoff distance of 0.9 nm for both electrostatic and van der Waals interactions and the particle-mesh Ewald (PME) method46 was used for long-range electrostatic interactions (as described in the GROMACS documentation). Temperature coupling used the velocity rescale algorithm47 with a reference temperature of 348 K. Pressure coupling used the Berendsen48 or ParrinelloRahman49 algorithms with a reference pressure of 1 bar and a compressibility of 4.5 105 bar1. All of the ionic surfactant MD simulations were established and performed using the following procedure: (1) The required number of sodium oleate and/or sodium laurate molecules was added using Silico script random_box (http:// silico.sourceforge.net) with a random position and orientation in the periodic cell. Water molecules were randomly added to fill the remaining volume, giving a total number of atoms in the range of 50 000 to 100 000. (2) A steepest-descent minimization of 500 steps was used to remove bad van der Waals contacts between atoms. (3) Three short simulations were used to prepare the system for production simulations: (a) The first simulation was performed without pressure coupling for 5000 steps with a time step of 2 fs. (b) The second simulation introduced Berendsen isotropic pressure coupling with a coupling time constant of 0.1 ps and was run for 10 000 steps with a time step of 2 fs. (c) The third simulation switched to isotropic Parrinello Rahman pressure coupling with a pressure coupling time constant of 0.1 ps and was run for 10 000 steps with a time step of 2 fs. (4) From this third coordinate file, the full MD simulation was commenced. A time step of 5 fs was used, along with a pressure reference of 1 bar, a pressure coupling time constant of 2 ps, a reference temperature of 348 K, and a temperature coupling time constant of 0.1 ps. The aggregate simulation time for the 33 simulations performed during this study is 4.4 μs. Molecular aggregation was analyzed using the Silico script find_aggregate, which combines molecules into aggregates by comparing distances between carbon atoms. Two molecules are considered to be part of the same aggregate if they have carbon atoms separated by a distance of less than 0.4 nm. The potential of mean force w(r) as a function of the distance r between atoms was calculated using the relationship w(r) = kT ln(p(r)), where k is Boltzmann’s constant, T is the absolute temperature, and p(r) is the probability density along that coordinate.50 An estimate of the probability density was obtained by calculating the radial distribution function with GROMACS utility g_rdf over the last 20 ns of each simulation. The sampling of the RDF is improved by including all surfactant molecules in the RDF calculation. A visualization of the simulation trajectories was performed using software packages VMD51 and PYMOL (http://www.pymol.org). Images for publication were produced with PYMOL. The fraction of trans dihedral angles for each of the simulations was calculated using the GROMACS program g_angle.
3. RESULTS AND DISCUSSION In this work, we have used molecular simulation to model the phase behavior of the ternary sodium laurate/sodium oleate/ water system. The phase behavior of such simple molecules is experimentally well characterized, and the principal factors 11382
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 2. Spontaneous structuring of surfactant/water mixtures over 100 ns of simulation time. (A) Formation of micelles in 20% w/w sodium laurate/ water. (B) Formation of a lamellar phase in 80% w/w sodium laurate/water. (C) Formation of micelles in 20% w/w sodium oleate/water. (D) Formation of a crystalline lamellar phase in 80% w/w sodium oleate/water. Atom coloring: Na+, dark blue; oxygen, red; and carbon, cyan (sodium laurate) or green (sodium oleate). Water molecules have been omitted for clarity.
dictating phase structure are well understood. The phase structure, whether it is micelle, hexagonal, cubic, or lamellar phase, is driven by the surfactant volume fraction and the surfactant packing parameters that control the curvature of the wateroil phase interface.52 The sodium laurate/sodium oleate/water system is an excellent model for investigating simulation and analysis methods that might be applied to more complex surfactant systems. This work extends our previous methods, which have been applied to the modeling of bile salt aggregation38 and glyceride lipid formulations,32 to a larger and more comprehensive study of a three-component phase system. Briefly, the modeling approach uses molecular dynamics simulations to generate well-equilibrated colloidal systems by starting from an unbiased, random arrangement of molecular components. This evolution is illustrated in Figure 2, showing the rapid and spontaneous formation of micelles
Table 1. Sodium Laurate/Water Systems Used to Investigate the Influence of Periodic Cell Size cell size no. (nm)
no. of laurate
no. of water
molecules
molecules
composition simulation (% w/w)
time (ns)
1
5
269
828
80.03
160
2
10
2152
6624
80.03
160
3
20
17216
52992
80.03
160
or bilayers from random sodium oleate/water and sodium laurate/ water mixtures over a simulation time of 100 ns. The rate at which the system reaches equilibrium is dependent on the amount of surfactant present and the resulting viscosity of the system. Simulations containing higher concentrations of surfactant equilibrate more slowly, with the simulation times required to reach 11383
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 3. Influence of periodic cell size on the simulation of sodium laurate, 80% w/w in water. (A) Radial distribution functions for the sodium laurate O1O1 distance: 5 (black), 10 (red), and 20 nm (green) cells. (B) Fraction of dihedrals in the trans orientation vs the dihedral number: 5 (b), 10 (9), and 20 nm (3) cells.
Table 2. Summary of the Simulations Used to Model the Ternary Sodium Laurate/Sodium Oleate/Water Phase System
a
number sodium
number sodium
number
actual %
simulation
median aggregate sizea/
final
laurate
oleate
water
w/w
time (ns)
number of aggregates
state
no
system
4
laurate 0.7%
18
31 000
0.72
100
3/2
small aggregate and isolated molecules
5
laurate 5%
136
31 272
5.1
100
14/8
spherical micelles and smaller aggregates
6
laurate 10%
272
29 360
10.3
100
27/11
spherical micelles
7
laurate 20%
536
26 496
20.0
100
33/17
spherical and prolate micelles
8 9
laurate 30% laurate 42%
808 1072
23 184 18 096
30.1 42.2
100 200
70/11 75/8
prolate micelles prolate and worm micelles
10
laurate 50%
1344
16 560
50.0
200
/1
worm micelles
11
laurate 55%
1472
14 850
55.0
100
/1
worm micelles
12
laurate 63%
1608
11 592
63.1
200
/1
worm micelles and lamellae
13
laurate 70%
1880
9936
70.0
200
/1
worm micelles and lamellae
14
laurate 80%
2152
6624
80.0
200
/1
lamellae
15
laurate 90%
2416
3312
90.0
200
/1
lamellae
16 17
oleate 5% oleate 10%
96 200
31 272 29 360
4.9 10.3
200 100
26.5/4 34/6
spherical micelles spherical micelles
18
oleate 20%
400
26 496
20.3
100
44/8
spherical and prolate micelles
19
oleate 30%
584
23 184
29.9
100
194/3
spherical, prolate and worm micelles
20
oleate 43%
800
18 096
42.8
100
/1
worm micelles
21
oleate 50%
984
16 560
50.1
100
/1
worm micelles
22
oleate 55%
1075
14 850
55.0
100
/1
worm micelles
23
oleate 62%
1272
13 248
61.7
100
/1
incomplete lamellae
24 25
oleate 70% oleate 80%
1368 1568
9936 6624
69.9 80.0
300 200
/1 /1
lamellae lamellae
26
oleate 90%
1760
3312
90.0
300
/1
lamellae
27
mixture 15%
201
146
28 050
15.0
100
38/9
spherical micelles
28
mixture 45%
602
440
18 150
45.0
100
/1
worm micelles
29
mixture 50%
669
488
16 500
50.0
100
/1
worm micelles
30
mixture 60%
803
586
13 200
60.0
100
/1
worm micelles
31
mixture 73%
970
708
9075
72.5
200
/1
lamellae
32 33
mixture 80% mixture 90%
1070 1204
782 879
6600 3300
80.0 90.0
100 200
/1 /1
lamellae lamellae
Aggregate numbers and median values exclude any isolated molecules.
equilibrium ranging from 100 to 300 ns. Plots of the system potential energy show that the more concentrated systems, such as 80% w/w sodium laurate, become stable after approximately 170 ns (Supporting Information Figure S1). The final equilibrated systems then give significant insight into the nature of the
colloidal structures that are generated and the intermolecular interactions that drive the formation of these structures. 3.1. Periodic Cell Size and Its Effect on Simulation Results. An important issue in the simulation of continuous phases is that of overall simulation size. What size system is necessary to 11384
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 4. Final frames of nine selected simulations of the sodium laurate/water system showing the progression from spherical micelles (simulations 5 and 6) through elongated micelles (simulations 7 and 8) and worm micelles (simulations 912) to lamellar phases (simulations 13 and 14). Simulation numbering is described in Table 1. Atom coloring: carbon atoms colored by aggregate, oxygen atoms are red, and sodium ions are dark blue. Water atoms have been omitted.
provide useful results? Larger systems enable the modeling of larger colloidal structures and should reduce packing and other artifacts imposed by periodic boundary conditions but are more computationally expensive; doubling the cell dimensions results in an 8-fold increase in the number of atoms. Periodic boundary conditions are well known to cause artifacts; the size of the structures formed cannot exceed the size of the cell and is constrained to integer repeat units in each of the cell dimensions. For example, in the 100 ns snapshots shown in Figure 2, sodium laurate bilayers in Figure 2B repeat three times in the Y dimension but only once in the X dimension, and sodium oleate bilayers in Figure 2D repeat twice in each of the X and Y dimensions. Furthermore, the behavior of each molecule can be perturbed by interacting with their own periodic images, either directly or indirectly (mediated through other molecules), with the potential of long-range electrostatic interactions causing significant perturbations. For example, Weber et al. found that the size of the simulation periodic cell strongly influenced the unfolding energies of charged peptides,53 where the unfolded length of the peptide approached or exceeded the edge length of the periodic cell. Before commencing the simulations that were the focus of this study, we investigated the influence of cell size on a sample system by performing three 160 ns simulations of an 80% w/w sodium laurate/water system using 5, 10, and 20 nm cubic cells. Full system specifications are shown in Table 1. The concentration of 80% w/w sodium laurate was chosen because it lies within the lamellar region of the phase diagram.1 An inspection of the
final frames showed that after 160 ns the 10 nm cell produced a strongly ordered system containing three planar bilayers and a single pool of water. The 20 nm cell also contained multiple lamellar regions and water pools. In contrast, the 5 nm cell produced a highly curved bilayer, inconsistent with the results of the 10 and 20 nm simulations. An analysis of the simulations revealed more measurable differences between the 5 nm and the 10/20 nm simulations. Figure 3A shows the radial distribution function (RDF) for laurate O1O1 interactions in the surfactant headgroups. The 5 nm cell simulation differs substantially from the other two, indicating that the small cell significantly influences the intermolecular interactions, in particular, the absence of splitting of the peak corresponding to the second coordination shell (∼0.45 nm). Additional RDFs are shown in Supporting Information Figure S2. The conformational preferences of the individual molecules are also affected, as shown in Figure 3B. The number of trans bonds in the 5 nm system is noticeably reduced compared to the number in the 10 and 20 nm cells. Accordingly, the 10 nm cubic cells were utilized for subsequent simulations. 3.2. Modeling of the Phase Behavior. If molecular simulations of surfactants (or any other system) are to be useful, they require good-quality molecular mechanics parameters that accurately model both internal molecular conformations and intermolecular interactions (surfactantsurfactant, surfactantsolvent, and solventsolvent). In this work, we have used the recently developed GROMOS 53A6 united-atom force field40 that has been parametrized to reproduce free energies of solvation in water and cyclohexane and has been used extensively to model proteins, micelles, 11385
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 5. Final frames from nine selected simulations of the sodium oleate/water system showing a progression from spherical micelles (simulations 16 and 17) through elongated and worm micelles (simulations 1823) to lamellar phases (simulations 24 and 25). Simulation numbering is described in Table 1. Atom coloring: carbon atoms are colored by aggregate, oxygen atoms are red, and sodium ions are dark blue. Water atoms have been omitted.
and membranes. To increase the simulation speed, we have used a united-atom force field with a rigid water model and bond length constraints, allowing a computationally efficient simulation time step of 5 fs.54 All simulations have been performed at a temperature of 348 K (75 °C) and a pressure of 1 atm (1 bar), corresponding to the experimentally determined phase diagram reported by Mongondry et al.1 Simulations utilized a 10 nm cubic periodic cell that, for these systems, provides an appropriate trade-off between periodic cell artifacts and CPU time. To cover the entire phase space of the ternary sodium laurate/ sodium oleate/water system, three sets of molecular dynamics simulations were performed. The first two sets investigated the binary sodium oleate/water and sodium laurate/water systems, covering the concentration range of 0.7 to 90% w/w sodium laurate and 5 to 90% w/w sodium oleate. A third, smaller series investigates the ternary sodium laurate/sodium oleate/water system with sodium laurate and sodium oleate maintained in a 1:1 mass ratio. Table 2 records the composition of each simulation, the total simulation time, and observations of the number of colloidal aggregates formed with a visual assessment of phases resulting at the end of the simulation. 3.2.1. Binary Systems Sodium Laurate/Water and Sodium Oleate/Water. Figures 4 and 5 show the molecular structures present in the final frames of several of the sodium laurate/water and sodium oleate/water simulations, with the surfactant molecules colored to distinguish separate surfactant aggregates. Two
molecules are defined as belonging to a common aggregate if they are two carbon atoms within 0.4 nm; that is, they are making at least one hydrophobic contact. Using this approach, individual aggregates can be identified at concentrations below approximately 3040% w/w surfactant. At higher concentrations, the molecular structures become continuous by this definition, and separate aggregates are not distinguished. The simulations clearly show progressive changes with increasing surfactant concentration. At the lowest surfactant concentration modeled, 0.7% w/w laurate (32 mM), which is near the reported critical micelle concentration55 of 30 mM, only very small aggregates were present at the end of the simulation. This simulation may not be fully equilibrated after 100 ns and may require multiple interaction events between separate aggregates before an equilibrated state is reached. Distinct individual aggregates can be readily identified at concentrations of between 5 and 40% w/w surfactant. In the 5% w/w simulations of both laurate and oleate, approximately spherical micelles form with diameters of ∼3 nm. As expected, the hydrophobic tails are buried in the core of the micelle whereas the headgroup and one or two methylene groups of the alkyl tail project from the micelle into the surrounding solvent. Many of the negatively charged carboxylate headgroups are closely associated with the sodium counterions. The micelles are dynamic, with a constant internal diffusion of surfactant molecules and an occasional exchange of surfactant molecules between micelles. The sodium laurate micelles contain a 11386
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 6. (A) Ordered packing in a lamellar region of 90% w/w sodium laurate with water. Two layers of sodium ions are present, colored dark and light blue. (B) The same region showing carboxylate groups. COO groups in the top layer are shown by thick bonds with green carbon atoms; those in the lower layer have yellow carbon atoms. The laurate alkyl chains have been omitted. Close contacts between sodium and oxygen atoms are shown as lines. (C) Detail showing the packing arrangement of sodium ions and carboxylate groups from the simulation. Atom coloring is the same as for B. (D) Packing arrangement of sodium ions and carboxylate groups in 2-oxooctanoate crystals.56
Figure 7. Final frames from six selected simulations of the sodium laurate/sodium oleate/water ternary system showing a progression from spherical micelles (simulation 27) through mixed micelles and worm micelles (simulations 28 and 30) to well-ordered lamellar phases (simulations 3133). Simulation numbering is described in Table 1. Atom coloring: carbon atoms in sodium laurate are green and those in sodium oleate are blue, oxygen atoms are red, and sodium ions are dark blue. Water atoms have been omitted. 11387
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 8. Radial distribution functions, g(r), for binary mixtures of (A, C, E, G) sodium laurate/water and (B, D, F, and H) sodium oleate/water systems. Line coloring for % w/w surfactant: 5 (black), 10 (red), 20 (light green), 30 (yellow), 40 (dark blue), 50 (pink), ∼60 (light blue), 70 (gray), 80 (brown), and 90 (dark green).
median of 14 molecules, and in sodium oleate systems the micelles have a median size of 27 molecules (Table 2). As the concentration is increased to 10% w/w, both the oleate and laurate micelles grow while remaining approximately spherical. Once the concentration reaches 20% w/w, some of the micelles are observed to become prolate, and at 30 (oleate) or 40% w/w (laurate), extended worm micelles form. The worm micelles have diameters of about 3 nm and are long enough to span the 10 nm periodic cell. In the concentration ranges of 4355% w/w oleate and 5055% w/w laurate, only worm micelles are present.
In both the oleate and laurate simulations, a transition to lamellar phases starts to occur above concentrations of 55% w/w and exclusively lamellar arrangements exist once the surfactant concentration reaches 70% w/w. At concentrations of between 55 and 70% w/w, worm micelles, flattened worm micelles (or incomplete lamellae), and lamellar regions are observed. The formation of these mixed structures is potentially governed by the incomplete equilibration of the molecular system and the small size of the periodic cell relative to the size of the colloidal structures being formed. It is likely that the small cell containing a limited number of surfactant molecules, combined with the 11388
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 9. Gibbs energies of atom pair association, w(rmin) (kJ mol1) for sodium laurate and sodium oleate simulations as a function of surfactant concentration: first minimum (b), second minimum (O), third minimum (1), fourth minimum (2), and fifth minimum (9). Blocks of color denote the experimentally observed formation of micelles (red), the hexagonal phase (green), the lamellar phase (blue), and crystals (black).1 Intermediate white regions consist of a mixture of the two adjoining phases.
restricted simulation times, does not allow proper partitioning into separate phases. At high surfactant concentrations, water is expelled from between the surfactant layers and carboxylate groups and sodium ions pack directly together in a double layer, driven by the strong electrostatic interaction. The excluded water then forms pools within the simulation cell. At the lamellar interface, the formation of wellorganized domains where sodium ions and surfactant carboxylate groups pack in highly regular cubic arrays is observed (Figure 6). Figure 6A shows the regular arrangement of sodium ions in the
close-packed regions of sodium laurate. Figure 6B shows the regular ordering of carboxylate groups within the ionic lattice, with each sodium ion coordinated to six oxygen atoms belonging to five carboxylate groups, with one carboxylate group providing two coordinating oxygen atoms. The remaining four carboxylate groups provide one oxygen atom each (Figure 6C). A search of the Cambridge crystallographic database reveals that although there are no reported crystal structures of simple fatty acid sodium salts, crystal structures of sodium 2-oxooctanoate (code NOXCAP)56 and sodium hydrogen α-ketoglutarate (code COTPEG)57 exhibit 11389
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
Figure 10. Average fraction of trans CC bonds in the sodium laurate/ water (top) and sodium oleate/water simulations (bottom). The bond number is marked on the right-hand side of each curve. Blocks of color denote the experimentally observed formation of micelles (red), the hexagonal phase (green), the lamellar phase (blue), and crystals (black).1 Intermediate white regions consist of a mixture of the two adjoining phases. Some bonds are omitted from the sodium oleate graph for clarity.
lamellar crystal packing with sodiumcarboxylate interactions that are very similar to those observed in our simulations (Figure 6D). 3.2.2. Ternary System Sodium Laurate/Sodium Oleate/Water. Following our simulations of the binary systems, we investigated ternary sodium laurate/sodium oleate/water by running simulations containing equal masses of the two surfactants and varying amounts of water (simulations 2733, Table 2). Figure 7 shows the final frames of these simulations with the oleate and laurate molecules colored differently to distinguish them. A very similar progression to the previous cases was observed, from spherical micelles to worm micelles and a sharp transition to lamellar phases. 3.3. Analysis of Radial Distribution Functions. The radial distribution function, g(r) or RDF, describes the relative probability of finding one atom a given distance from another. This probability is directly related to the free energy of interaction between the two particles. Consequently, RDFs computed from the simulations are sensitive tools for the analysis of simulated systems. Figure 8 shows RDFs for key atom pairs in the sodium oleate/water and sodium laurate/water systems: carboxylate oxygen atoms (O1O1), carboxylate and water oxygen atoms (O1OW), sodium ions (NaNa), and terminal-alkane-chain carbon atoms (C12C12 for laurate or C18C18 for oleate). The sodium oleate and sodium laurate RDFs are quite similar, showing only minor differences in peak shape and peak magnitude (Supporting Information Table S1). It is readily apparent that all of the RDFs change with surfactant concentration, with
ARTICLE
the most striking changes belonging to the O1O1 and NaNa distribution functions. These functions change markedly upon the formation of lamellar phases. In the micellar phases (below ∼50% w/w surfactant), the O1O1 RDFs (Figure 8A,B) display broad maxima at about 0.46 and 0.65 nm and a low probability in the range of 0.20.4 nm. The peak at 0.65 nm reflects the general intercarboxylate spacing on the micelle surface. An examination of the simulation structures shows that the small peak at 0.43 nm is caused by the close approach of two carboxylate groups promoted by an intermediary sodium ion. The low probability of close approach reflects the electrostatic repulsion of the carboxylate headgroups. However, once the surfactant concentration is above ∼70% w/w and the system becomes lamellar, three well-defined peaks appear in the laurate O1O1 RDF at short distances of 0.31, 0.43, and 0.46 nm. Two peaks at distances of 0.31 and 0.44 nm become apparent in the oleate RDFs. An examination of the molecular structures (Figure 6B,C) reveals that the 0.31 nm peak results from pairs of oxygen atoms that are bonded to one or more common sodium atoms; for example, oxygen pairs ac, be, and ef in Figure 6C are separated by distances ranging from 0.31 to 0.36 nm. The peak at 0.44 nm is due to oxygen atoms bonded to a single sodium ion, with the three atoms making an angle of about 120° (pairs ad and bf). The 0.46 nm peak arises from oxygen atoms bonded to sodium, making an angle of close to 160° (pairs df and ec). The separate 0.44 and 0.46 nm peaks, present in laurate RDF, are not apparent in the oleate RDFs because the lamellar oleate structures are less ordered. Like the O1O1 RDFs, the NaNa RDF (Figure 8E,F) also undergoes a dramatic change on transition from micellar to lamellar phases in both oleate and laurate simulations. The lamellar phases show a dominant peak at an interatomic distance of 0.36 nm, corresponding to the distance between adjacent sodium ions (pair st). A shoulder at 0.42 corresponds to pair rt (Figure 6C). The carboxylatewater oxygen RDFs (O1OW) for both oleate and laurate show the expected first, second, and third solvation shells at distances of 0.3, 0.5, and 0.6 nm (Figure 8C,D). The RDFs for the terminal methyl groups (C12C12 or C18C18) are broad envelopes, reflecting the nonspecific nature of the hydrophobic interaction between alkyl chains and the fluid nature of the alkyl terminal region. Neither the O1OW or terminal methyl RDFs change shape with surfactant concentration, but they do change in intensity. These changes are discussed further below. 3.4. Analysis of Phase Change. Our key area of interest is the ability of MD simulations to predict the phase of a system at a given point on the phase diagram. Performing multiple series of MD simulations traversing the sodium laurate/sodium oleate/ water phase diagram allows a comparison of the changes in molecular structure and dynamic behavior across these series. An analysis of the chemical properties and the 3D structure of the equilibrated system allows the detection of phase changes. By using the RDFs shown in Figure 8, potentials of mean force (PMFs) have been derived for the interaction of the key atoms (Supporting Information Figures S3 and S4). The PMFs describe the average Gibbs energy (G) for the interaction of those atom pairs, with the Gibbs energy at large separation being zero. Figure 9 plots the changes in Gibbs energy with concentration at local energy minima in each PMF w(rmin); these distances correspond to maxima in the RDFs shown in Figure 8. The most striking feature of these PMF curves is the large change that occurs on transition from the micellar/hexagonal to the lamellar phase. This change is observed in both systems and is in agreement with the experimental data. 11390
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
Figure 11. Comparison of experimental1 and molecular dynamics-derived phase diagrams of sodium laurate, sodium oleate, and water at 348 K (75 °C). Single phase regions are shaded gray, and intermediate white regions consist of a mixture of the two adjoining phases. L1 denotes the micellar, H1 denotes the hexagonal, and Lα denotes the lamellar phase. Points indicate the compositions of each simulation and the phase observed in the final frame: isolated surfactant molecules and small aggregates (f), spherical micelles (b), prolate micelles (b), wormy micelles (2), and lamellae (9). The coloring represents micelles (red), the hexagonal phase (green), and the lamellar phase (blue).
The O1OW and O1O1 coordination numbers were calculated using RDFs of the binary systems by numerically integrating g(r) multiplied by the atom number density with respect to the radius (Supporting Information Figure S5). The number of molecules in the first water solvation shell around the carboxylate oxygens remains constant at approximately five until the experimentally observed formation of the lamellar phase for both laurate and oleate (Supporting Information Figures S5A and S5B). In the case of the O1O1 coordination number, a distinct change in the coordination number with concentration is observed near the experimentally observed transition from the hexagonal to the lamellar phase (Supporting Information Figures S5C and S5D) Figure 10 shows the shows the fraction of trans dihedrals in the sodium laurate, sodium oleate, and mixture simulations. This data was calculated using the final 20 ns of each simulation. In each case, an increase in the fraction of trans dihedrals is observed, coinciding with the experimentally observed phase transition from hexagonal to lamellar. This is an expected trend because the lamellar phase is a more ordered structure with the alkyl chains aligning themselves with each other. Figure 11 shows the experimental phase diagram of the sodium laurate/sodium oleate/water system at 348 K (75 °C)1 overlaid with the phases observed at the end of each simulation (Table 2). In general, there is very good agreement between the predicted and observed phases; observed prolate and spherical micelles are present in the micellar region of the plot, wormy micelles are observed in the hexagonal phase regions, and lamellar structures are found in the lamellar region. The simulations perform particularly well in the
micellar region, accurately modeling the differences in the propensity for micelle formation of sodium laurate and sodium oleate; the simulations show that sodium oleate forms micelles up to a maximum concentration of 20% w/w whereas some micelles are still observed at 40% w/w sodium laurate. The transition between the hexagonal and lamellar phases is also predicted in the correct region although the experimentally observed differences between oleate and laurate are not as clear. The hexagonal phase regions of the diagram proved to be the least well modeled portion in all three sets of simulations. Extended wormy micelles were observed in each case, rather than a well-ordered array of rods. The disorder in this phase region may be due to the influence of the cubic periodic cell on the formation of the internal structures—the rod diameter is approximately 30% of the cell size. Additionally, the transition from wormy micelles to a well-ordered hexagonal phase requires the breakage and reassembly of the wormy micelles and is expected to be slow.
4. SUMMARY AND CONCLUSIONS In this study, we have performed a series of simulations modeling the phase behavior of the ternary sodium laurate/sodium oleate/ water system. These simulations accurately model the experimentally observed phase behavior of this system. Molecular simulation is a useful and reasonably widely used approach for investigation of the molecular interactions that govern phase behavior; however, the large majority of published simulations do not attempt to cover entire phase diagrams, generally being interested in a single point or a small region. This study demonstrates that a relatively 11391
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir straightforward, although computationally intensive, approach provides a good-quality model of the phase diagram and insight into the molecular processes of phase formation. Such predictive, molecularscale models of phase processes are an important tool in many areas of research, for example, the modeling of biological surfactants in cell wall lipids or digestive processes and the study of synthetic systems such as drug formulations. In the reverse sense, the ability to predict phase behavior accurately is also an important test of current molecular modeling methods. This study uses the “off-the-shelf” GROMOS 53A6 united-atom force field, which was developed to model liquid phase behavior, particularly for biological systems. The parametrization of the force field is critical to the correct modeling of the system, and in this case, the simulations recapitulate the important features of the phase diagram. The comparison between experimental and theoretical results shows that the carboxylate and lipid parameters of the 53A6 force field give very good quality results, even though the united-atom force field omits the hydrogen atoms attached to carbon. A large body of phase data is available in the literature, and in general, we expect comparisons such as this one between the phase data and theoretical studies to be generally applicable to force field validation. Having investigated this surfactant system, we intend to extend this modeling approach to the investigation of systems relevant to drug formulation and absorption.
’ ASSOCIATED CONTENT
bS Supporting Information. Additional figures and tables as described in the text. GROMACS format topology files for sodium laurate and sodium oleate. Final frames for each of the 33 simulations listed in Table 2. This material is available free of charge via the Internet at http://pubs.acs.org. ’ AUTHOR INFORMATION Corresponding Author
*Phone: +61 3 9903 9110. E-mail:
[email protected].
’ ACKNOWLEDGMENT We thank Dr. Hassan Benameur for his encouragement of this work and financial support from Capsugel. This work was also generously supported by the Victorian Life Sciences Computational Initiative (VLSCI), which provided computer time, systems support, and a scholarship for Dylan King, and the Victorian Partnership for Advanced Computing (VPAC), which provided computer time and systems support. ’ REFERENCES (1) Mongondry, P.; Macosko, C. W.; Moaddel, T. Rheol. Acta 2006, 45, 891. (2) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. J. Comput. Chem. 2005, 26, 1781. (3) Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. J Chem. Theory Comput. 2008, 4, 435. (4) Christen, M.; Hunenberger, P. H.; Bakowies, D.; Baron, R.; Burgi, R.; Geerke, D. P.; Heinz, T. N.; Kastenholz, M. A.; Krautler, V.; Oostenbrink, C.; Peter, C.; Trzesniak, D.; van Gunsteren, W. F. J. Comput. Chem. 2005, 26, 1719. (5) Case, D. A.; Cheatham, T. E., III; Darden, T.; Gohlke, H.; Luo, R.; Merz, K. M., Jr.; Onufriev, A.; Simmerling, C.; Wang, B.; Woods, R. J. J. Comput. Chem. 2005, 26, 1668.
ARTICLE
(6) Kranenburg, M.; Smit, B. J. Phys. Chem. B 2005, 109, 6553. (7) Hogberg, C. J.; Lyubartsev, A. P. J. Phys. Chem. B 2006, 110, 14326. (8) Zidar, J.; Merzel, F.; Hodoscek, M.; Rebolj, K.; Sepcic, K.; Macek, P.; Janezic, D. J. Phys. Chem. B 2009, 113, 15795. (9) Marrink, S. J.; de Vries, A. H.; Tieleman, D. P. Biochim. Biophys. Acta 2009, 1788, 149. (10) Faller, R.; Marrink, S. J. Langmuir 2004, 20, 7686. (11) Cascales, J. J. L.; Otero, T. F.; Romero, A. J. F.; Camacho, L. Langmuir 2006, 22, 5818. (12) Leekumjorn, S.; Sum, A. K. Biochim. Biophys. Acta, Biomembr. 2007, 1768, 354. (13) Qin, S. S.; Yu, Z. W.; Yu, Y. X. J. Phys. Chem. B 2009, 113, 8114. (14) Coppock, P. S.; Kindt, J. T. Langmuir 2009, 25, 352. (15) Coppock, P. S.; Kindt, J. T. J. Phys. Chem. B 2010, 114, 11468. (16) Horta, B. A. C.; de Vries, A. H.; Hunenberger, P. H. J. Chem. Theory Comput. 2010, 6, 2488. (17) Tiberio, G.; Muccioli, L.; Berardi, R.; Zannoni, C. ChemPhysChem 2009, 10, 125. (18) Lintuvuori, J. S.; Wilson, M. R. J. Chem. Phys. 2010, 132, 105103. (19) Maillet, J. B.; Lachet, V.; Coveney, P. V. Phys. Chem. Chem. Phys. 1999, 1, 5277. (20) Yakovlev, D. S.; Boek, E. S. Langmuir 2007, 23, 6588. (21) Turner, D. C.; Yin, F. C.; Kindt, J. T.; Zhang, H. L. Langmuir 2010, 26, 4687. (22) Rosenfeld, D. E.; Schmuttenmaer, C. A. J. Phys. Chem. B 2006, 110, 14304. (23) Sammalkorpi, M.; Karttunen, M.; Haataja, M. J. Phys. Chem. B 2007, 111, 11722. (24) Chowdhary, J.; Ladanyi, B. M. J. Phys. Chem. B 2009, 113, 15029. (25) Schuler, L. D.; Walde, P.; Luisi, P. L.; van Gunsteren, W. F. Eur. Biophys. J. 2001, 30, 330. (26) Marrink, S. J.; Tieleman, D. P. Biophys. J. 2002, 83, 2386. (27) de Meyer, F. J.; Benjamini, A.; Rodgers, J. M.; Misteli, Y.; Smit, B. J. Phys. Chem. B 2010, 114, 10451. (28) Cao, Z.; Lin, Z.; Wang, J.; Liu, H. J. Comput. Chem. 2009, 30, 645. (29) McDonald, A. J.; Hanna, S. J. Chem. Phys. 2006, 124, 164906. (30) Misquitta, Y.; Cherezov, V.; Havas, F.; Patterson, S.; Mohan, J. M.; Wells, A. J.; Hart, D. J.; Caffrey, M. J. Struct. Biol. 2004, 148, 169. (31) Pouton, C. W. Eur. J. Pharm. Sci. 2006, 29, 278. (32) Warren, D. B.; Chalmers, D. K.; Hutchison, K.; Dang, W. B.; Pouton, C. W. Colloids Surf., A 2006, 280, 182. (33) Cuine, J. F.; Charman, W. N.; Pouton, C. W.; Edwards, G. A.; Porter, C. J. Pharm. Res. 2007, 24, 748. (34) Cuine, J. F.; McEvoy, C. L.; Charman, W. N.; Pouton, C. W.; Edwards, G. A.; Benameur, H.; Porter, C. J. J. Pharm. Sci. 2008, 97, 995. (35) Porter, C. J.; Pouton, C. W.; Cuine, J. F.; Charman, W. N. Adv. Drug Delivery Rev. 2008, 60, 673. (36) Pouton, C. W.; Porter, C. J. Adv. Drug Delivery Rev. 2008, 60, 625. (37) Mohsin, K.; Long, M. A.; Pouton, C. W. J. Pharm. Sci. 2009, 98, 3582. (38) Warren, D. B.; Chalmers, D. K.; Pouton, C. W. Mol. Pharmaceutics 2009, 6, 604. (39) Warren, D. B.; Benameur, H.; Porter, C. J.; Pouton, C. W. J. Drug Targeting 2010, 18, 704. (40) Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. J. Comput. Chem. 2004, 25, 1656. (41) Bachar, M.; Brunelle, P.; Tieleman, D. P.; Rauk, A. J. Phys. Chem. B 2004, 108, 7170. (42) Martinez-Seara, H.; Rog, T.; Karttunen, M.; Reigada, R.; Vattulainen, I. J. Chem. Phys. 2008, 129, 105103. (43) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; Hermans, J. In Intermolecular Forces; Pullman, B., Ed.; Reidel: Dordrecht, The Netherlands, 1981; p 331. 11392
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
Langmuir
ARTICLE
(44) Miyamoto, S.; Kollman, P. A. J. Comput. Chem. 1992, 13, 952. (45) Hess, B.; Bekker, H.; Berendsen, H. J. C.; Fraaije, J. G. E. M. J. Comput. Chem. 1997, 18, 1463. (46) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. J. Chem. Phys. 1995, 103, 8577. (47) Bussi, G.; Donadio, D.; Parrinello, M. J. Chem. Phys. 2007, 126, 014101. (48) Berendsen, H. J. C.; Postma, J. P. M.; Vangunsteren, W. F.; Dinola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684. (49) Parrinello, M.; Rahman, A. J. Appl. Phys. 1981, 52, 7182. (50) Kirkwood, J. G. J. Chem. Phys. 1935, 3, 300. (51) Humphrey, W.; Dalke, A.; Schulten, K. J. Mol. Graphics 1996, 14, 33. (52) Mitchell, D. J.; Tiddy, G. J. T.; Waring, L.; Bostock, T.; Mcdonald, M. P. J. Chem. Soc., Faraday Trans. 1 1983, 79, 975. (53) Weber, W.; Hunenberger, P. H.; McCammon, J. A. J. Phys. Chem. B 2000, 104, 3668. (54) Marrink, S. J.; Tieleman, D. P. Biophys. J. 2002, 83, 2386. (55) Stanley, F.; Warner, A.; Schneiderman, E.; Stalcup, A. J. Chromatogr., A 2009, 1216, 8431. (56) Tavale, S. S.; Pant, L. M.; Biswas, A. B. Acta Crystallogr. 1964, 17, 215. (57) Lis, T.; Matuszewski, J. Acta Crystallogr., Sect. C: Cryst. Struct. Commun. 1984, 40, 2016.
11393
dx.doi.org/10.1021/la2022903 |Langmuir 2011, 27, 11381–11393
