Subscriber access provided by UNIV MASSACHUSETTS BOSTON

Article

Coarse-Grained Molecular Dynamics Simulation of Self-Assembly and Surface Adsorption of Ionic Surfactants Using an Implicit Water Model Shihu Wang, and Ronald G. Larson Langmuir, Just Accepted Manuscript • Publication Date (Web): 07 Jan 2015 Downloaded from http://pubs.acs.org on January 7, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Langmuir is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

Coarse-Grained Molecular Dynamics Simulation of Self-Assembly and Surface Adsorption of Ionic Surfactants Using an Implicit Water Model Shihu Wang and Ronald G. Larson* Department of Chemical Engineering, University of Michigan 2300 Hayward Street, Ann Arbor, Michigan 48109-2136 ABSTRACT We perform coarse-grained molecular dynamics simulations for sodium dodecyl sulfate (SDS) surfactant using a modification of the Dry Martini force field (Arnarez et al. 2014) with implicit water. After inclusion of particle mesh Ewald (PME) electrostatics, an artificially high dielectric constant for water (εr = 150), and re-parameterization, we obtain structural and thermodynamic properties of SDS micelles that are close to those obtained from the standard Martini force field with explicit water, which in turn matches those of atomistic simulations. The gains in computational efficiency obtained by removing explicit water allow direct simulations of the self-assembly of SDS in solution. We observe surfactant exchange among micelles and micelle fission and fusion and obtain realistic, equilibrated micelle size distributions at modest computational cost, as well as a transition to cylindrical micelles at high surfactant concentration or with added salt. We further apply this parameterized force field to study the adsorption of SDS onto hydrophobic surfaces and calculate the adsorption kinetics and equilibrium adsorption isotherm. The greatly increased speed of computation of surfactant self-assembly made possible by this Dry Martini method should allow future simulation of competitive adsorption of multiple surfactant species to surfaces, as well as simulation of micellar shape transitions. INTRODUCTION

ACS Paragon Plus 1 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Sodium dodecyl sulfate (SDS) is widely used in many commercial products, such as toothpastes, shampoos, soaps etc. It is also used for DNA extraction and polyacrylamide gel electrophoresis for protein separation. Due to its common use, SDS surfactant has been extensively studied experimentally over the past century via small angle neutron scattering,1,2 small angle X-ray scattering,3 light scattering,4–7 calorimetry and voltammetry8 and others. At 25°C, SDS self-assembles first into spherical micelles above its first critical micelle concentration (CMC) at ~8mM and into rod-like micelles above its second CMC at ~70mM.4,8 The morphology and size of SDS aggregates and SDS phase behavior are also strongly affected by other parameters, such as temperature,2,5,6 and salt concentrations1,5,6 and counterions.1,7 In addition to experimental studies, understanding of SDS physical chemistry has been advanced tremendously by computer simulations. Many different atomistic force fields (FFs), such as GROMOS,9,10 OPLS,10 CHARMM,11,12 and AMBER,13 have been applied to the study of SDS. The sizes, shape, conformation, and internal structures of SDS micelles of different aggregate numbers have been calculated and compared with experimental results.10–12,14–16 The restricted penetration of water into sulfate dodecyl micelle leaving a 1.2nm water-free hydrocarbon core and disruption of water-water hydrogen bonding networks near SDS headgroups,17 the distribution of ions around sulfate dodecyl micelles,15,17 and influences of the types of ions, such as Li+ and NH4+, on sulfate dodecyl micelles15 have also been determined computationally. Most of the FFs mentioned above produce quite similar results for the overall structures of micelles with aggregation numbers 60 or 100, such as the radial distribution of SDS dodecyl tail.22 However, some minor differences, for example, ion distribution around the micelle, have also been observed and were found to strongly affect micellar structures with aggregation number of 300 or higher.10 An extensive comparison of the influence of different FFs on SDS aggregate structure has recently been published.10 Additionally, many previous atomistic simulations were performed using a

ACS Paragon Plus 2 Environment

Page 2 of 31

Page 3 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

single pre-assembled micelle, either a spherical micelle with an aggregation number around 609,11,14–16 or elongated micelles with much larger aggregation numbers.10,12 Direct atomistic simulation of SDS self-assembly into multiple micelles involves simulation of millions of atoms for microseconds of time, and hence incurs high computational cost,18 although cost might be mitigated by using a united atom (UA) model19 or a thermodynamic integration method.12 To mitigate these limitations, many recent surfactant simulations have resorted to coarsegrained (CG) models, where several atoms (and their associated hydrogens) are grouped into one CG bead,20–24 thereby making larger spatial and temporal scales computationally accessible. Examples include the SDK CG model23,25 which uses roughly 3 heavy atoms for each CG bead (i.e., a 3:1 mapping), and the Martini CG model which mainly employs a 4:1 mapping.21,26 While CG FFs are computationally faster than atomistic FFs, simulations of self-assembly still involve a large number of degrees of freedom because of the number of waters contained in a relatively large simulation box required to observe this process. In addition, micelle fission and fusion are too slow to fully equilibrate even with CG FFs at a reasonable time scale.24 For example, for the Martini FF, only aggregates of discrete sizes rather than a continuous size distribution were observed.33 These results indicate that a much longer time is required to successfully simulate the self-assembly process and to obtain the equilibrated micelle size distribution. Since removing explicit water reduces considerably the degrees of freedom, a further enhancement in computational speed is possible by using implicit water models.27–29 Panagiotopoulos and coworkers have parameterized an implicit solvent model for SDS by focusing on matching the CMC and the micellar aggregation number at CMC to experimental values.27 In their model, SDS sulfate group is treated as one CG bead while each tail CH2 group and the terminal CH3 are modeled using united atom (UA) model. An effective ion interaction that considers the hydration layering and a distance-dependent electric permittivity are used. To match

ACS Paragon Plus 3 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

both the experimentally observed CMC and micelle aggregation number, the solvophobic attraction is concentrated only on the terminal tail bead, while all other tail beads still interact with an LJ potential corresponding to that of the OPLS-UA force field. Using modified potentials and grand canonical Monte Carlo simulations they obtained a realistic CMC (9.3±0.3mM) and a realistic average aggregate size and standard deviation about the average (57±6) at the CMC at 298K. Panagiotopoulos et al. also showed that this implicit water model reduces CPU time by two orders of magnitude relative to the UA model with SPC water.27 However, its transferability to other systems, such as interactions of SDS surfactants with other molecules in multicomponent mixtures, might be limited because of the complicated potential functions and the non-standard coarse-graining methods used, i.e. only the sulfate group was treated as a CG bead while the dodecyl tail was described by an UA model. Recently Marrink et al. developed Dry Martini FF with implicit water for biomolecules.28 Similar to the standard Martini FF,26 Dry Martini is also a CG force field and considers several different types of CG beads. The non-bonded interactions among different bead types are described by a simple 6-12 Lennard-Jones (LJ) potential and Coulombic interactions if the beads are charged. The LJ potentials for these different bead types are parameterized in the Dry Martini FF using different levels of interactions to match experimental data such as the partition free energy between water and organic solvents. In the Dry Martini FF, each bead type is similar to a segment in a chemical compound and thus serves as a building block for a molecule. With these different bead types in Dry Martini, building CG models for new molecules is relatively easy and the parameters are usually transferable. In fact, the Martini FF with explicit water was originally developed for lipids,26,30 and was later extended to carbohydrates,31 proteins32 and polymers,33–35 and a similar strategy may work for Dry Martini.

ACS Paragon Plus 4 Environment

Page 4 of 31

Page 5 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

In this paper, we parameterize the SDS surfactant by selecting appropriate bead types in the Dry Martini FF, add particle mesh Ewald (PME) electrostatics with a high dielectric constant (εr=150) instead of a shifted electrostatic cutoff, and apply our parameterization to simulate directly SDS self-assembly and surface adsorption. We show that the Dry Martini FF greatly increases the speed of computation. As a result, surfactant self-assembly, surfactant exchange between micelles, and micelle fission and fusion are observed and the equilibrium cluster size distribution obtained. We also calculate the equilibrium adsorption of SDS surfactants onto hydrocarbon surfaces at several surfactant concentrations. We note that because of the computational speed enhancement produced by use of an implicit solvent, the parameters designed for simulating SDS micellization can also be used to simulate competitive adsorption of multiple surfactant species to hydrocarbon surfaces as well as micellar shape transitions, both of which we explore here. As discussed above, these parameters are also transferable and thus should make possible future simulations of SDS surfactants with other molecules developed using Dry Martini FF. SIMULATION DETAILS We perform the following simulations using the Dry Martini FF recently developed by Marrink et al. for biomolecules.28 Dry Martini FF is very similar to the standard Martini FF.26 It maps roughly four heavy atoms to one CG bead. Five main types of beads including neutral, polar, apolar, charged, and a special bead (for example, for ring compounds) are considered. Each particle type has a number of subtypes to fine-tune the interactions, such as hydrogen-bonding capabilities and degree of polarity. The non-bonded interactions among the different bead types are described by a simple 6-12 Lennard-Jones (LJ) potential and Coulombic interaction if the beads are charged. Different levels of LJ potentials are assigned to each bead type pair to match thermodynamic experimental data. The bonded interactions are modeled by weak harmonic

ACS Paragon Plus 5 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

stretching and bending potentials. Relative to standard Martini, the potentials in Dry Martini are adjusted to take into account the implicit water.

Figure 1: a) Mapping scheme for dodecyl sulfate in the Dry Martini FF: sulfate group corresponds to one CG bead of type Qa and 4 carbon tail atoms to one CG bead of type C1. C1 is sub-scripted with numbers that run sequentially from the terminal C13 bead to the C11 bead that is attached to the head group. b) Tabulated parameters for non-bonded and bonded potentials for different beads in the Dry Martini FF. The non-bonded LJ interaction potentials between each bead type pair are taken from the original Dry Martini FF28 and the bonded parameters are assigned to match the RDF of SDS micelles as discussed in the text. Sodium dodecyl sulfate (SDS) consists of a 12-carbon tail, a sulfate group and a sodium ion. As shown in Figure 1a, the mapping of atomistic model to CG model in Dry Martini is the same as that in standard Martini.21 Specifically, the 12-carbon tail is modeled using three C1 beads; Qa bead with a charge -1 is selected to represent the sulfate group. The hydrated chloride ion in Dry Martini is represented by a Qa bead with a charge of -1 and hydrated sodium ions by Qd beads each with a charge of 1. The parameters for the non-bonded LJ potentials between

ACS Paragon Plus 6 Environment

Page 6 of 31

Page 7 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

different bead types from Dry Martini FF28 are listed in Table 1 for reference. The bonded parameters for bond and angle potentials are selected here in order to reproduce the radial distribution function of SDS micelles from the standard Martini FF as shown below. All simulations and analyses in the following are performed with the GROMACS 4.5.636 simulation package. Simulations using standard Martini FF are performed in an NPT ensemble; the temperature is controlled by an Berendsen thermostat37 and pressure is held at 1atm by an Berendsen barostat37 both with a coupling constant of 4ps. Simulations using the Dry Martini FF are performed in an NVT ensemble with temperature coupled to an Berendsen thermostat37 with a coupling constant of 4ps. In both FFs, the LJ potential is smoothly shifted to zero between 9Å and 12Å. Please note that in standard Martini and Dry Martini, the Coulomb potential is shifted smoothly to zero between 0 and 12Å to achieve high simulation speeds. However, as discussed below, we find that the shifted cutoff scheme leads to unrealistic clustering of SDS (and sodium octyl sulfate, SOS) aggregates. So, in our simulations, Particle Mesh Ewald (PME) summation38 is adopted with a relative permittivity (εr) equal to 80 for simulations with standard Martini and εr=150 with Dry Martini. The real-space cutoff is 1.2nm and the Fourier grid spacing is 0.20nm in both FFs. It is worth mentioning that radial distribution function (RDF) of SDS micelles calculated using PME with εr=80 agrees well with that using original SDS Martini FF (with shifted cutoff),21 as shown in Figure S1. Even though εr=150 seems unreasonably high, RDF of SDS micelles calculated using Dry Martini is also close to that calculated using Martini with shifted cutoff and εr=15 or Martini with PME and εr=80. More importantly, using εr=150 with Dry Martini also removes clustering of surfactant aggregates that is observed with either the shifted cutoff scheme or PME with smaller values of εr such as 15 or 80. The second order stochastic dynamics integrator in GROMCS is used for Dry Martini. We use a time step of 40 femtoseconds for both FFs.

ACS Paragon Plus 7 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

In this study, we consider different scenarios, including pre-assembled SDS micelles, selfassembly of SDS (with and without salts) and adsorption of SDS onto hydrocarbon slabs. All systems are charged neutral. The preassembled SDS micelle is composed of 60 SDS surfactants in a box of 9nm×9nm×9nm and is simulated using both standard Martini and Dry Martini at 300K. The structural properties of SDS micelle, such as radial distribution functions and micelle size, and the potential of mean force (PMF) of an SDS molecule pulled from the micelle are compared to validate the Dry Martini FF parameters given in Figure 1. After these comparisons, we investigate the self-assembly of SDS surfactants from a randomly distributed initial configurations using Dry Martini. Two different SDS surfactant concentrations, namely 72.8mM and 182mM (corresponding to 200 and 500 SDS in a box of 16.6nm×16.6nm×16.6nm), are simulated at 300K. A 12µs simulation is performed for each concentration and the last 6µs of the trajectory is used for analysis. We calculate the cluster size distribution using a cutoff distance between C1 beads of 0.75nm24 within which the sulfate dodecyl molecules are taken to belong to the same cluster. We also study the adsorption of SDS surfactants on to hydrocarbon slabs at 323K. We firstly compare the potential of mean force (PMF) for pulling a single SDS surfactant from a hydrocarbon slab using the standard Martini FF and Dry Martini FF in simulation box of 5.8nm×5.8nm×15nm. In these simulations, the hydrocarbon slab is composed of 208 linear alkyl chains, each of 36 carbon atoms in length. In both standard Martini and Dry Martini FFs, the alkyl chain is represented by 9 C1 beads. The default parameters for the C1 beads in standard Martini are used. For the Dry Martini FF, the bonded and non-bonded parameters between C1 beads are shown in Table 1. We then study the adsorption of various concentrations of SDS surfactants onto the hydrocarbon slab with Dry Martini at 323K using the following procedures: 1) We firstly equilibrate a slab of hydrocarbons composed of linear alkyl chains of 36 carbon atoms at 323K. To do this, 400 chains are placed in a simulation box of dimensions 7.5nm×7.5nm×8nm. The ACS Paragon Plus 8 Environment

Page 8 of 31

Page 9 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

hydrocarbon slab is then equilibrated for 100ns and the average thickness of the hydrocarbon slab after reaching equilibrium is ~7.4nm. 2) We then place this equilibrated hydrocarbon slab in a longer box of dimensions 7.5nm×7.5nm×23.8nm. The top and bottom of the hydrocarbon slab are two surfaces that are separated by ~16.4nm in the periodic geometry. 3) Finally, between 50 and 400 SDS surfactants are randomly added to the gallery space between two surfaces so that SDS are homogenously distributed at the beginning of simulations. Correspondingly, the initial concentration of SDS surfactants in the gallery varies from 90mM to 717mM and the volume fraction of them ranges from 3.25vol% to 25vol%. We monitor the adsorption of SDS surfactants onto the hydrocarbon slab and the self-assembly in the gallery. In our simulations, a dodecyl sulfate anion is considered adsorbed if the distance between any C1 bead in the surfactant and the C1 bead in the hydrocarbon slab is less than a cutoff distance of 0.75nm. A 3µs simulation is performed for each concentration and the last 1µs of the trajectory is averaged to get the equilibrium amount of adsorption and to calculate the interfacial tension γ using pressure anisotropy by the following equation:39

1 2

γ = LZ PZ −

PX + PY 2

where ½ accounts for the presence of two interfaces; LZ is the simulation box length normal to the interface and PZ and (PX+PY)/2 are the normal and lateral components of the pressure tensor. The potential of mean force (PMF) is calculated using umbrella sampling40 in order to determine free energy profiles for removal of a surfactant either from a micelle or from a hydrophobic surface. Specifically, we generate a series of configurations by pulling one single target SDS surfactant away from the SDS micelle or hydrocarbon surface with a pulling force constant of 1000kJ/mol/nm2 and a pulling rate of 1nm/ns. We then extract 60-80 frames from the configurations that correspond to a center of mass spacing of 0.6Å between the target SDS

ACS Paragon Plus 9 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 31

surfactant and the SDS micelle or the hydrocarbon slab. Each configuration is pre-equilibrated for 2ns, followed by restrained molecular dynamics simulations for 150ns in an NPT ensemble for Martini and an NVT ensemble for Dry Martini. The restrained molecular dynamics simulations are performed using a harmonic position restraint with a force constant of 1000kJ/mol/nm2 applied to the center of mass of the target surfactant molecule. Note for the SDS micelle simulations, the position restraint is applied radially in the x, y, and z directions, while for the hydrocarbon surface simulations, the position restraint is only applied in the z direction. The PMF is then obtained using g_wham, i.e., the weighted histogram analysis method (WHAM) implementation in GROMCS.41 No correction for the geometric effect is performed for PMF calculations.42 It has been shown that compared to atomistic simulations the standard Martini FF results in a speed up of 3-4 orders of magnitude due to the use of large time step, reduced number of atoms, increased dynamics and the consideration of only short-range interactions.30 The Dry Martini FF is expected to further increase the simulation speed by removing waters. Besides, it also dramatically speeds the movement of molecules in the implicit water environment. Although the computational speed enhancement of Dry Martini over standard Martini depends greatly on the system size and the number of waters that can be removed, we find that for a system containing 60 SDS in a box size of 9nm×9nm×9nm at 300K, Dry Martini after removing ~6000 water beads is more than 10 times faster in nominal simulation time than standard Martini if the shifted cutoff is used in both FFs, not counting additional speed-up due to the faster transport that occurs in the absence of explicit water. The enhancement is expected to be much greater for larger systems or over Martini FFs with polarizable water beads.43 The implementation of PME in GROMACS is rather machine dependent because the allocation of processors to calculate the short-range (in real space) and long-range (in reciprocal space) interactions36 does not consider, for example, the underlying network or process clock

ACS Paragon Plus 10 Environment

Page 11 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

rates. For a specific system with 60 SDS in a box size of 9nm×9nm×9nm at 300K, Dry Martini with PME decreases the computational speed ~15 times compared to that with shifted cutoff. So Dry Martini with PME is actually comparable to that of standard Martini with shifted cut-off. Nevertheless, for the larger system we considered in this paper, i.e, 500 SDS in a box size of 16.6nm×16.6nm×16.6nm (182mM) at 300K, we could achieve around 2µs/day with 12 2.53 GHz Intel Xeon E5549 processors using Dry Martini with PME. Lastly, to gain optimal performance with PME in GROMACS, one might want to use g_tune_pme tool with short simulation runs to find the number of PME-only processors for the specific system.44 RESULTS AND DISCUSSIONS 1. Structural and Energetic Properties of Micelle Table 1. Structural properties of SDS micelles Rs (Å) Rg (Å) Dry Martini 15.52±0.31 19.82±0.27 Standard Martini 15.80±0.30 19.98±0.25 16 All-Atom (Palazzesi et al.) 15.0 19.55 United-Atom (Shang et al.)9 15.8 19.0 14 All-Atom (Bruce and Berkowitz) 16.2±0.12 19.6 11 All-Atom (Mackerell) 16.0±0.06 19.7±0.08 All-Atom (Rakitin and Pack)15 16.4 19.9

I1/I3 1.39+0.09 1.69±0.10 na 1.21 1.05 1.03±0.03 1.09

We simulate an SDS micelle composed of 60 surfactants in a box of 9nm×9nm×9nm at 300K using both Dry Martini and standard Martini. We calculate the micelle size measured by the average radius of gyration (Rg) and by the root-mean-squared distance between the head group and micelle center of mass (Rs) over a 100ns simulation. The results from Dry Martini are compared to those from standard Martini and from other atomistic force fields11,14–16 in Table 1. As shown, both Rg (15.52±0.31Å) and Rs (19.82±0.27Å) using Dry Martini are very close to those using standard Martini. The Rg value (15.52±0.31Å) is very close to the experimental value of 15.4Å from X-ray scattering45 and is at the lower end of results obtained from other simulations using

ACS Paragon Plus 11 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 12 of 31

atomistic FFs. On the other hand, Rs calculated from Dry Martini FF is fairly close to the result from atomistic simulations. We also characterize the shape of the micelle by its sphericity defined using the moments of inertia, I1, I2 and I3, in the order of decreasing magnitude, evaluated along each of the principal axes of the radius of gyration tensor. The ratios I1/I3, I2/I3, and I1/I2, averaged over 100ns, are 1.39:1.25:1.14. These values are in reasonable agreement with those obtained from the standard Martini FF, namely 1.69:1.10:1.10. Both FFs indicate that the micelle deviates somewhat from a spherical shape. We note that in Table 1 the ratio of maximum to minimum moments of inertia, I1/I3 calculated from Dry Martini is slightly smaller than that from standard Martini but is closer to other simulation studies.

Figure 2. Radial number density of various CG beads, including three C1 (alkane tail) beads as designated in Figure 1a, all C1 beads, and sulfate group, in a micelle composed of 60 SDS surfactants as a function of distance from micelle center of mass (COM). The radial density of hydrated sodium ions (Na) beads is multiplied by 5 for clarity. Lines represent calculations using standard Martini FF and open circles are for Dry Martini. The simulation is performed at 300K. We show in Figure 2 that the radial number densities of different CG beads calculated from Dry Martini agree well with those from standard Martini. For both FFs, the micellar interior is completely composed of hydrophobic tail and free of sulfate groups and sodium ions (as well as

ACS Paragon Plus 12 Environment

Page 13 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

waters in the Martini FF). Note here again that Dry Martini uses a dielectric constant of 150 so that the radial density of its hydrated sodium ion beads matches that from standard Martini. Using a smaller dielectric constant in Dry Martini FF, such as 80, steepens the radial density of hydrated sodium ions and thus worsens the comparison. One might wish to test the accuracy of CG results against atomistic predictions of ion density distributions, but unfortunately different atomistic FFs report different ion density profiles and different degrees of counterion binding.10

Figure 3. Potential of mean force (PMF) as a function of distance of center of mass (COM) of a “target” surfactant from the micelle COM at 300K from both regular and Dry Martini. Two representative snapshots are shown: on the left the target surfactant, with tail beads shown in red spheres, is inside the micelle and on the right it is outside. The tail beads of the other surfactants in the micelle are shown in cyan. Surfactant head groups are in yellow, and sodium ion beads in blue. The potential of mean force (PMF) for pulling one single SDS surfactant out of a micelle composed of 60 surfactants in a box 9nm×9nm×9nm at 300K is shown in Figure 3 as a function of distance from the micellar center of mass (COM). The PMF shows the change in potential energy when a SDS surfactant is pulled gradually away from a micelle. As is seen in Figure 3, PMF features a minimum at around 1.5nm from micellar COM, followed by a sharp increase to a maximum at around 3.1nm, and a gradual decrease at distances beyond 3.1nm. The maximum

ACS Paragon Plus 13 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 14 of 31

change in PMF is around 9.4kT for standard Martini and 9.7kT for Dry Martini, both of which are close to the PMF calculated using atomistic CHARMM FF (~9.0 kT).42 The slightly sharper increase of PMF in Dry Martini is probably related to the smaller micelle size as shown in Table 1 and a slightly thinner interfacial region between the hydrocarbon micellar core and water. The sharper PMF profile in Dry Martini is actually closer to the results using atomistic CHARMM and CG SDK FFs.42 Additionally, we observe that the decrease in the PMF at separation distances greater than 3.1nm is not as pronounced in Dry Martini as in Martini. The more gradual decrease of PMF in Martini occurs at the periphery of the micelle and is probably due to the weaker electrostatic interaction between charged sulfate groups and sodium ions. Even though we showed above that the density distribution of ions is almost identical for Dry Martini and standard Martini, the PMF profile is different in this region where water dominates. 2. Self-Assembly of SDS surfactants

Figure 4. Cluster size distribution for two different SDS surfactant concentrations, i.e. 72.8mM and 182mM, using Dry Martini. The temperature is 300K. To show that Dry Martini FF allows us to simulate the self-assembly process successfully with available computational resources, we perform 12µs simulations using Dry Martini at two different SDS surfactant concentrations: 72.8mM and 182mM (corresponding to 200 and 500 SDS

ACS Paragon Plus 14 Environment

Page 15 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

in a box of 16.6nm×16.6nm×16.6nm) at 300K. SDS surfactants are randomly placed in a simulation box at the beginning. Over time, small micelles form and they gradually merge to form big micelles. Using a cutoff distance of 0.75nm24 between C1 beads, we calculate the cluster size distribution using the last 6µs of the simulation trajectory. As shown in Figure 4, at 72.8mM, the range of micellar aggregation number is between 55 and 80 with a preferred micelle aggregation number of 67. The range of micelle sizes calculated from the simulation thus lies in the experimentally reported range of mean aggregation numbers between 45 and 77.46–48 There is also a certain fraction of free SDS surfactants in solutions, as evidenced by the peak at cluster size of 1. These single surfactants are in the process of being exchanged from one micelle to another. Additionally, if we label all the surfactants in an individual micelle at t=6000ns with a color specific to that micelle as shown in Figure 5a, we find that the micelles present at a later time of t=9000ns are well mixed with surfactants of different colors, indicating a thorough exchange of micellar content over 3µs time period. We emphasize here that 3µs nominal simulation time should correspond to more than 12µs of effective time because of the 4:1 mapping as well as the removal of water which speeds the diffusion of surfactants in implicit water. This confirms that using Dry Martini FF, we can achieve reasonable equilibration for SDS self-assembly. In fact, starting simulations with a different initial configuration where SDS surfactants are placed in a long cylindrical micelle, we obtain the same results. Furthermore, cylindrical micelles with larger cluster sizes form at a high surfactant concentration of 182mM. Besides the exchange of surfactants among different clusters, we also observe the formation and breakage of cylindrical micelles as seen in Figure 5b. The long cylindrical micelle observed at t=6000ns breaks at t=7000ns. The SDS surfactants in the original long cylindrical micelle at t=6000ns is distributed in several different micelles at t=7000ns. Additionally, a new long cylindrical micelle is seen at

ACS Paragon Plus 15 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 16 of 31

t=7000ns. At this high surfactant concentration, large clusters are part of the cluster size distribution. As seen in Figure 4, the cluster size distribution shifts right and covers a much larger range of cluster sizes from 50 to 100 surfactants. More importantly, as we discussed, very long cylindrical micelles also appear from time to time in our simulations as a result of micelle fusion. Thus small peaks around cluster sizes of 150 are also seen in the cluster size distribution. Further increasing surfactant concentration beyond 182mM results in more frequent formation of cylindrical micelles and an increase in their length. For example, at 364mM (i.e., 1000 SDS in a box of 16.6nm×16.6nm×16.6nm), a very long cylindrical micelle that spans the periodic box is seen. It forms, breaks and reforms over the course of simulation.

Figure 5. Simulation snapshots for SDS concentrations at a) 72.8mM and b) 182mM at 300K. In a) dodecyl sulfate anions in different micelles at t=6000ns are labeled with different colors and the colors of the individual surfactants are preserved until t=9000ns. In b), dodecyl sulfate anions in the largest micelle at t=6000ns are labeled in red and the color is preserved until t=7000ns. The tails of the other surfactants are in cyan, head groups in yellow, and sodium ions in blue.

ACS Paragon Plus 16 Environment

Page 17 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

3. Effect of PME on the Clustering of SDS Aggregates We find that the electrostatic interactions and how they are implemented, such as via a shifted cutoff or by a Particle-Mesh Ewald (PME) scheme, in either standard Martini or Dry Martini FF affects strongly the interactions among charged beads. When SDS surfactants assemble into aggregates, not only are there complicated intra-aggregate electrostatic interactions, but in addition the aggregates themselves are charged and thus are also affected by inter-aggregate electrostatic interactions. In addition, a recent study10 has shown that micellar structure is affected by the Lennard-Jones parameters for sodium ions, the ionic oxygens on sulfate group, and water models. Thus, depending on the atomistic FF, aggregates with 300 SDS surfactants or more could form either cylindrical micelles or unphysical lamellar bicelles. For standard Martini FF, we find that when the Coulombic interaction is shifted to zero from 0 to 1.2nm with a relative permittivity (εr) of 15 as suggested in the original publications,21,26 at low concentration, SDS surfactants form long cylindrical micelles instead of the expected spherical micelles, and when the concentration is high, SDS micelles (both spherical and cylindrical) tend to cluster. Using the PME electrostatics with εr=1522 also results in clustering of SDS micelles although the clustering is less obvious than it is with shifted cutoff electrostatics. For another similar surfactant with the same sulfate headgroup and a shorter alky tail, sodium octyl sulfate (SOS), clustering of spherical micelles are also observed using both shifted cutoff electrostatics and PME electrostatics with εr=15. However, using PME with εr=80, the clustering of SDS (and SOS) surfactant aggregates disappears in regular Martini. The radial densities of different CG beads using PEM with εr=80 are also comparable to that using shifted cutoff with εr=15 as seen in Figure S1. Additionally, as shown by another study,42 with shifted cutoff electrostatics and εr=15, the maximum change in the PMF for pulling a single SDS surfactant from its micelle is significantly higher (~10.5 kT) than that for CHARMM FF (~9.0 kT). The maximum

ACS Paragon Plus 17 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 18 of 31

change in the PMF using PME with εr=80 is 9.4 kT in Figure 3 and is actually much closer to that for the CHARMM FF than is the PMF for regular Martini with shifted cutoff instead of PME. Some previous studies using the standard Martini FF also indicated that PME electrostatics provides a more realistic description of the interactions between charged dendrimers and lipid membranes or between charged antimicrobial peptides and lipid membranes. 31,43,49–51

Figure 6. Snapshots from simulations using Dry Martini for different implementations of electrostatic interactions at 300K: a) shifted cutoff with a relative permittivity (εr) of 15 and the Coulombic interaction shifted to zero from 0 to 1.2nm for a concentration of 182mM SDS and no added salt; b) Particle-Mesh Ewald (PME) electrostatics with εr=80 and c) the same as b) but with εr=150, for a concentration of 182mM SDS surfactant and 72.8mM sodium chloride (NaCl) in both b) and c). The tail beads of the SDS are shown in cyan and the head groups in yellow, sodium ions in blue and chloride ions in red. For Dry Martini, we performed many simulations with various cutoff schemes for the electrostatic interaction, including cutoff, shifted cutoff, encad-shifted cutoff, and reaction fields, all with various cutoff distances ranging from 1.2nm to 2.0nm with εr=15 at 300K in a box size of 16.6nm×16.6nm×16.6nm. We find that different cutoff schemes always lead to clustering of SDS aggregates at high surfactant concentrations (for example, at 364 mM SDS concentration, i.e., 1000 SDS surfactants). An example for the shifted cutoff scheme with εr=15 is shown in Figure 6a

ACS Paragon Plus 18 Environment

Page 19 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

for an SDS solution at a lower concentration of 182mM (500 SDS). We also ran simulations with PME and varying value of εr in Dry Martini as seen in the supporting information. With εr at 80 or 120, at a surfactant concentration below 100mM, the clustering of SDS aggregates is not obvious, but the clustering becomes obvious at higher SDS concentration (for example at 364mM of SDS surfactants, i.e., 1000 SDS) or with added salt (for example, in Figure 6b with 182mM of SDS surfactant and 72mM of sodium chloride (NaCl), i.e. 500 SDS and 200 NaCl). Interestingly, we find that a further increase of εr to 150 with PME eliminates this clustering. At a high SDS concentration of 364mM, we observe formation of long cylindrical micelles containing more than 200 surfactants, which is much longer than the cylinder that appears at 182mM. In fact the size of the cylindrical micelles can sometimes reach 800 surfactants and still no clustering is seen at this concentration. Additionally, as seen in Figure 6c, adding 72mM sodium chlorides to a system containing only 182mM SDS surfactants, we also observe the formation of long cylindrical micelles and most importantly no clustering occurs. Compared to the size distribution at 182mM in Figure 4, adding salt induces more cylindrical micelles and their sizes can be as large as 450 surfactants. The sphere-to-rod transition due to increasing surfactant concentration or salt concentration has also been seen previously.1,2,5,6,52 Lastly, we have found through a number of simulations that the same scheme of using PME and modifying the value of relative permittivity does

not

prevent

aggregation

of

micelles

of

zwitterionic

surfactants,

such

as

dodecylphosphocholine (DPC), in Dry Martini and further modifications such as changing the charge bead interactions might be necessary in order to obtain reasonable results. 4. Adsorption of SDS Surfactants on Hydrocarbon Slab Surfactants can also be adsorbed onto colloidal surfaces or oil droplets in aqueous solutions. Such important adsorption processes are found in many applications including lubrication, detergency, oil recovery, wetting and emulsification. Compared to the self-assembly

ACS Paragon Plus 19 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 20 of 31

of surfactants in solutions, adsorption is more complicated due to the presence of surfaces and usually involves multiple steps, interesting structural transitions, and different thermodynamics and kinetics.53,54 For example, not only do surfactants self-assemble and form aggregates in solution, but there is also a competition between these aggregates and the surfaces for free surfactants. Simulations of self-assembly of surfactants at surfaces are also susceptible to equilibration problems similar to those encountered in self-assembly of micelles in solution.

Figure 7. Potential of mean force (PMF) as a function of distance from the surface of a hydrocarbon slab composed of linear alkyl chains of 36 carbon atoms each. Two representative snapshots from the simulations are shown: one with the target surfactant tail inside, and the other with it outside, the hydrocarbon slab. The tail beads of the surfactant are shown in red, head beads in yellow, sodium ions in blue, and the alkyl chains in pink. The temperature is 323K. We firstly compare the PMF for pulling a single SDS surfactant from a hydrocarbon slab at 323K calculated using standard Martini and Dry Martini FF. In these simulations, the hydrocarbon slab is composed of 208 linear alkyl chains each of 36 carbon atoms in length. The simulation box is 5.8nm×5.8nm×15nm. We plot the PMF as a function of distance between the center of mass (COM) of the SDS surfactant and the hydrocarbon surface in Figure 7. The location of the hydrocarbon surface is defined as the inflection point of the hydrocarbon density at the

ACS Paragon Plus 20 Environment

Page 21 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

hydrocarbon/water interface and is set at 0 in Figure 7. As is seen, the PMF exhibits a minimum at the hydrocarbon surface and increases quickly to reach a plateau when the separation distance exceeds 1.4nm, at which point the SDS molecule has completely broken contact with the hydrocarbon surface. We find that PMF profile calculated from Dry Martini agrees very well with that from standard Martini and we only observe a small difference in the total potential energy difference, i.e., 17kT for Dry Martini vs. 17.7kT for standard Martini. We further note that the PMF for pulling an SDS from the hydrocarbon slab (even at 323K) is much larger than for pulling the SDS surfactant from a micelle at 300K.

Figure 8. a) Equilibrium surface coverage as a function of SDS surfactant concentration. Shown in the inset is a simulation snapshot at one specific surfactant concentration. The SDS tail beads are shown in cyan, head group in yellow, sodium ions in blue, and the alkyl chains in pink. b) Interfacial tension as a function of equilibrium surface coverage. The temperature is 323K. We then simulate the adsorption of SDS surfactants onto a slab of hydrocarbons composed of 400 linear alkyl chains, each of a length of 36 carbon atoms, in a 7.5nm×7.5nm×23.8nm box. The temperature is set to an elevated value of 323K to speed the simulations. We found that some sulfate dodecyl anions are adsorbed onto the surface almost immediately and the rest of surfactants in the gallery space form aggregates. As discussed in Figure S5 in the supporting information, the surfactant in the gallery are adsorbed to the surface via two different pathways: 1) they are ACS Paragon Plus 21 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 22 of 31

adsorbed to the hydrocarbon surface as a micelle and correspondingly a sharp increase in the number of adsorbed surfactants at a particular time is seen; or 2) some surfactants escape their micelles formed in the gallery and are then adsorbed onto the hydrocarbon surfaces. In the latter case, the number of adsorbed surfactants increases gradually with time. Please note that some adsorbed surfactants could still leave the respective surface and become free in the solutions for a short period of time before they are reabsorbed by the surfaces again. When a surface is fully saturated with surfactants at high concentrations, the non-adsorbed surfactants form spherical micelles or cylindrical micelles in the gallery as seen in Figure S4. Similar to what we discuss in previous sections, the micelles exist dynamically with micelle fission and fusion occurring. When the hydrocarbon slab is present, the surfactants can also exchange among surfactant aggregates in the gallery space and the surfaces. Lastly, we note that once the surfactants are adsorbed on surfaces, the hydrophobic alkyl tails are mostly in contact with the hydrocarbon slab and the sulfate groups are distributed uniformly along the surfaces mainly due to the electrostatic interactions between them, and no clustering of SDS surfactants on surfaces is observed. In Figure 8a, we plot the equilibrium surface coverage as a function of SDS concentration, both total and bulk (or non-adsorbed) surfactant concentration. The equilibrium surface coverage is calculated by dividing the equilibrium number of adsorbed surfactants by the surface area of the hydrocarbon slab, which includes both surfaces of the slab. The total surfactant concentration refers to the initial concentration in our simulation and is calculated by dividing the total number of surfactants by the volume of the gallery space, while the bulk concentration is obtained by dividing the average number of surfactants not adsorbed at equilibrium by the volume of gallery space. As is seen, with an increase in total surfactant concentration, surface coverage increases linearly before reaching a peak surface coverage around 1.68 sulfate dodecyl per nm2. In the region of low surfactant concentration, nearly all surfactants are absorbed on the hydrocarbon slab

ACS Paragon Plus 22 Environment

Page 23 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

surfaces, leaving almost no surfactants in the gallery (although single surfactants do occasionally leave the surfaces and enter the gallery before they are reabsorbed by a surface again). In another words, surface coverage increases without a detectable change in the bulk concentration, as seen in Figure 8a. Interestingly, we find that a further increase in surfactant concentration leads to a slight drop in the surface coverage. This is because at the concentration near the drop, the surface is already crowded with surfactants, and a slight increase in the surfactant concentration makes possible the formation of a single micelle in the gallery, as seen in Figure 8a inset. The micelle is only favored if its aggregation number is above a critical value. Hence, if the concentration of surfactant is high enough to saturate the surface with only enough left over to form a micelle fragment, rather than a whole micelle, the micelle fragment will scavenge some surfactants from the surface in order to form a complete micelle. This causes a dip in surface concentration over this limited range of total concentration. Although some previous experimental studies also reported a decrease in the SDS adsorption at high concentration and the occurrence of a local maximum in adsorbed concentration instead of a plateau,55,56 the experimental maximum is presumably unrelated to the maximum we observed. In our simulations, the maximum is an artifact of the small size of our simulation box, which contains only enough surfactant to form a single micelle near the concentration where the plateau appears. In our simulation, after the dip in surface coverage, as the surfactant concentration increases further, it increases again and reaches its saturation level around 1.70 sulfate dodecyl per nm2. Above this level, a further increase in surfactant concentration only results in more surfactants in bulk, as evidenced by the increase in bulk concentration. These surfactants in the bulk form spherical micelles or cylindrical micelles at high concentration. Lastly, the increasing surface coverage is accompanied by a decrease in interfacial tension, as seen in Figure 8b. In the absence of surfactants, the interfacial tension calculated in Dry Martini FF is 41mN/m at 50C. Even though the interfacial tension exhibits a

ACS Paragon Plus 23 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 24 of 31

complex dependence on initial surfactant concentration as shown in Figure S6, all points fall on one simple curve in Figure 8b as a function of surface coverage, i.e., the interfacial tension decreases in a concave down fashion when the surfaces are covered by more and more dodecyl sulfate anions. It reaches a minimum value around 5mN/m at the highest surfactant concentration. CONCLUSIONS We have parameterized an ionic surfactant, sodium dodecyl sulfate (SDS), within the coarse-grained (CG) Dry Martini force field (FF) with implicit water so that the structural properties of a micelle composed of 60 surfactants obtained using Dry Martini are in good agreement with those using standard Martini with explicit water. The results are also similar to many other atomistic simulations and experimental studies. We found that replacing the shifted cutoff electrostatic interaction with long-range particle-mesh Ewald (PME) electrostatics with a high dielectric constant (εr=150), was necessary to prevent micelles from unphysically aggregating. The potential of mean force (PMF) calculated for Dry Martini by pulling a single surfactant from a micelle exhibits a similar free energy profile to that obtained using standard Martini especially at small distance from micelle center of mass when the surfactant is inside the micelle, and the maximum free energy change on extraction of the surfactant is around 9.7kT for Dry Martini versus 9.4kT for standard Martini with PME. When the surfactant is outside the micelle, a slower decrease in the PMF beyond the local maximum is obtained for Dry Martini than for standard Martini, possibly due to the neglect of explicit water in the former. By removing the explicit water molecules, the CG Dry Martini FF dramatically reduces the number of particles in the simulation and quickens molecular dynamics. As a result, it becomes possible to obtain equilibrium cluster size distributions for surfactant self-assembly in solutions. At low surfactant concentration, spherical micelles with preferable size around 67 are attained. Increasing the surfactant concentration leads to merger of spherical micelles and dynamic formation of long

ACS Paragon Plus 24 Environment

Page 25 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

cylindrical micelles that grow, break and re-form with time. Thus the distribution of micelle sizes extends to larger aggregation numbers. Similarly, addition of salts also induces the transition of spherical micelles to cylindrical micelles, in qualitative agreement with experiments. We also studied the adsorption of SDS surfactants onto hydrocarbon slab surfaces. The PMF for pulling a single surfactant from the hydrocarbon slab using Dry Martini shows excellent agreement with that using standard Martini, with a maximum free energy change of 17kT for Dry Martini vs. 17.7kT for standard Martini, about 75% higher than the maximum PMF for extraction of an SDS surfactant from an SDS micelle. The adsorption of surfactants to the hydrocarbon surface is studied as a function of the surfactant concentration. Two different absorption processes are observed: 1) surfactants form aggregates in the gallery and then are adsorbed to surfaces as an aggregate; 2) surfactants can also escape as individual molecules from the aggregates in the gallery and be adsorbed to surfaces. The former leads to a jump in surface adsorption as a function of time, while the latter results in a gradual increase in the number of adsorbed surfactants. When more and more surfactants are added in the gallery, the surface coverage increases almost linearly with surfactant concentration before reaching a peak value around 1.68 surfactant per nm2, below which nearly all the surfactants are adsorbed on surfaces. However, the surfaces are crowded with surfactants so further increase in surfactant concentration leads to an increase in bulk surfactant concentration and formation of aggregates in the gallery. Interestingly there is a local dip in the surface coverage when the first micelle forms as overall concentration increases. Lastly, we showed that the interfacial tension decreases from 41mN/m to 5mN/m with an increase in surface coverage following a concave down fashion. SUPPORTING INFORMATION Comparison of micelle structures using different force fields; Effect of electrostatics on SDS selfassembly; Snapshots from adsorption simulations for three different SDS concentrations; Number

ACS Paragon Plus 25 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 26 of 31

of adsorbed surfactants with time; Interfacial tension as a function of initial surfactant concentration. This information is available free of charge via the Internet at http://pubs.acs.org. AUTHOR INFORMATION *E-mail: [email protected] ACKNOWLEDGEMENTS We are grateful to Dow Chemical Company for funding this research, to Valeriy Ginzburg and Antony van Dyke of Dow Chemical Company for valuable discussions and inspiration that guided the research, and to Siewert-Jan Marrink for helpful comments. REFERENCES

(1)

Berr, S.; Jones, R. Effect of Added Sodium and Lithium Chlorides on Intermicellar Interactions and Micellar Size of Aqueous Dodecyl Sulfate Aggregates as Determined by Small-Angle Neutron Scattering. Langmuir 1988, 4, 1247–1251.

(2)

Bezzobotnov, V.; Borbely, S. Temperature and Concentration Dependence of Properties of Sodium Dodecyl Sulfate Micelles Determined from Small Angle Neutron Scattering Experiments. J. Phys. Chem. 1988, 92, 5738–5743.

(3)

Reiss-Husson, F.; Luzzati, V. The Structure of the Micellar Solutions of Some Amphiphilic Compounds in Pure Water as Determined by Absolute Small-Angle X-Ray Scattering Techniques. J. Phys. Chem. 1964, 844, 3504–3511.

(4)

Kodama, M.; Kubota, Y.; Miura, M. The Second CMC of the Aqueous Solution of Sodium Dodecyl Sulfate. III. Light-Scattering. Bull. Chem. Soc. Jpn. 1972, 45, 2953–2955.

(5)

Mazer, N.; Benedek, G.; Carey, M. An Investigation of the Micellar Phase of Sodium Dodecyl Sulfate in Aqueous Sodium Chloride Solutions Using Quasielastic Light Scattering Spectroscopy. J. Phys. Chem. 1976, 80, 1075–1085.

(6)

Missel, P.; Mazer, N.; Benedek, G.; Young, C.; Carey, M. Thermodynamic Analysis of the Growth of Sodium Dodecyl Sulfate Micelles. J. Phys. Chem. 1980, 84, 1044–1057.

(7)

Missel, P.; Mazer, N.; Carey, M.; Benedek, G. Influence of Alkali-Metal Counterion Identity on the Sphere-to-Rod Transition in Alkyl Sulfate Micelles. J. Phys. Chem. 1989, 93, 8354–8366.

(8)

Ma, C.; Li, G.; Xu, Y.; Wang, H.; Ye, X. Determination of the First and Second CMCs of Surfactants by Adsorptive Voltammetry. Colloids Surf., A 1998, 143, 89–94.

ACS Paragon Plus 26 Environment

Page 27 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(9)

Shang, B. Z.; Wang, Z.; Larson, R. G. Molecular Dynamics Simulation of Interactions between a Sodium Dodecyl Sulfate Micelle and a Poly(Ethylene Oxide) Polymer. J. Phys. Chem. B 2008, 112, 2888–2900.

(10)

Tang, X.; Koenig, P. H.; Larson, R. G. Molecular Dynamics Simulations of Sodium Dodecyl Sulfate Micelles in Water-the Effect of the Force Field. J. Phys. Chem. B 2014, 118, 3864–3880.

(11)

MacKerell Jr, A. Molecular Dynamics Simulation Analysis of a Sodium Dodecyl Sulfate Micelle in Aqueous Solution: Decreased Fluidity of the Micelle Hydrocarbon Interior. J. Phys. Chem. 1995, 99, 1846–1855.

(12)

Yoshii, N.; Iwahashi, K.; Okazaki, S. A Molecular Dynamics Study of Free Energy of Micelle Formation for Sodium Dodecyl Sulfate in Water and Its Size Distribution. J. Chem. Phys. 2006, 124, 184901.

(13)

Dominguez, H.; Berkowitz, M. L. Molecular Dynamics Simulation of Pyrene Solubilized in a Sodium Dodecyl Sulfate Micelle. J. Phys. Chem. B 2000, 104, 5302–5308.

(14)

Bruce, C.; Berkowitz, M. Molecular Dynamics Simulation of Sodium Dodecyl Sulfate Micelle in Water: Micellar Structural Characteristics and Counterion Distribution. J. Phys. Chem. B 2002, 106, 3788–3793.

(15)

Rakitin, A. R.; Pack, G. R. Molecular Dynamics Simulations of Ionic Interactions with Dodecyl Sulfate Micelles. J. Phys. Chem. B 2004, 108, 2712–2716.

(16)

Palazzesi, F.; Calvaresi, M.; Zerbetto, F. A Molecular Dynamics Investigation of Structure and Dynamics of SDS and SDBS Micelles. Soft Matter 2011, 7, 9148–9156.

(17)

Bruce, C.; Senapati, S.; Berkowitz, M.; Perera, L.; Forbes, M. Molecular Dynamics Simulations of Sodium Dodecyl Sulfate Micelle in Water: The Behavior of Water. J. Phys. Chem. B 2002, 106, 10902–10907.

(18)

Sammalkorpi, M.; Karttunen, M.; Haataja, M. Structural Properties of Ionic Detergent Aggregates: A Large-Scale Molecular Dynamics Study of Sodium Dodecyl Sulfate. J. Phys. Chem. B 2007, 111, 11722–11733.

(19)

Gao, J.; Ge, W.; Hu, G.; Li, J. From Homogeneous Dispersion to Micelles a Molecular Dynamics Simulation on the Compromise of the Hydrophilic and Hydrophobic Effects of Sodium Dodecyl Sulfate. Langmuir 2005, 21, 5223–5229.

(20)

Sanders, S. a; Panagiotopoulos, A. Z. Micellization Behavior of Coarse Grained Surfactant Models. J. Chem. Phys. 2010, 132, 114902.

(21)

Jalili, S.; Akhavan, M. A Coarse-Grained Molecular Dynamics Simulation of a Sodium Dodecyl Sulfate Micelle in Aqueous Solution. Colloids Surf., A 2009, 352, 99–102.

ACS Paragon Plus 27 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 28 of 31

(22)

Pires, J.; Moura, A.; Freitas, L. Investigating the Spontaneous Formation of SDS Micelle in Aqueous Solution Using a Coarse-Grained Force Field. Quím. Nova. 2012, 35, 978–981.

(23)

Shinoda, W.; DeVane, R.; Klein, M. L. Coarse-Grained Force Field for Ionic Surfactants. Soft Matter 2011, 7, 6178.

(24)

LeBard, D. N.; Levine, B. G.; Mertmann, P.; Barr, S. A.; Jusufi, A.; Sanders, S.; Klein, M. L.; Panagiotopoulos, A. Z. Self-Assembly of Coarse-Grained Ionic Surfactants Accelerated by Graphics Processing Units. Soft Matter 2012, 8, 2385–2397.

(25)

Shinoda, W.; Devane, R.; Klein, M. Multi-Property Fitting and Parameterization of a Coarse Grained Model for Aqueous Surfactants. Mol. Sim. 2007, 33, 27–36.

(26)

Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. The MARTINI Force Field: Coarse Grained Model for Biomolecular Simulations. J. Phys. Chem. B 2007, 111, 7812–7824.

(27)

Jusufi, A.; Hynninen, A.-P.; Panagiotopoulos, A. Z. Implicit Solvent Models for Micellization of Ionic Surfactants. J. Phys. Chem. B 2008, 112, 13783–13792.

(28)

Arnarez, C.; Uusitalo, J.; Masman, M.; Ingólfsson, H. I.; de Jong, D.; Melo, M.; Periole, X.; de Vries, A.; Marrink, S. Dry Martini, a Coarse Grained Implicit Water Force Field for (bio)molecular Simulations. J. Chem. Theory Comput. 2015, DOI: 10.1021/ct500477k.

(29)

Mirzoev, A.; Lyubartsev, A. P. Systematic Implicit Solvent Coarse Graining of Dimyristoylphosphatidylcholine Lipids. J. Comput. Chem. 2014, 35, 1208–1218.

(30)

Marrink, S.; Vries, A. de; Mark, A. Coarse Grained Model for Semiquantitative Lipid Simulations. J. Phys. Chem. B 2004, 750–760.

(31)

Lo, C. A.; Rzepiela, A. J.; de Vries, A. H.; Dijkhuizen, L.; Hu, P. H.; Marrink, S. J. Martini Coarse-Grained Force Field : Extension to Carbohydrates. J. Chem. Theory Comput. 2009, 5, 3195–3210.

(32)

Monticelli, L.; Kandasamy, S. The MARTINI Coarse-Grained Force Field: Extension to Proteins. J. Chem. Theory Comput. 2008, 4, 819–834.

(33)

Lee, H.; Pastor, R. W. Coarse-Grained Model for PEGylated Lipids: Effect of PEGylation on the Size and Shape of Self-Assembled Structures. J. Phys. Chem. B 2011, 115, 7830– 7837.

(34)

Rossi, G.; Monticelli, L.; Puisto, S. R.; Vattulainen, I.; Ala-Nissila, T. Coarse-Graining Polymers with the MARTINI Force-Field: Polystyrene as a Benchmark Case. Soft Matter 2011, 7, 698–708.

ACS Paragon Plus 28 Environment

Page 29 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

(35)

Lee, H.; de Vries, A. H.; Marrink, S. J.; Pastor, R. W. A Coarse-Grained Model for Polyethylene Oxide and Polyethylene Glycol: Conformation and Hydrodynamics. J. Phys. Chem. B 2009, 113, 13186–13194.

(36)

Hess, B.; Kutzner, C.; van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435–447.

(37)

Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular Dynamics with Coupling to an External Bath. J. Chem. Phys. 1984, 81, 3684.

(38)

Essmann, U.; Perera, L. A Smooth Particle Mesh Ewald Method. J. Chem. Phys. 1995, 103, 8577–8593.

(39)

Kirkwood, J. G.; Buff, F. P. The Statistical Mechanical Theory of Surface Tension. J. Chem. Phys. 1949, 17, 338–343.

(40)

Torrie, G.; Valleau, J. Nonphysical Sampling Distributions in Monte Carlo Free-Energy Estimation: Umbrella Sampling. J. Comput. Phys. 1977, 23.

(41)

Hub, J. S.; de Groot, B. L.; van der Spoel, D. g_wham - A Free Weighted Histogram Analysis Implementation Including Robust Error and Autocorrelation Estimates. J. Chem. Theory Comput. 2010, 6, 3713–3720.

(42)

Yuan, F.; Wang, S.; Larson, R. G. Potentials of Mean Force and Escape Times of Surfactants from Micelles Using MD Simulations. Langmuir 2015, accepted.

(43)

Yesylevskyy, S. O.; Schäfer, L. V; Sengupta, D.; Marrink, S. J. Polarizable Water Model for the Coarse-Grained MARTINI Force Field. PLoS Comput. Biol. 2010, 6, e1000810.

(44)

Kutzner, C.; Apostolov, R.; Hess, B.; Grubmüller, H. Scaling of the GROMACS 4.6 Molecular Dynamics Code on SuperMUC. In Parallel Computing: Accelerating Computational Science and Engineering (CSE); Bader, M.; Bode, A.; Bungartz, H.-J.; Gerndt, M.; Joubert, G. R.; Peters, F., Eds.; ISO Press BV, 2014; pp. 722–730.

(45)

Itri, R.; Amaral, L. Distance Distribution Function of Sodium Dodecyl Sulfate Micelles by X-Ray Scattering. J. Phys. Chem. 1991, 95, 423–427.

(46)

Hayter, J. B.; Penfold, J. Determination of Micelle Structure and Charge by Neutron SmallAngle Scattering. Colloid. Polym. Sci. 1983, 261, 1022–1030.

(47)

Lianos, P.; Zana, R. Fluorescence Probe Studies of the Effect of Concentration on the State of Aggregation of Surfactants in Aqueous Solution. J. Colloid Interface Sci. 1981, 84, 100– 107.

ACS Paragon Plus 29 Environment

Langmuir

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 30 of 31

(48)

Bales, B. L.; Messina, L.; Vidal, A.; Peric, M.; Nascimento, O. R. Precision Relative Aggregation Number Determinations of SDS Micelles Using a Spin Probe. A Model of Micelle Surface Hydration. J. Phys. Chem. B 1998, 102, 10347–10358.

(49)

Lee, H.; Larson, R. G. Coarse-Grained Molecular Dynamics Studies of the Concentration and Size Dependence of Fifth- and Seventh-Generation PAMAM Dendrimers on Pore Formation in DMPC Bilayer. J. Phys. Chem. B 2008, 112, 7778–7784.

(50)

Lee, H.; Larson, R. G. Molecular Dynamics Simulations of PAMAM Dendrimer-Induced Pore Formation in DPPC Bilayers with a Coarse-Grained Model. J. Phys. Chem. B 2006, 110, 18204–18211.

(51)

Rzepiela, A. J.; Sengupta, D.; Goga, N.; Marrink, S. J. Membrane Poration by Antimicrobial Peptides Combining Atomistic and Coarse-Grained Descriptions. Faraday Discuss. 2010, 144, 431–443.

(52)

Velinova, M.; Sengupta, D.; Tadjer, A. V; Marrink, S.-J. Sphere-to-Rod Transitions of Nonionic Surfactant Micelles in Aqueous Solution Modeled by Molecular Dynamics Simulations. Langmuir 2011, 27, 14071–14077.

(53)

Somasundaran, P.; Huang, L. Adsorption/aggregation of Surfactants and Their Mixtures at Solid-Liquid Interfaces. Adv. Colloid Interface Sci. 2000, 88, 179–208.

(54)

Zhang, R.; Somasundaran, P. Advances in Adsorption of Surfactants and Their Mixtures at Solid/solution Interfaces. Adv. Colloid Interface Sci. 2006, 123-126, 213–229.

(55)

Arnebrant, T.; Bäckström, K. An Ellipsometry Study of Ionic Surfactant Adsorption on Chromium Surfaces. J. Colloid Interface Sci. 1989, 128, 303–312.

(56)

Vold, R.; Sivaramakrishnan, N. The Origin of the Maximum in the Adsorption Isotherms of Association Colloids. J. Phys. Chem. 1958, 62, 984–989.

ACS Paragon Plus 30 Environment

Page 31 of 31

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Langmuir

TOC Graphic

Coarse-Grained Molecular Dynamics Simulation of Self-Assembly and Surface Adsorption of Ionic Surfactants Using an Implicit Water Model Shihu Wang and Ronald G. Larson* Department of Chemical Engineering, University of Michigan 2300 Hayward Street, Ann Arbor, Michigan 48109-2136

ACS Paragon Plus 31 Environment