Molecular Dynamics Simulation of SDS and CTAB ... - ACS Publications

Article pubs.acs.org/Langmuir

Molecular Dynamics Simulation of SDS and CTAB Micellization and Prediction of Partition Equilibria with COSMOmic Sandra Storm,*,† Sven Jakobtorweihen,‡ Irina Smirnova,† and Athanassios Z. Panagiotopoulos§ †

Institute of Thermal Separation Processes, ‡Institute of Chemical Reaction Engineering, Hamburg University of Technology, D-21073 Hamburg, Germany § Department of Chemical and Biological Engineering, Princeton University, Princeton, New Jersey 08544, United States ABSTRACT: Molecular dynamics (MD) simulations of the self-assembly of different ionic surfactants have been performed in order to obtain representative micellar structures. Subsequently, these structures were used to predict the partition behavior of various solutes in these micelles with COSMOmic, an extension of COSMO-RS. This paper includes multiple self-assembled micelles of SDS (sodium dodecyl sulfate, anionic surfactant) and CTAB (cetyltrimethylammoniumbromide, cationic surfactant) at different concentrations. Micellar size, density profiles, and shape (eccentricity) have been investigated. However, the size strongly depends on the functional definition of a micelle. For this reason, we present a method based on the free monomer concentration in aqueous solution as an optimization criterion for the micelle definition. The combination of MD with COSMOmic has the benefit of combining detailed atomistic information from MD with fast calculations of COSMOmic. For the first time the influence of micelle structure on pratition equilibria, predicted with COSMOmic, were investigated. In case of SDS more than 4600 and for CTAB more than 800 single micelles have been studied. The predictions of the partition coefficients with COSMOmic are in good agreement with experimental data. Additionally, the most favorable locations of selected molecules in the micelles as well as probable energy barriers are determined even for complex solutes such as toluene, propanolol, ephedrine, acetone, phenol, lidocaine, syringic acid, coumarin, isovanillin, ferulic acid, and vanillic acid. This method can therefore be applied as a potential screening tool for solutes (e.g., drugs) to find the optimal solute−surfactant combination.

1. INTRODUCTION Amphiphilic molecules consist of polar groups (head groups) and nonpolar groups (tail groups). Depending on the character of the head group, amphiphiles can be distinguished between anionic, cationic, zwitterionic, and nonionic. In this paper, SDS and CTAB have been investigated. The main reason for selecting these two surfactants is due to their high relevance and the large basis of experimental data available.1 Additionally, one is an anionic (SDS) and the other a cationic (CTAB) surfactant, but the tail groups only differ in size and the core of the micelles is comparable. The contribution of the different head groups on the partition behavior of different solute molecules was investigated. Micelles formed by these surfactants can be used for extraction processes, because different components (solutes) can be solubilized depending on their chemical structure. On the one hand, micelles (or liposomes) can be used to mimic the cell membrane in the human body to study the drug uptake,2 and on the other hand they can serve as solubility enhancers in the food industry to solubilize different food additives like vitamins.3 The understanding of the self-assembly process of surfactants to form micelles is of great interest not only for a fundamental understanding of this process, but also for technical © 2013 American Chemical Society

applications. Since the self-assembly process operates in microscopic scale, it is difficult to study via experiments. Experimental techniques, which are mostly applied to study the micellization, are, for example, light scattering or neutron scattering and are time-consuming and expensive. However, computer simulations provide an alternative way to calculate micelle properties and to understand the process of micelle formation. There are several different simulation methods and models available, with different typical applications for each. The most detailed modeling approach, and thus computationally most demanding, is based on all-atom molecular dynamics (MD) simulations. Many groups, who have studied micellar systems with MD, have chosen preassembled structures to shorten the simulation time.4−6 Consequently, in these studies, only one micelle size is investigated and they are based on the assumption that the preassembled structures are close to the global energy minimum. When starting the MD simulations from a random configuration of surfactant molecules in water, Received: July 2, 2013 Revised: August 13, 2013 Published: August 13, 2013 11582

dx.doi.org/10.1021/la402415b | Langmuir 2013, 29, 11582−11592

Langmuir

Article

which is afterward filled with water to match a predetermined concentration by using the tool genbox of GROMACS. An insertion of a water molecule is successful, if the van der Waals criterion is fulfilled, otherwise, the solvent molecule will be removed. The energies of the initial configurations were minimized with the steepest descent method. After the initialization, a short simulation in the NVT ensemble has been performed for 600 ps, using the Nosé-Hoover thermostat28 with the coupling time constant of 1 ps and T = 298 K. After these initialization steps, all simulations were performed in the NPT ensemble at p = 1 bar and T = 298 K, using the Nosé-Hoover thermostat and the Parrinello-Rahman barostat29,30 (coupling constant τp = 2 ps). The temperatures of the surfactants, counterions, and solvent were controlled independently with the coupling constant of 1 ps. For the NPT simulations, a time step of 0.002 ps has been applied. The Lennard−Jones interactions were cut off at 1.2 nm and switched from 0.8 nm on.31,32 For the Coulomb interactions, a real space cutoff of 1.0 nm was used, and the particle mesh Ewald method33,34 for longranged electrostatic interactions was chosen. All bonds to hydrogen atoms were constrained using the LINCS algorithm.35 Recently, Piggot el al.36 compared different force fields for lipids. Although their recommendations for CHARMM36 force field (CHARMM TIP3P water model, 0.8 nm Lennard−Jones switching) are based on bilayer simulations, they were followed here. 2.2. The Model COSMOmic. In the COSMOmic model, which was developed by Klamt et al.14 as an extension of COSMO-RS15−20 a micelle is considered as a layered liquid. Each layer has a specific composition with respect to the atoms of the surfactants, ions, and water molecules. The composition within a layer is assumed to be homogeneous. Therefore, the atomic distribution in the micelle/water system has to be known either from experimental data or from simulations. The program COSMOtherm (including the COSMOmic package, version C30_130137) can so far distinguish between spherical, cylindrical, and laminar micelles. Besides the composition matrix, COSMOmic also requires DFT/COSMO files for all components of the system and all solutes of interest. These DFT/ COSMO files were calculated with Turbomole 5.1038 on the BPTZVP 39−41 level and the RI (resolution of the identity) approximation.42 The DFT/COSMO files include the histogram of the surface area with respect to the screening charge density, which is called σ-profile in the framework of COSMO-RS. A conformer analysis with HyperChem 8.043 of the solute molecules has been performed, if necessary. The explicit consideration of conformers (optimized in their correct environment) in COSMO-RS has been proven to give better predictions when compared to experiments.44 A representative surfactant monomer has to be selected for the assignment of the surface charges to the composition matrix of the micelle/water system. In our previous work,45 we demonstrated to use the lipid conformer, which has a solvent accessible surface (SAS) close to the average SAS within the considered vesicle. For the sodium ion, which acts as a counterion for SDS, the optimized radius for COSMO-RS was taken from ref 46. Molecular interactions are expressed in the model COSMOmic as pairwise surface interactions based on the polarization densities σ and σ′ of the surfaces (likewise COSMO-RS).14 Statistical thermodynamics is then applied to calculate free energies and chemical potentials according to the interaction of the surface segments. The resulting chemical potential for a solute X at a fixed position and orientation in a micellar system M can be written as:14

this assumption is not necessarily correct, as different sizes and structures of micelles are formed within the simulation.7−10 However, one drawback of the start from a random configuration is the difficulty of reaching equilibrium conditions. As shown by Sammalkorpi et al.10 even large-scale simulations with a time scale of about 200 ns are not sufficient to reach the equilibrium for SDS at the concentration of 1 mol/ L at ambient conditions. There is also good evidence, that even microsecond simulations are not sufficient.11 Simulations of even shorter time have been performed so far for CTAB. The surfactant CTAB has been studied so far starting from different sizes of preassembled micelles (aggregation numbers Nagg: 50, 60, 70, 80, 90, 100, 110, and 120) by Catá et al.12 and by selfassembly by Wang et al.13 starting from 151 CTAB molecules in a simulation box filled with 32 400 water molecules (cCTAB = 0.22 mol/L). In this work, the self-assembly of SDS and CTAB has been studied starting from a nonaggregated configuration in water at different concentrations to retrieve different sizes and shapes of micelles, which are realistic in comparison to experiments. The simulations have been performed for 100 ns. For higher surfactant concentrations, also multiple independent simulations have been accomplished in order to achieve better statistics. The micelles obtained by MD simulations have been analyzed concerning their structure, which is of special importance for the predictions with COSMOmic.14 The model COSMOmic is an extension of COSMO-RS,15−20 the conductor-like screening model for realistic solvation. For COSMOmic, the three-dimensional structure of a micelle from MD simulations is used as an input. Thereby, the free energy profile along the micelle of a wide variety of solutes can be predicted as well as partition coefficients. Recently, we have compared COSMOmic with COSMO-RS and the “pseudophase” approach,2 which relies on the assumption of two homogeneous phases: one micellar phase and one water phase.21 Anisotropy is thereby neglected. The purpose of this study is to prove the sensitivity of COSMOmic14 toward different micellar structures provided by MD simulations. Therefore, in the case of SDS more than 4600 and for CTAB more than 800 single micelles have been studied.

2. MODELS AND METHODS In this work, all-atom MD simulations have been used to study surfactant self-assembly. Subsequently, the obtained micelle structures were used for partition behavior predictions with the model COSMOmic, thereby the two methods were combined. The fundamentals of the models are described in the following sections as well as the analysis methods. 2.1. Molecular Dynamics Simulations. The open-source molecular dynamic package GROMACS 4.5.322,23 was chosen. The CHARMM3624 force field was used due to its good parametrization (e.g., agreement with experimental surface area per lipid is on average within 2% and the density profiles agree well with neutron and X-ray diffraction experiments) of amphiphilic molecules like lipids.25,26 The surfactant SDS is already parametrized in the force field CHARMM36 and was used without any changes. The surfactant CTAB was not directly included in the force field and was build from parameters of related lipid structures (DSPE and LPPC), which were already included in the force field. In the case of SDS and CTAB, the counterions sodium and bromide have been added to the system for reaching charge neutrality. Since bromide is not included in CHARMM36, the parametrization has been taken from Horinek et al.,27 who used the CHARMM general force field. To obtain an initial configuration for the MD simulations, a particular amount of surfactant monomers was placed into a cubic box,

X X μMX (r X , dX) = μM,CRS (r X , dX) + μM,comb (r X , dX)

(1)

X

The radial coordinate of X is denoted as r , and the orientation of X with respect to the radial vector of the micellar system by a unit vector dX. In all COSMOmic calculations, an infinite dilution of the solute is assumed. Thereby no interaction between solute molecules is taken into account. The part μXM,CRS of the chemical potential is usually the dominating part, especially for polar compounds. By the summation of the local segment chemical potentials, this part is responsible for favoring such 11583


Langmuir

Article

Figure 1. Surfactant structures of SDS (left) and CTAB (right). Counterions are presented only along with the chemical structure. Blue circles in the chemical structure represent the carbon atoms selected for the micelle definition via distances.

Table 1. System Configurations of SDS and CTAB at Two Different Concentrations no. of molecules

a

surfactant system

surfactants/ions

water

cmcexperimental [mol/L]

cSurfactant [mol/L]

no. of simulations

NPT simulation [ns]

SDS SDS CTAB CTAB

216 216 216 216

8535 25 605 12 000 120 000

8.0 × 10−3ab/9.0 × 10−3c

1.00 0.41 0.73 0.10

4 2 2 1

100 100 100 100

9.0 × 10−4a/10.0 × 10−4de

Reference 48. bReference 49. cReference 50. dReference 51. eReference 52.

positions and orientations of a solute X, in which the more polar parts are located in more polar regions of the micellar system and less polar parts in less polar regions, respectively. The combinatorial part μXM,comb takes into account size and shape ratios of the solute and solvent molecules. The parts describing the elastic deformation free energy or the long-ranged electrostatic potentials have been studied by Klamt el al.14 as well and were proven to be negligible for neutral species. The total chemical potential, calculated with eq 1, can then be used to calculate the partition function of a solute X with respect to all n layers (center positions) and all m orientations within one layer as follows: n X ZM =

n

m

∑ ZMX(ri) = ∑ ∑ i=1

i=1 j=1

X ⎧ (r , d j) ⎫ ⎪ −μ ⎪ M i ⎬ exp⎨ ⎪ ⎪ kT ⎭ ⎩

KMX =

V (rn) n

pMX (ri) =

w (rn)

pMX (rn)

(5)

w

(2)

P MW =

ng NAV (rn) KMX nsurfactantMsurfactant n w (rn)

(6)

NA is the Avogadro constant, nsurfactant is the amount of surfactant molecules of the micelle, and Msurfactant the molecular mass of a surfactant. For further details about the models COSMO-RS and COSMOmic, please refer to the corresponding publications of Klamt and co-workers.14−20 2.3. Micellar Systems. In Figure 1, surfactant structures and their conformations, used in COSMOmic, as well as the sigma-surfaces are presented. In case of aqueous surfactant systems, micelles are assembling and disassembling and therefore the surfactant conformers undergo changes. To select a representative surfactant conformer for COSMOmic calculations, the SAS of all surfactant molecules within a micelle were calculated and thus an average SAS was obtained. Then the conformers with a SAS close to the average were chosen.45 The tool g_sas of GROMACS was used for the SAS calculations.47 The

(3)

From the probability, the free energy profile of the solute in the micellar system with respect to the bulk water phase can be calculated: GMX(ri) = − kT ln pMX (ri)

ngw

Whereby V(ri) denotes the volume of layer i, nw(ri) is the amount of water molecules in layer i, and nwg represents the total amount of water molecules in the system. The partition coefficient KXM gives the relation between the probabilities of a solute to be located in the micelle or the water phase. The system size dependent coefficient can be converted into a system size independent partition coefficient PMW [Lwater/kgsurfactant], that can be directly compared to experimental data:

In eq 2, k denotes the Boltzmann constant and T the temperature. The probability to find solute X in layer i can be written as

Z MX(ri) Z MX

⎛ n (r ) ⎞⎞ n ⎛ ∑i = 1 ⎜V (ri)pMX (ri) − V (rn)pMX (rn)⎜ nw i ⎟⎟ ⎝ w(rn) ⎠⎠ ⎝

(4)

In this work, the whole free energy profile along the micelle has been taken into account for the calculation of the partition coefficient of the solutes. The system size dependent partition coefficient KXM [molSurfactant/molwater] can be calculated as follows:2,14 11584


Langmuir

Article

SAS of the selected SDS conformer is 4.59 nm2 and for the CTAB conformer 5.66 nm2. For a MD simulation of a specific concentration, a defined amount of the surfactant monomers has to be placed into a cubic waterbox for reaching a specific concentration. All systems studied within this paper are given in Table 1. The self-assembly of SDS micelles has been simulated at concentrations of c1,SDS = 1.00 mol/L and c2,SDS = 0.41 mol/L. The high concentration was also used in the simulation study of Sammalkorpi et al.10 Also the self-assembly of CTAB was studied at two different concentrations: c1,CTAB = 0.73 mol/L and c2,CTAB = 0.10 mol/L. The first relatively high concentrations of SDS and CTAB have been chosen in order to reduce the simulation time until stable micelle sizes are reached, neglecting the poor solubility of the surfactants in water. The second lower concentrations have been selected to investigate the concentration influence on the self-assembly of the surfactants. 2.4. Analysis Methods. In order to distinguish between several micelles, a definition method for the micelles boundary is necessary. One efficient way is the definition via distances between the center of mass of surfactant monomers and selected carbon atoms. If the monomers come close enough to each other during the simulation, they are considered to be a micelle. The micelles are defined via distances between certain carbon atoms and the center of mass of single surfactants. This procedure was developed and validated by Sammalkorpi et al.10 and is adapted for the surfactants and models used in this work. For each surfactant, three to four different distances and cutoffs (rc1−rc4 [nm]) were selected, depending on the molecule length (SDS, three distances and cutoffs; CTAB, four distances and cutoffs). As described by Sammalkorpi et al.,10 the distance between the center of mass of two surfactants is always the first distance and the other distances are between carbon atoms. The last distance is the one between the terminal carbon atoms of the hydrophobic tails. Micelles are identified, if at least one of the distances is smaller than the first and smallest cutoff rc1 or, if two distances are smaller than the second cutoff rc2 and so on, while rc1 < rc2 < rc3 < rc4. It is important to note that the different distances are not related to a specific cutoff, so that all distances are checked against all cut-offs. In the case of SDS, the distances between the center of mass (COM), C5, and C12 (terminal carbon atom) of all surfactants have been calculated and in case of CTAB the distances between the COM, C5, C9, and C16 (terminal carbon atom) have been determined (compare Figure 1). The choice of the cutoffs has an influence on the size of the micelles and is discussed in more detail in section 3.1. The radial density profiles of the micelles are calculated with respect to the center of mass of the micelles. From these profiles, as well as from the radius of gyration Rg, the radius of the micelles Rs can be estimated as proposed by Bogusz et el.:53 Rs =

5 Rg 3

of SDS and 16 C-atoms for CTAB and comparable small head groups (see Figure 1). Sammalkorpi et al.10 found the three cutoff distances of rc1 = 0.45, rc2 = 0.50, and rc3 = 0.60 nm for SDS at cSDS= 1 mol/L for an united atom parametrization. Using the cutoffs of Sammalkorpi et al.,10 many micelles, which can visually easily be identified as those, are not considered to be aggregated according to these cutoffs (using the parametrization of the CHARMM36 force field with an all-atom parametrization). For this reason, we have adapted this method for different cutoffs. However, the reference atoms and the micelle center of mass are kept the same as defined by Sammalkorpi et al.10 In Figure 2, results using three different cutoffs are presented in order to show the sensitivity of the aggregation numbers on

Figure 2. Probabilities of aggregation numbers (micelle sizes) for different cutoff definitions rc1−rc3 (red, rc1 = 0.45, rc2 = 0.50, rc3 = 0.60; blue, rc1 = 0.52, rc2 = 0.57, rc3 = 0.67; orange, rc1 = 0.61, rc2 = 0.66,rc3 = 0.76); all cutoffs are given in [nm].

the cutoff definition. The red curve results from the smallest cut-offs investigated (rc1 = 0.45; rc2 = 050; rc3 = 0.60 [nm]), which leads to a broad distribution of the probability of different aggregation numbers and the free monomer concentration is the highest. These are the same cutoffs as used by Sammalkorpi et al.10 The blue curve results from larger cutoffs and thereby more defined peaks of the aggregation numbers and a lower free monomer concentration. The third curve (orange) represents the largest cutoff analyzed, the most pronounced peaks and lowest free monomer concentration are obtained. In order to decide which cutoffs are most suitable for the definition of the micelles, an optimization criteria had to be defined. In this paper, we propose the critical micelle concentration (cmc) as a relevant parameter to estimate reasonable cut-offs. All free surfactants that are not incorporated within a micelle are taken into account in the cmc calculation. In literature, sometimes even oligomers of surfactants are taken into account for the cmc calculation.10,11 The cutoffs for SDS and CTAB optimized in this way are presented in Table 2. In Figure 3, the free monomer concentration of two NPT simulations of SDS at c1,SDS = 1.00 mol/L and c2,SDS = 0.41 mol/L and one NPT simulation of CTAB at c1CTAB = 0.73 mol/ L are presented and the experimental range of the cmc (compare Table 1) is marked. The free monomer concentrations result from applying the cutoffs presented in Table 2. The graph illustrates that the highest amount of free monomers in the solution is present at the beginning of the NPT simulations. After ∼20 ns the concentration fluctuates

(7)

The shape and stability of micelles can be further analyzed by examing the eccentricity ε,54 which is defined as:

ε=1−

Imin Iavg

(8)

Imin is the moment of inertia along the principal axes with the smallest magnitude, whereas Iavg is the average of all three moments of inertia. If ε = 0 or Imax/Imin = 1, the shape of the micelle is a perfect sphere, and if ε = 1, the shape is ellipsoid.4,54

3. RESULTS AND DISCUSSION 3.1. Definition of Micelles Using the Free Monomer Concentration. As mentioned earlier in section 2.4, the definition of the micelles is a crucial step for the analysis process. Therefore, the sensitivity of the aggregation sizes on the cut-offs has been determined and a universal optimization method has been developed for SDS and CTAB. Both surfactants have a similar tail length, 12 C-atoms in the case 11585


Langmuir

Article

Table 2. Optimized Cutoffs for the Different Surfactants SDS and CTAB cutoffs [nm] sufactant system SDS CTAB

distances COM COM

C5 C5

C12 C9

C16

rc1

rc2

rc3

rc4

0.61 0.70

0.66 0.74

0.76 0.89

0.89

Figure 3. (a) Free monomer concentration over simulation time for different simulations of SDS (green, c1,SDS= 1.00 mol/L; blue, c2,SDS = 0.41 mol/L; red, range of experimental cmc48−50). (b) Free monomer concentration over simulation time for CTAB (gray, c1,CTAB = 0.73 mol/L; red, range of experimental cmc).48,51,52

around the experimental cmc. Obviously, the free monomer concentration in the solution is almost independent of the total concentration of surfactant, so that the cmc is reached both for c1,SDS and for c2,SDS within the simulation time of 100 ns. This leads to the conclusion, that the cmc can be used as a criterion for the cutoff definition. The method has also been successfully applied to the system of CTAB (Figure 3b) and optimized cutoffs have been applied for all analysis (see Table 2). In contrast to our findings, LeBard et al.11 found a decrease in free surfactant at increasing surfactant concentration using coarse grained models in the concentration range from cSDS = 0.13−1.005 mol/L, which can be explained by the increase of ionic strength as the overall surfactant concentration is increased.55,56 However, the authors state, that the effect is surprisingly strong. This concentration effect on the free monomer concentration was not observed during our work and may be explained by the different methods and models applied (CHARMM36, all-atom MD, etc.) or the use of the optimized cut-offs for the micelle definition. 3.2. Self-Assembly of Surfactants. In Figure 4, the selfassembly of SDS, namely, the progress of the number of micelles (N), the maximum size (maximum aggregation number), and the average size (average aggregation number) are shown. In the system with the higher concentration (top figure,) the aggregation process is much faster than the process in the system with the lower concentration (bottom figure). Furthermore, the number of micelles is higher for the lower concentration, whereas the average and maximum aggregation numbers (Nagg) are much smaller. The average aggregation number in the system with the higher SDS concentration c1,SDS = 1.00 mol/L fluctuates around 45−55 and the maximum aggregation number is around 80− 85. In comparison, the values for the SDS system with c2,SDS = 0.41 mol/L are 20−25 and 30−38, respectively. In the SDS system with the higher concentration, the values for size and number of micelles seem to become almost constant after the

Figure 4. Self-assembly of SDS at different concentrations: c1,SDS = 1.00 mol/L averaged over multiple simulations (MR1-MR4) (top figure) and c2,SDS = 0.41 mol/L averaged over multiple simulations (MR1-MR2) (bottom figure); red, maximum aggregation number; green, average aggregation number; blue, number of aggregates.

first 20 ns, whereas it takes around 45 ns for the system with the lower concentration. Our findings are in correspondence with literature data,11 where Nagg = 43 ± 8 for cSDS = 0.363 mol/L and Nagg = 113 ± 14 for cSDS = 1.005 mol/L have been observed with a coarse grained model and simulation times in the range of μs. Sammalkorpi et al.10 studied SDS at 1 mol/L as well, but with a united atom force field and retrieved average sizes of Nagg = 60−70 for simulations of 200 ns. Still it is possible that even in large-scale simulations the micelles could evolve to form bigger aggregates. For the self-assembly process of the cationic CTAB a similar concentration dependence as for the anionic SDS is observed as presented in Figure 5. In case of the high concentration with c1,CTAB = 0.73 mol/L the maximum aggregation number is already large (Nagg ≥ 110) after a very short time, resulting in a small amount of aggregates (25 aggregates after 5 ns) with a maximum aggregation number of Nagg = 23 are present. The self-assembly of CTAB is faster for higher 11586


Langmuir

Article

1.00 mol/L as well as for c2,SDS = 0.41 mol/L. Figure 6 illustrates the results of four independent simulations each 100 ns long for c1,SDS = 1.00 mol/L and for a single run at c2,SDS = 0.41 mol/L. Obviously, the aggregation numbers are not the same for all simulations (named as MR1-MR4), which were started from different start configurations (black lines). The blue line indicates the average over all simulations at c1,SDS = 1.00 mol/L. Comparing the results of SDS at different concentrations, it can be noted that the aggregation numbers at higher concentration are larger (for the simulation time of 100 ns) than at lower concentration. Several distinguished peaks can be observed at higher concentration, whereby the probability distribution at the lower concentration is much more clustered in the region of smaller aggregation numbers. Recently, Anachkov et al.57 found experimentally an aggregation number of Nagg = 65 at cSDS = 0.10 mol/L, which corresponds to our data and as well to other literature data.10 The authors57 also investigated the change of the aggregation number with increasing surfactant concentration and could not find any significant changes above cSDS = 0.05 mol/L. Similar findings as for SDS can also be observed for CTAB, as shown in Figure 7. Due to the fact, that at the lower

Figure 5. Self-assembly of CTAB at different concentrations: c1,CTAB = 0.73 mol/L averaged over multiple simulations (MR1-MR2) (top figure) and c2,CTAB = 0.10 mol/L (bottom figure); red, maximum aggregation number; green, average aggregation number; blue, number of aggregates.

concentrations compared to lower concentrations, analogous to SDS. Thus, in order to achieve realistic micelle sizes in a reasonable time-scale for predictions with COSMOmic, higher concentrations are advantageous. In section 3.4, it will be shown that the COSMOmic predictions are quite insensitive to the micelle size. From the MD simulations, the probabilities of the aggregation numbers can be calculated. In Figure 6, the distribution of those probabilities is shown for SDS at c1,SDS =

Figure 7. Probabilities of the aggregation numbers for CTAB at different concentrations: black, probabilities for c1,CTAB = 0.73 mol/L of 2 independent simulations (MR1-MR2); red, probability for c2,CTAB = 0.10 mol/L of a single simulation.

concentration (c2,CTAB = 0.10 mol/L) 10 times more water is present, the aggregation process is much slower resulting in much smaller aggregates in the same simulation time compared to results at higher concentration. In the literature, the aggregation numbers for CTAB in dependence of surfactant concentration were published as follows: Padalkar et al.58 received for cCTAB = 0.10 mol/L Nagg = 156, Anachkov et al.57 found for cCTAB = 0.05 mol/L Nagg = 135, and Sánchez and Ruiz59 for cCTAB = 0.02 mol/L Nagg = 81. Hence, the equilibrium conditions may not have been reached for CTAB in the first 100 ns of the simulation, because the micelles observed during our MD simulation are smaller than the values reported. However, these micelles were taken for the calculation of the partition coefficient to prove the influence of the size of the micelles on the predictions of the partition behavior of different solutes (section 3.4). 3.3. Density Profiles. For further analysis of the micelle structures, several micelles have been selected according to their highest probabilities over the simulation time (see Figures 6 and 7). Only heavy atoms are taken into account when

Figure 6. Probabilities of the aggregation number for SDS at different concentrations: black, probabilities for c1,SDS = 1.00 mol/L of 4 independent simulations (MR1-MR4); blue, average probabilities for c1,SDS = 1.00 mol/L of 4 independent simulations (MR1-MR4); red, probability for c2,SDS = 0.41 mol/L of single simulation. 11587


Langmuir

Article

agreement to literature,6,10,60 although in most publications the profiles are given for SDS micelles with an aggregation number of 60 (mean aggregation number observed in experiments). In Figure 8, it can be recognized, that the density of different micelle sizes and simulations do not differ significantly, but the profiles become broader for larger micelles with Nagg ≥ 70 compared to smaller micelles. The highest density of the tail groups can be found between 0.9 and 1.3 nm and for the head groups in the range of 1.8−2.2 nm, whereby the head groups are accompanied by counterions (not shown in the graph). In case of SDS individual simulations were analyzed, for CTAB we analyzed different parts of one simulation to investigate structural changes over simulation time. In Figure 9, the density profiles of the CTAB micelles, which are formed in the intervals 0−10, 40−50, and 90−100 ns are presented in comparison for Nagg = 97−98 and Nagg = 109−110. The profiles indicate that micelles formed at the beginning of the MD simulation are similar to those formed at 100 ns. This is due to the fact that micelles, once formed during the simulation, are unlikely to fall apart in the time scale of 100 ns. In the case of CTAB, the head groups are even smaller in relation to the tail groups, since the tail groups consist of 16 carbon atoms instead of 12 carbon atoms for SDS. The trend of the profiles is comparable to those of SDS. The highest density for the tail groups can be found at 1.3−1.7 nm and for the head groups at 2.3−2.8 nm. As mentioned in section 3.2, on the one hand, the equilibrium conditions are not reached during the 100 ns simulation time, but on the other hand no significant changes occur in terms of micelle structure either. 3.4. Influence of the Micelle Structure on the Predictions with COSMOmic. The micelles, which have been obtained during the MD simulations, have been selected as input structures for calculations with COSMOmic to predict the partition behavior of different solute molecules. For SDS and CTAB the solute molecules toluene, propanolol, ephedrine, acetone, phenol, lidocaine, syringic acid, coumarin, isovanillin, ferulic acid, and vanillic acid are chosen in their neutral state. All experimental data of the partition coefficients has been taken from the study of Mehling et al.1 The parameters of SDS and CTAB for the eccentricity ε, the radius of gyration Rg, and micelle radius Rs as well as the RMSE (root mean square error, as defined in eq 9) of all log PMW for micelles of different sizes are presented in Table 3. The RMSE is defined as

calculating the radial density profiles with respect to the center of mass of the micelle (COM) as presented in Figures 8 and 9.

Figure 8. Density profile of SDS micelles at c1,SDS = 1.00 mol/L for 4 independent simulations, MR1−MR4 over 100 ns each. The profiles are averaged over different micelle size groups (MR1: Nagg = 57−59 and Nagg = 71−73; MR2: Nagg = 59−61 and Nagg = 72−74; MR3: Nagg = 54−56; MR4, Nagg = 54−56).

Figure 9. Density profile of CTAB micelles at c1 = 0.73 mol/L for the simulation time intervals 0−10, 40−50, and 90−100 ns.

The profiles of SDS are averaged over three different micelle size groups according to their probability (MR1, Nagg = 57−59 and Nagg = 71−73; MR2, Nagg = 59−61 and Nagg = 72−74; MR3, Nagg = 54−56; MR4, Nagg = 54−56). The tail groups are pointing toward the center of the micelle, whereas the head groups are pointing into the water phase as expected. The structure of the tail group of SDS is very close to dodecane. So the radial density profile of SDS has been validated by comparison of the tail group density with the density of a dodecane droplet in water prior to further analysis (data not shown). The tail group of SDS is much larger than the headgroup, which results also in a higher density of the tail groups compared to the head groups. This trend is in

RMSE =

1 n

n

∑ (log P MW,COSMOmic − log P MW,exp)2 i=1

(9)

In case of SDS more than 4600 and for CTAB more than 800 single micelles have been selected in order to achieve

Table 3. Structure Parameters of SDS and CTAB Micelles and RMSE of log PMW surfactant system SDS

CTAB

no. of micelles

shape

ε

Rg [nm]

Rs [nm]

RMSE of log PMW

55−57 71−73 92−94

1470 1229 1970

spheric spheric spheric

0.13 ± 0.06 0.12 ± 0.04 0.18 ± 0.06

1.60 ± 0.06 1.73 ± 0.02 1.92 ± 0.07

2.07 ± 0.07 2.23 ± 0.03 2.48 ± 0.09

0.44 0.43 0.42

97−98 97−98 109−110

23 575 222

cylindric spheric spheric

0.59 ± 0.04 0.13 ± 0.06 0.18 ± 0.06

2.36 ± 0.08 1.99 ± 0.03 2.08 ± 0.06

3.04 ± 0.10 2.57 ± 0.03 2.68 ± 0.08

0.61 0.43 0.41

Nagg

11588


Langmuir

Article

COSMOmic, the micelles formed during the MD simulations, need to be analyzed and selected carefully according to shape and structure. The log PMW,COSMOmic values of SDS and CTAB micelles for individual solutes are presented in Tables 4 and 5 for different

statistically reliable results. Since the COSMOtherm software (including the COSMOmic package) is designed for only one input structure, we have developed an approach for automatic calculations of the partition coefficients for a large number of micelles. Hence, every micelle was used in a single calculation and the results were averaged afterwards. Note, that the averagese were calculated for different groups of micelle sizes (see Table 3). In order to reach charge neutrality for the systems used for the COSMOmic calculations for each selected micelle/water system, the amount of counterions corresponds to the amount of surfactants in the micelle. Note that some of the layers of the micelle may be charged during the COSMOmic calculation and the last layer is a pure bulk water layer in order to calculate the micelle/water partition coefficient. This assumption has been validated by further MD simulations, showing that the highest radial density of the counterions is located close to the micelles and nearly no ions can be found in the bulk water phase (data not shown). The structure parameters ε, Rg, and Rs of the SDS micelles of size Nagg = 55−57 obtained in this work are similar to those values found in literature: Bruce et al.5 reported an average value of Imax/Imin = 1.05 (spherical micelle) and for Rg = 1.62 ± 0.012 nm (Rs = 2.09 ± 0.015) for a simulation of 60 SDS molecules and 7579 water molecules (cSDS = 0.4 mol/L), MacKerell4 studied a similar system with 60 SDS molecules and 4398 water molecules and retrieved Imax/Imin = 1.02 (spherical micelle) and Rg = 1.602 nm. Palazzesi et al.60 found an eccentricity factor for a SDS micelle with an aggregation number of Nagg = 60 of ε = 0.154 averaged over the simulation time, which is close to our values of ε = 0.12−0.18. Comparing the different sizes of SDS micelles (see Table 3), it becomes obvious that Rg and Rs increase significantly with size, as expected. In case of larger SDS micelles with an aggregation number of Nagg ≥ 92, the micelles become more ellipsoid, what is indicated by the increase in the eccentricity factor from ε = 0.12 to ε = 0.18. The results presented in Table 3 for spherical CTAB micelles are in line with the data shown for SDS, since the CTAB micelles have larger aggregation numbers and bigger monomers, the radii of the micelles are larger. Anachkov et al.57 published a hydrodynamic radius of 2.85 nm for CTAB with Nagg = 135, which corresponds to our values. Comparing the RMSE of all log PMW values presented in Table 3, no dependency on the aggregation number can be observed for spherical micelles with ε 0.5).

Comparing the results of the spherical micelles with RMSE = 0.41−0.43 with the cylindrical micelles with RMSE = 0.61, the prediction quality appears to decrease for cylindrical micelles. However, the larger error for cylindrical micelles can also be assigned to the much smaller amount of cylindrical micelles and the larger influence of a single micelle on the overall RMSE. It has to be noted that not all solute molecules can be predicted with the same accuracy with COSMOmic for SDS and CTAB (see both Tables 4 and 5). In Figure 10, the RMSE of the log PMW,COSMOmic for different groups of micelle sizes is presented (Nagg = 55 (blue), Nagg = 71 (red), and Nagg = 94 (black)). Thereby more than 1600 micelles have been analyzed 11589


Langmuir

Article

Figure 10. Prediction quality of COSMOmic for individual SDS micelles with Nagg = 55 (blue), Nagg = 71 (red), and Nagg = 94 (black).

in this graph. Analyzing the micelles in detail, each micelle differs in structure and shape, therefore resulting in different partition coefficients of the solutes. According to Figure 10, it seems that small changes in the structure of the micelle can lead to significant changes in the predicted partition behavior. However, tests, in which the number of layers (COSMOmic input, see section 2.2) are varied, show that the prediction quality of the outliers can be improved to reach similar RMSE values as for the majority of the micelles. Based on this analysis, it becomes obvious that, due to the high probability of the outliers, it is not sufficient to use only one micelle structure as an input for COSMOmic. In order to get representative values of partition coefficients, several structures from MD simulations should be taken. However, they can be randomly chosen, since no specific trends are observed. Nevertheless, COSMOmic not only predicts partition coefficients, but the whole free energy profiles for each solute in the micelle/water system at infinite dilution (no interaction of solute molecules) as well. Note, that also MD simulations can be used to predict free energy profiles. However, as these MD simulations (e.g., Umbrella Sampling) are extremly demanding, COSMOmic is the only model up to now that can calculate these profiles for a large amount of solutes on a single standard computer in reasonable time. The free energy profiles predicted by COSMOmic are presented in Figure 11. The profiles indicate potential energy barriers at the head group region of the micelles and the most favored location of the solutes in the micelle/water system, which is the location of the free energy minimum. The favored location of the solute molecules in SDS micelles is shifted toward the center of mass of the micelle (COM) for smaller micelles (left graphs of Figure 11). Therefore, these free energy profiles give additional information on the most probable location of the solute in the system, which is not accessible by analyzing only the partition coefficients. In case of CTAB no significant differences between Nagg = 97 and Nagg = 110 can be observed, indicating no significant influence of the difference of size in this range on the predictions with COSMOmic. However, the difference in size was much smaller for the CTAB micelles investigated compared to the SDS micelles. Hence, the relative locations and trends of the profiles of all three solutes phenol, lidocaine and syringic acid are comparable in SDS and CTAB micelles, since all solutes are of similar hydrophobicity.

Figure 11. Free energy profiles calculated with COSMOmic for SDS micelles with Nagg = 55 and Nagg = 94 and CTAB micelles with Nagg = 97 and Nagg = 110 for the solutes phenol, lidocaine, and syringic acid; all solutes are in their neutral charge state.

4. CONCLUSIONS In this paper, the self-assembly of the amphiphilic molecules SDS and CTAB in aqueous solutions at different concentrations were studied with all-atom MD simulations. For the definition of micelles we presented a method based on the free monomer concentration in the aqueous solution as an optimization criterion, which was validated for SDS and CTAB at high and low concentrations. During the MD simulations micelles are formed spontaneously in the first nanoseconds of the simulations, but equilibrium conditions are not reached in the time scale of 100 ns. The results are discussed concerning the size, shape, and radial density profiles, since different sizes of micelles are formed with high probabilities. All results are in good agreement with the simulations and experiments reported in literature.6,10,57−60 Furthermore, the micelles produced by MD simulations were used as an input for the model COSMOmic in order to predict the partition behavior of various solutes. Thereby the influence of the micellar structures on the partition behavior as described by COSMOmic was studied. In the case of SDS more than 4600 and for CTAB more than 800 single micelles have been selected in order to achieve statistical reliable results. Therefore, an approach to calculate automatically the partition coefficients of large numbers of micelles has been developed. By comparing the predictions of the partition coefficients of various solutes with experimental data,1 no strong influence of the micelle size was observed; however, outliers are probable. Thus, we can advice to use not a single micelle but several different micelle structures as an input in COSMOmic and to use an average as a 11590


Langmuir

Article

coarse-grained ionic surfactants accelerated by graphics processing units. Soft Matter 2012, 8, 2385. (12) Catá, G. F.; Rojas, H. C.; Gramatges, A. P.; Zicovich-Wilson, C. M.; Á lvarez, L. J.; Searle, C. Initial structure of cetyltrimethylammonium bromide micelles in aqueous solution from molecular dynamics simulations. Soft Matter 2011, 7, 8508. (13) Wang, Y.; Larsson, D. S. D.; van der Spoel, D. Encapsulation of Myoglobin in a Cetyl Trimethylammonium Bromide Micelle in Vacuo: A Simulation Study. Biochemistry 2009, 48, 1006−1015. (14) Klamt, A.; Huniar, U.; Spycher, S.; Keldenich, J. COSMOmic: A Mechanistic Approach to the Calculation of Membrane−Water Partition Coefficients and Internal Distributions within Membranes and Micelles. J. Phys. Chem. B 2008, 112, 12148−12157. (15) Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. J. Phys. Chem. 1995, 99, 2224−2235. (16) Klamt, A.; Jonas, V.; Bürger, T.; Lohrenz, J. C. W. Refinement and Parametrization of COSMO-RS. J. Phys. Chem. A 1998, 102, 5074−5085. (17) Klamt, A.; Eckert, F. COSMO-RS: a novel and efficient method for the a priori prediction of thermophysical data of liquids. Fluid Phase Equilib. 2000, 172, 43−72. (18) Klamt, A.; Eckert, F.; Diedenhofen, M. Prediction or partition coefficients and activity coefficients of two branched compounds using COSMOtherm. Fluid Phase Equilib. 2009, 285, 15−18. (19) Klamt, A.; Eckert, F.; Arlt, W. COSMO-RS: An Alternative to Simulation for Calculating Thermodynamic Properties of Liquid Mixtures. Annu. Rev. Chem. Biomol. Eng. 2010, 1, 101−122. (20) Klamt, A. The COSMO and COSMO-RS solvation models. Wiley Interdiscip. Rev.: Comput. Mol. Sci. 2011, 1, 699−709. (21) Mokrushina, L.; Buggert, M.; Smirnova, I.; Arlt, W.; Schomäcker, R. COSMO-RS and UNIFAC in Prediction of Micelle/ Water Partition Coefficients. Ind. Eng. Chem. Res. 2007, 46, 6501− 6509. (22) Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith, J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a high-throughput and highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845−854. (23) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. GROMACS: Fast, flexible, and free. J. Comput. Chem. 2005, 26, 1701−1718. (24) MacKerell, A. D.; Brooks, B.; Brooks, C. L., III; Nilsson, L.; Roux, B.; Won, Y.; Karplus, M. CHARMM: The Energy Function and Its Parameterization. In Encyclopedia of Computational Chemistry; John Wiley and Sons, Ltd.: 2002. (25) Pastor, R. W.; Mackerell, A. D. Development of the CHARMM Force Field for Lipids. J. Phys. Chem. Lett. 2011, 2, 1526−1532. (26) Klauda, J. B.; Venable, R. M.; Freites, J. A.; O’Connor, J. W.; Tobias, D. J.; Mondragon-Ramirez, C.; Vorobyov, I.; MacKerell, A. D.; Pastor, R. W. Update of the CHARMM All-Atom Additive Force Field for Lipids: Validation on Six Lipid Types. J. Phys. Chem. B 2010, 114, 7830−7843. (27) Horinek, D.; Mamatkulov, S. I.; Netz, R. R. Rational design of ion force fields based on thermodynamic solvation properties. J. Chem. Phys. 2009, 130, 124507. (28) Martyna, G. J.; Tuckerman, M. E.; Tobias, D. J.; Klein, M. L. Explicit reversible integrators for extended systems dynamics. Mol. Phys. 1996, 87, 1117−1157. (29) Nosé, S.; Klein, M. Constant pressure molecular dynamics for molecular systems. Mol. Phys. 1983, 50, 1055−1076. (30) Parrinello, M. Polymorphic transitions in single crystals: A new molecular dynamics method. J. Appl. Phys. 1981, 52, 7182. (31) Durell, S. R.; Brooks, B. R.; Ben-Naim, A. Solvent-Induced Forces between Two Hydrophilic Groups. J. Phys. Chem. 1994, 98, 2198−2202. (32) Mackerell, A. D.; Bashford, D.; Bellott; Dunbrack, R. L.; Evanseck, J. D.; Field, M. J.; Fischer, S.; Gao, J.; Guo, H.; Ha, S.; Joseph-McCarthy, D.; Kuchnir, L.; Kuczera, K.; Lau, F. T. K.; Mattos,

representative result. Outlier analysis might be also needed in each specific case. COSMOmic has been successfully used to calculate free energy profiles of solutes in micelles within reasonable calculation time. Free energy profiles gave additional information on the most probable location of the solute in the system, which is not accessible by analyzing only the partition coefficients. These facts underline the benefits of the combination of micelles obtained from MD simulations with COSMOmic and its suitability as a screening tool for solute partition behavior.

■

AUTHOR INFORMATION

Corresponding Author

*Telephone: +49 (40) 42878 2988. Fax: +49 (40) 42878 4072. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors appreciate the financial support of the DFG (Project SM 82/4-2) and the Department of Energy, Office of Basic Energy Sciences, Award DE-SC0002128. Computational resources have been provided by the North-German Supercomputing Alliance (HLRN).

■

REFERENCES

(1) Mehling, T.; Kloss, L.; Ingram, T.; Smirnova, I. Partition Coefficients of Ionizable Solutes in Mixed Nonionic/Ionic Micellar Systems. Langmuir 2012, 29, 1035−1044. (2) Ingram, T.; Storm, S.; Kloss, L.; Mehling, T.; Jakobtorweihen, S.; Smirnova, I. Prediction of Micelle/Water and Liposome/Water Partition Coefficients Based on Molecular Dynamics Simulations, COSMO-RS, and COSMOmic. Langmuir 2013, 29, 3527−3537. (3) Mehling, T.; Ingram, T.; Storm, S.; Bobe, U.; Liu, F.; Michel, M.; Smirnova, I. Estimation of LPC/water partition coefficients using molecular modeling and micellar liquid chromatography. Colloids Surf., A 2013, 431, 105−113. (4) MacKerell, A. D. 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. (5) Bruce, C. D.; Berkowitz, M. L.; Perera, L.; Forbes, M. D. E. Molecular Dynamics Simulation of Sodium Dodecyl Sulfate Micelle in Water: Micellar Structural Characteristics and Counterion Distribution. J. Phys. Chem. B 2002, 106, 3788−3793. (6) Gao, J.; Ge, W.; Hu, G.; Li, J. From Homogeneous Dispersion to MicellesA Molecular Dynamics Simulation on the Compromise of the Hydrophilic and Hydrophobic Effects of Sodium Dodecyl Sulfate in Aqueous Solution. Langmuir 2005, 21, 5223−5229. (7) Sammalkorpi, M.; Sanders, S.; Panagiotopoulos, A. Z.; Karttunen, M.; Haataja, M. Simulations of Micellization of Sodium Hexyl Sulfate. J. Phys. Chem. B 2011, 115, 1403−1410. (8) Sanders, S. A.; Sammalkorpi, M.; Panagiotopoulos, A. Z. Atomistic Simulations of Micellization of Sodium Hexyl, Heptyl, Octyl, and Nonyl Sulfates. J. Phys. Chem. B 2012, 116, 2430−2437. (9) Piotrovskaya, E. M.; Vanin, A. A.; Smirnova, N. A. Molecular dynamics simulation of micellar aggregates in aqueous solution of hexadecyl trimethylammonium chloride with different additives. Mol. Phys. 2006, 104, 3645−3651. (10) 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. (11) LeBard, D. N.; Levine, B. G.; Mertmann, P.; Barr, S. A.; Jusufi, A.; Sanders, S.; Klein, M. L.; Panagiotopoulos, A. Z. Self-assembly of 11591


Langmuir

Article

(53) Bogusz, S.; Venable, R. M.; Pastor, R. W. Molecular Dynamics Simulations of Octyl Glucoside Micelles: Structural Properties. J. Phys. Chem. B 2000, 104, 5462−5470. (54) Salaniwal, S.; Cui, S. T.; Cochran, H. D.; Cummings, P. T. Molecular Simulation of a Dichain Surfactant/Water/Carbon Dioxide System. 1. Structural Properties of Aggregates. Langmuir 2001, 17, 1773−1783. (55) Moroi, Y. Micelles: Theoretical and applied aspects; Plenum Press: New York, 1992. (56) Zana, R. Dynamics of surfactant self-assemblies: Micelles, microemulsions, vesicles, and lyotropic phases; Taylor & Francis/CRC Press: Boca Raton, FL, 2005. (57) Anachkov, S. E.; Danov, K. D.; Basheva, E. S.; Kralchevsky, P. A.; Ananthapadmanabhan, K. P. Determination of the aggregation number and charge of ionic surfactant micelles from the stepwise thinning of foam films. Adv. Colloid Interface Sci. 2012, 183−184, 55− 67. (58) Padalkar, K.; Gaikar, V.; Aswal, V. Small angle neutron scattering studies of mixed micelles of sodium cumene sulphonate with cetyl trimethylammonium bromide and sodium dodecyl sulphate. Pramana 2008, 71, 953−957. (59) Sánchez, F. G.; Ruiz, C. C. Intramicellar energy transfer in aqueous CTAB solutions. J. Lumin. 1996, 69, 179−186. (60) Palazzesi, F.; Calvaresi, M.; Zerbetto, F. A molecular dynamics investigation of structure and dynamics of SDS and SDBS micelles. Soft Matter 2011, 7, 9148.

C.; Michnick, S.; Ngo, T.; Nguyen, D. T.; Prodhom, B.; Reiher, W. E.; Roux, B.; Schlenkrich, M.; Smith, J. C.; Stote, R.; Straub, J.; Watanabe, M.; Wiórkiewicz-Kuczera, J.; Yin, D.; Karplus, M. All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins †. J. Phys. Chem. B 1998, 102, 3586−3616. (33) Essmann, U.; Perera, L.; Berkowitz, M. L.; Darden, T.; Lee, H.; Pedersen, L. G. A smooth particle mesh Ewald method. J. Chem. Phys. 1995, 103, 8577. (34) Darden, T.; York, D.; Pedersen, L. Particle mesh Ewald: An N· log(N) method for Ewald sums in large systems. J. Chem. Phys. 1993, 98, 10089. (35) Hess, B. P-LINCS: A Parallel Linear Constraint Solver for Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 116−122. (36) Piggot, T. J.; Piñeiro, Á .; Khalid, S. Molecular Dynamics Simulations of Phosphatidylcholine Membranes: A Comparative Force Field Study. J. Chem. Theory Comput. 2012, 8, 4593−4609. (37) COSMOlogic GmbH & Co. KG - cosmologic.de. COSMOlogic, computational chemistry and fluid phase thermodynamics from chemical engineering to life science; http://www.cosmologic.de/index.php (accessed 27 May 2013). (38) Ahlrichs, R.; Bär, M.; Häser, M.; Horn, H.; Kölmel, C. Electronic structure calculations on workstation computers: The program system turbomole. Chem. Phys. Lett. 1989, 162, 165−169. (39) Schäfer, A.; Huber, C.; Ahlrichs, R. Fully optimized contracted Gaussian basis sets of triple zeta valence quality for atoms Li to Kr. J. Chem. Phys. 1994, 100, 5829. (40) Perdew, J. P. Density-functional approximation for the correlation energy of the inhomogeneous electron gas. Phys. Rev. B: Condens. Matter Mater. Phys. 1986, 8822−8824. (41) Becke, A. D. Density-functional exchange-energy approximation with correct asymptotic behavior. Phys. Rev. A: At., Mol., Opt. Phys. 1988, 3098−3100. (42) Eichkorn, K.; Treutler, O.; Oehm, H.; Haeser, M.; Ahlrichs, R. TZVP. Chem. Phys. Lett. 1995, 652−660. (43) HyperChem. HyperChem Profesional Overview; http://www. hyper.com/?TabId=361 (accessed 29 May 2013). (44) Buggert, M.; Cadena, C.; Mokrushina, L.; Smirnova, I.; Maginn, E. J.; Arlt, W. COSMO-RS Calculations of Partition Coefficients: Different Tools for Conformation Search. Chem. Eng. Technol. 2009, 32, 977−986. (45) Jakobtorweihen, S.; Ingram, T.; Smirnova, I. Combination of COSMOmic and molecular dynamics simulations for the calculation of membrane-water partition coefficients. J. Comput. Chem. 2013, 34 (15), 1332−1340. (46) Ingram, T.; Gerlach, T.; Mehling, T.; Smirnova, I. Extension of COSMO-RS for monoatomic electrolytes: Modeling of liquid−liquid equilibria in presence of salts. Fluid Phase Equilib. 2012, 314, 29−37. (47) Eisenhaber, F.; Lijnzaad, P.; Argos, P.; Sander, C.; Scharf, M. The double cubic lattice method: Efficient approaches to numerical integration of surface area and volume and to dot surface contouring of molecular assemblies. J. Comput. Chem. 1995, 16, 273−284. (48) Aguiar, H. B. de; Strader, M. L.; Beer, A. G. F. de; Roke, S. Surface Structure of Sodium Dodecyl Sulfate Surfactant and Oil at the Oil-in-Water Droplet Liquid/Liquid Interface: A Manifestation of a Nonequilibrium Surface State. J. Phys. Chem. B 2011, 115 (12), 2970− 2978. (49) Umlong, I.; Ismail, K. Micellization behaviour of sodium dodecyl sulfate in different electrolyte media. Colloids Surf., A 2007, 299 (1−3), 8−14. (50) Wu, C.; Zhou, S. Effects of surfactants on the phase transition of poly(N-isopropylacrylamide) in water. J. Polym. Sci., Part B: Polym. Phys. 1996, 34 (9), 1597−1604. (51) Ruiz, C. C.; Aguiar, J. Interaction, Stability, and Microenvironmental Properties of Mixed Micelles of Triton X100 and n -Alkyltrimethylammonium Bromides: Influence of Alkyl Chain Length. Langmuir 2000, 16, 7946−7953. (52) Mata, J.; Varade, D.; Bahadur, P. Aggregation behavior of quaternary salt based cationic surfactants. Thermochim. Acta 2005, 428 (1−2), 147−155. 11592


Molecular Dynamics Simulation of SDS and CTAB ... - ACS Publications

Recommend Documents