A Coarse-Grained Implicit Solvent Model for Poly(ethylene oxide

Oct 16, 2015 - A Coarse-Grained Implicit Solvent Model for Poly(ethylene oxide), CnEm Surfactants, and Hydrophobically End-Capped Poly(ethylene oxide)...
2 downloads 8 Views 4MB Size
Download PDF
Article pubs.acs.org/Macromolecules

A Coarse-Grained Implicit Solvent Model for Poly(ethylene oxide), CnEm Surfactants, and Hydrophobically End-Capped Poly(ethylene oxide) and Its Application to Micelle Self-Assembly and Phase Behavior Shihu Wang and Ronald G. Larson* Department of Chemical Engineering, University of Michigan, Ann Arbor, Michigan 48109, United States ABSTRACT: We parametrize a molecular dynamics force field for poly(ethylene oxide) (PEO) within the Dry Martini coarse-grained (CG) implicit solvent model and show that it yields similar distributions of bond lengths, angles and torsional angles, and a similar dependence of radius of gyration on PEO chain length as those obtained from atomistic and explicit-solvent CG simulations. We apply these new parameters to simulate self-assembly of long-chain single-endcapped CnEm surfactants and double-end-capped telechelic polymers CnEmCn, where n is the number of alkane carbons and m the number of ethylene oxide monomers. We determine that the average equilibrium micelle aggregation number at high concentrations is around 23 hydrophobes for C12E100, C12E100C12, and C12E200C12, consistent with experimental/ theoretical results. However, for EO block less than 200 monomers at 300 K, for example, for C12E45C12 and C12E136C12, telechelic polymers phase separate to form a polymer dense phase at low concentrations due to the formation of bridged flowerlike micelles, with critical concentrations for phase separation also in agreement with experimental values.

1. INTRODUCTION Poly(ethylene glycol) (PEG) or poly(ethylene oxide) (PEO) is widely used in many biomedical and technological applications due to its excellent solubility in aqueous solutions, low toxicity, and low immunogenicity. PEO substituents are commonly used to modify many biological molecules, such as oligonucleotides, peptides, and proteins, and to chemically decorate drug/gene delivery vehicles, to reduce their innate immune response, via the PEGylation process.1,2 PEO is also an important building block for many chemical compounds of industrial significance. For instance, short PEG-based surfactants, such as alkyl poly(ethylene glycol) ethers, are among the most important types of nonionic surfactants and are extensively utilized as detergents, solubilizers, and emulsifiers.3 Long-chain poly(ethylene oxide)s of 1−100 kg/mol that are hydrophobically end-capped by C8 to C18 alkyl groups via urethane (isocynate) groups, known as hydrophobically modified ethylene oxide urethanes (i.e., “HEURs”), on the other hand, are widely used as rheology modifiers in water-borne coatings, paper coatings, oil recovery, and so on.4,5 The properties of PEO arise from its unique chemical structure with both hydrophobic and hydrophilic character and yet much more soluble than analogous polymers like poly(methylene oxide) and poly(propylene oxide).6 Because of its practical importance and scientific interest, the structural and functional properties of PEO and its derivatives have been extensively studied.6−16 From a computational perspective, different force fields have been put forward in © 2015 American Chemical Society

order to investigate the physical properties of PEO in solution.8−14 Although atomistic force fields9 provide the highest resolution, most atomistic simulations are limited to PEO with a low degree of polymerization in small simulation boxes. To simulate long chains, one normally resorts to coarsegrained force fields with explicit coarse-grained solvent to circumvent the time- and length-scale limitations of atomistic simulations. A common strategy in coarse-graining is to replace each ethylene oxide group with a single CG bead.10−12 For example, the Martini CG force fields (FF) parametrizes poly(ethylene oxide) in this way, and its predictions of conformations PEO with different molecular weights at various polymer concentrations show good agreement with atomistic simulations.10 The same mapping scheme, using a single CG bead for each PEO monomer, was also employed by Shinoda et al. (in the “SDK” FF),12 by Prasitnok et al.,13 and by Rossi et al.,11 all using different specific parametrizations. An even larger degree of coarse-graining was also tested by Chen8 and Wang,14 who grouped two EO repeat units into one CG bead. These different coarse-grained force fields differ from each other not only by their coarse graining and parametrization of PEO but and also by their coarse-grained mapping of other molecules, such as water. For example, in the Martini FF, one CG water Received: July 16, 2015 Revised: October 8, 2015 Published: October 16, 2015 7709

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules

makes transferability of its force field parameters easy. Specifically, Dry Martini FF adopts mainly a mapping of four heavy atoms to one CG bead. Different segments or repeat units of a chemical compound are modeled using different bead types that correspond to different levels of polarity, hydrogenbond capabilities, etc. As a result, new molecules can be easily constructed using these CG bead types as simple building blocks. By mixing different molecules, multicomponent systems are thereby readily constructed. In fact, the original Martini FF with explicit solvent that adopts the same strategy has been successfully applied to study a variety of molecules from lipids29 to carbonhydrates,30 proteins,31 and polymers.10,11,19,32 In this paper, we develop new bead parameters for PEO and combine them with some other bead types in a Dry Martini FF. We apply these parameters to simulate the self-assembly CnEm surfactants and HEURs. Using this implicit solvent model, we obtain coil dimensions for different lengths of PEO in dilute solution and compositional profiles of CnEm micelles that are similar to those obtained using explicit solvent models. We also investigate the self-assembly of long-chain CnEm surfactants and HEURs and show that the experimentally observed phase transition for HEURs is due to the formation of interconnected flower-like micelles. We note that the PEO parameters are transferable and can be used to simulate the adsorption of HEURs to surfaces or to model competitive adsorption of molecules to surfaces, for example, in the presence of sodium dodecyl sulfate (SDS) surfactants and HEURs.

bead represents four water molecules while it corresponds to only three water molecules in the SDK FF. The main advantage of coarse-graining is in reducing the numbers of degrees of freedom and the number of interactions. With softer potentials, and thus larger time steps and accelerated dynamics, the simulation speed can be significantly enhanced over that of atomistic force fields, by as much as 3 or even 4 orders of magnitude. As a result, larger system sizes and longer simulation times become affordable. Some of these CG force fields also combine well with other force field parameters.17−20 For example, PEO parameters from both Martini and SDK have been successfully combined with CG alkane FFs to study the self-assembly of CnEm surfactants17,18 (where n is the number of carbons and m is the number of ethylene oxides in the surfactant) and PEGylated lipids.19−21 But CG force fields with explicit solvent are still not fast enough to allow affordable simulations of many systems of interest, including many PEO-containing polymers and biomolecules. In fact, most CG simulations are still performed on short CnEm surfactants or (block) polymers of short lengths.17,18,20 Even within these restrictions simulations are expensive; it has been shown that using the SDK model on GPUs,18 obtaining an equilibrated micelle size distribution for C9E6, is still difficult. Additionally, as mentioned above, the backbones of typical HEURs consist of hundreds to thousands of repeating units of ethylene oxide, which is usually out of reach of even coarsegrained simulations with explicit solvent. To study the structural and dynamic properties of long PEO and PEO-based polymers and biomolecules, a further computational acceleration is needed. This can be achieved by removing the explicit solvent molecules and replacing their influence by implicit solvent potentials for the nonsolvent beads. The enhancement in speed thereby attained is especially great when the volume fraction of solvent is high, as in typical surfactant and polymer solutions. By omitting explicit solvent molecules, simulating dynamic processes essential to the equilibration, such as intermicelle surfactant exchange, micellar fission, and fusion, becomes possible.22 Both atomistic and coarse-grained PEO models with implicit solvent have been developed for some specific systems.23−27 For example, Chen et al.23 developed an implicit solvent model for polystyrene-b-poly(ethylene oxide) and studied its self-assembly kinetics using Brownian dynamics simulations over a millisecond time scale. Bedrov et al.24 investigated the structural properties of longchain poly(ethylene oxide)−poly(propylene oxide)−poly(ethylene oxide) micelles using coarse-grained implicit solvent simulations with parameters obtained from an inverse Boltzmann method. Fischer et al.25 developed coarse-grained potentials for PEO in implicit solvents using iterative Boltzmann inversion techniques and reproduced the experimentally measured radii of gyration for some long PEO chains containing up to 1500 monomers. Jusufi et al.27 employed an implicit solvent model and grand canonical Monte Carlo simulations to study micellization behavior and thereby determined the critical micelle concentration and aggregation number for short CnEm surfactants with several lengths of ethoxy and hydrocarbon tails at different temperatures. Although simulations with these force fields are not computationally demanding, the transferability of these implicit solvent models to new molecules or multicomponent mixtures is either unclear or limited. The recently developed Dry Martini FF28 with implicit solvent for biomolecules adopts an interesting approach that

2. SIMULATION DETAILS We perform simulations and analyses using the GROMACS33 simulation package 4.5.5. The force fields parameters are based on the Dry Martini force field developed by Marrink et al.28 However, we have added some new bead interactions as shown in Table 1 to obtain reasonable results involving ethylene oxide (EO), as discussed below. Specifically, we consider three different types of beads: C1, EO, and OA. Bead-type C1 represents four heavy carbon atoms (CCCC) and has a mass of 72 amu, and its interaction parameters are those of the established Dry Martini FF of the same name. EO and OA are newly parametrized beads in this study. The EO bead represents one ethylene oxide (COC) and OA one hydroxymethyl (COH) group. EO and OA each has a mass of 45 amu. The nonbonded interaction between any two CG beads is described using a simple 6−12 Lennard-Jones potential that is shifted to 0 from 0.9 to 1.2 nm. Two different values of the LJ distance parameter σ are used: σ = 4.3 Å for interactions between any two special types of beads, i.e., EO and OA, and σ = 4.7 Å for any interactions involving one or two regular CG beads, i.e., C1. We vary the energetic parameter ε for interactions involving EO and OA beads to get the best fit of the CG model with implicit solvent to the simulation results for the Martini CG model with explicit solvent for different CG bead−bead interactions. The final set of parameters is shown in Table 1. The systems we consider contain no charges, so no electrostatic interaction is considered. The bonded interactions are described by weak harmonic stretching and bending potentials.28 In the standard explicitsolvent Martini FF for PEO, different angle parameters are available based on either cosine or theta types of angle potential functions,10,17,19 and the theta-type angle potential function is preferred mainly because of improved numerical stability. Here for our implicit water model, we adopt the cosine-type angle potential function to be consistent with that used in the Dry 7710

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules

note that the results are not very sensitive to the particular piecewise linear/harmonic distance restraints. Here we select a simple linear harmonic distance restraint starting from 0 (r1 = 0). Below we compare results from our implicit water model with those from the standard Martini FF10,29 or the SDK FF,12,20 both of which are coarse-grained (CG) force fields with explicit water. For simulations using implicit water, an NVT ensemble is employed with temperature coupled to a Berendsen thermostat34 with a coupling constant of 4 ps. The second-order stochastic dynamic integrator in GROMACS is used, and all simulations are run with a time step of 40 fs. For simulations with explicit water, the default simulation parameters specific to Martini and SDK CG force fields are used, including the potential functions and cutoff schemes, which differ between the two FFs. For both FFs, an NPT ensemble is used to maintain both pressure and temperature. Other parameters, such as temperature and simulation box size, are kept the same as those used with the implicit water model. The time step used for simulations with explicit water is 10 fs with a leapfrog integration algorithm. The following describes the specific simulations carried out. 2.1. Simulations of Single PEO Chains. Single PEO chains with 10−600 EO beads, i.e., PEO10 to PEO600, are simulated for 1 μs each at 296 K, with the last 0.6 μs used for analysis. Boxes of sizes ranging from 12 nm × 12 nm × 12 nm to 60 nm × 60 nm × 60 nm are used for different chain lengths to make sure there is no interaction between PEO chains through the periodic boundary conditions. 2.2. Simulations of Multiple Chains of PEO77. From 1 to 72 chains of PEO77 each with 77 CG EO beads, corresponding to 2−148 mg/cm3, are randomly placed in a simulation box of 14 nm × 14 nm × 14 nm. A 1 μs of simulation is performed at 296 K for each concentration with the last 0.6 μs used for analysis. 2.3. Potential of Mean Force for Pulling a Short PEO Chain from a Hydrocarbon Slab. We use umbrella sampling35 to calculate the potential of mean force (PMF) for pulling a PEO8 molecule from a dodecane slab at 300 K. The procedure for calculating the PMF has been described elsewhere.22 Specifically, we generate 60−80 frames of configurations with spacing of 0.4 Å between the center of mass of the dodecane slab and the PEO chain. After a short preequilibration, a restrained molecular dynamics simulation is performed for 100 ns for each frame with a force constant of 1000 kJ/mol/nm2 applied to the center of mass of the PEO chain. We used g_wham, i.e., the weighted histogram analysis method (WHAM) implementation in GROMACS,36 to calculate the potential of mean force. 2.4. Adsorption of PEOs onto Hydrocarbon Surfaces. 2.4.1. Planar Surface. A planar hydrocarbon slab consisting of 300 hexadecanes (with each molecule represented by four bonded C1 beads) is first simulated for 100 ns in a box of size 6.0 nm × 6.0 nm × 9.0 nm. After equilibrium, either 1 or 20 copies of PEO20 or PEO40 are added to the continuum (i.e., water phase). An additional 400 ns simulation is then performed at 300 K. The last 200 ns, by which time the number of adsorbed EO beads has reached equilibrium, are used to calculate the equilibrium number of EO beads adsorbed onto the hexadecane slab. An EO bead is considered adsorbed if its distance to any C1 bead in hexadecane is less than 0.7 nm. 2.4.2. Adsorption of PEOs onto a Droplet. These droplet simulations are very similar to those for the planar surface. A

Table 1. Tabulated Parameters for Nonbonded (LennardJones) and Bonded (Harmonic) Potentials for Different Beads Used in This Studya interaction potentials LJ C1 C1 EO EO OA OA C1 EO C1 OA EO OA bond C1−C1 EO−EO C1−EO EO−OA angle C1−C1−C1 EO−EO−EO C1−C1−EO C1−EO−EO EO−EO−OA distance restraints EOi EOi+3

σ (nm) 0.47 0.43 0.43 0.47 0.47 0.43 d (nm) 0.480 0.330 0.390 0.280 θ (deg) 180 150 180 180 150 cutoff (r1, r2) r1 = 0.00 r2 = 0.60

ε (kJ/mol) 4.50 0.20 0.20 1.90 1.50 0.20 kb (kJ/mol/nm2) 1250.0 16000.0 3000.0 22000.0 kθ (kJ/mol/rad2) 35.0 80.0 20.0 20.0 80.0 Kdr (kJ/mol/nm2) 60

a

The parameters for C1 beads are borrowed from Dry Martini, and the rest of the parameters are derived in this study.

Martini FF, although use of the theta-type potential yields nearly identical results to those from the cosine-type potential. Additionally, in standard Martini FF a dihedral potential with four individual terms is used to get reasonable end-to-end distances and radii of gyration for PEO.10 However, use of a dihedral limits the time step to no larger than 10 fs. The Ryckaert−Bellemans dihedral function was proposed by Rossi et al.11 in a later study with smaller equilibrium values for bond lengths and bond angles to improve further the numerical instability of simulations, which thereby allows a time step of 20 fs. Using our implicit water model, we could obtain very good agreement with regular Martini results, such as the probability distribution of bond lengths, angles, and dihedrals, the radii of gyration, and end-to-end distance for PEO of different lengths, by using the same set of parameters for the dihedral potentials as originally put forward by Lee et al.10 for the standard Martini FF. However, as in the standard Martini FF, with these dihedral potentials the implicit water model can only use a small time step of 10 fs. To prevent numerical instabilities when using a larger time step, instead of using a dihedral potential function, we here impose piecewise linear/harmonic distance restraints (bond type 10 in GROMACS) on the separation of the ith and i+3rd EO beads in a PEO chain, as shown below: ⎧0 ri , i + 3 < r1 ⎪ ⎪ r1 ≤ ri , i + 3 < r2 ⎪ 0.5Kdr(ri , i + 3 − r1)2 V (ri , i + 3) = ⎨ ⎪ 0.5K (r − r )(2r dr 2 1 i , i + 3 − r2 r2 ≤ ri , i + 3 ⎪ ⎪ − r) ⎩ 1

Using these distance restraints, we obtain very good agreement with the results from the standard Martini FF as discussed below. More importantly, we can use a time step of 40 fs by using distance restraints rather than a dihedral potential. Please 7711

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules

water beads. For example, using 12 cores of 2.60 GHz Intel Xeon E5-2670 processors, we could obtain ∼4.8 μs/day for 150 C12E45C12 with box sizes ranging from 14 to 20 nm and ∼1.8 μs/day for 150 C12E136C12 with box sizes ranging from 25 to 50 nm. As discussed by ref 28, the speed-up is usually system dependent, and other issues, such as homogeneity of the system and the domain decomposition for parallel computation (due to the presence of vacuum, i.e., water phase), also affect the efficiency of the implicit solvent model. In addition, due to coarse-graining and implicit water, the dynamic exchange of hydrophobes among clusters also happens fast. Simulations usually reach equilibrium in less than 1.0 μs. For example, for 150 C12E45C12 in a box of 13 nm × 13 nm × 13 nm, we observe exchange of hydrophobes among different clusters in less than 0.5 μs (see Figure 6). Systems with long chains reach equilibrium faster due to smaller hydrophobe cluster sizes.

spherical droplet of 60 hexadecane molecules is simulated for 100 ns in a box of size 6.6 nm × 6.6 nm × 6.6 nm. After equilibrium, either 1 or 20 copies of PEO20 or PEO40 are added to the continuum (i.e., water phase). An additional 400 ns simulation is then performed at 300 K, and the last 200 ns is used to calculate the equilibrium number of EO beads adsorbed using the same 0.7 nm cut off distance stated above. Note that we also perform simulations of the planar surface and droplet scenarios using SDK and Martini FFs. The procedure in each case is similar to that described above except that hydrocarbon slab or droplet is simulated in explicit water initially before PEO20 or PEO40 are added to the water phase. 2.5. C12E5 Micelle Simulations. Pentaethylene glycol monododecyl ether (C12E5) is modeled using our implicit water model using the standard Dry Martini C1 bead to represent four carbons, an EO bead for each ethylene monomer, and an OA bead to represent the terminal hydroxymethyl group. As shown in Table 1, the OA group is very similar to an EO bead. These parameters were obtained by comparing the results from implicit water model with those from the atomistic GROMOS 53a6OXY+D force field.9 A preassembled micelle with an aggregation number of 40, 54, or 80 is simulated in a box of 9 nm × 9 nm × 9 nm at 300 K. A 150 ns simulation is performed for each case, and the radial distribution functions, g(r), of C1 and EO beads are calculated using the last 100 ns. 2.6. Simulation of Self-Assembly of Surfactants with Long EO Blocks. The self-assembly polyoxyethylene alkyl ether with a long EO block, C12E100, is studied at 300 K using the implicit water model. At the beginning of the simulation, 100 chains are placed randomly in boxes of sizes ranging from 14 nm × 14 nm × 14 nm to 40 nm × 40 nm × 40 nm to obtain a range of concentrations of interest. A simulation of 4 μs duration is performed for each concentration, and the last 2 μs is used to calculate the cluster size distribution based on a cutoff distance of 0.7 nm between C1 beads: if any C1 bead in a hydrophobe is closer to a cluster than the cutoff distance of 0.7 nm, it is included in the cluster. 2.7. Self-Assembly and Phase Separation of Telechelic Polymers. We study the self-assembly of 100 telechelic molecules C12E100C12 and C12E200C12 and 150 telechelic molecules C12E45C12, C12E136C12, C12E227C12, and C12E300C12 at 300 K. To maintain the same concentration, the box sizes of C12E100C12 and C12E200C12 are chosen such that the concentrations of the telechelic molecules are kept the same as that for the simulations of C12E100 discussed in section 2.6. The box sizes for C12E45C12 and C12E136C12 are varied from 12 nm × 12 nm × 12 nm to 50 nm × 50 nm × 50 nm to achieve the concentrations where a phase transition can be observed. For comparison, we also simulate C12E227C12 with box sizes of 60 nm × 60 nm × 60 and 80 nm × 80 nm × 80 nm and C12E300C12 with box sizes of 70 nm × 70 nm × 70 and 90 nm × 90 nm × 90 nm. These polymer chain lengths are selected in order to compare directly with experimental results.37 Starting from a random configuration, a 4 μs simulation is performed for each case at 300 K, and the last 2 μs is used to calculate the cluster size distribution or to analyze phase behavior. 2.8. Simulation Speed. The main advantage of implicit solvent model is decrease of computational cost due to a reduced number of beads. In our simulations, especially for long chain molecules, the improvement over explicit solvent model is rather significant due to reduction of large amount of

3. RESULTS AND DISCUSSION 3.1. Parametrization of EO−EO Interactions. Shown in Figure 1 are the comparisons of probability distributions of bond lengths, angles, and dihedrals between this study using implicit water and a previous study by Lee et al.10 using both coarse-grained (CG) Martini FF and an atomistic CHARMM FF with explicit water. As seen, we can get bond and angle distributions nearly identical to those of the standard Martini FF. The probability distribution of the dihedrals are all very close to each other around 180°, but differences are seen around 0°. This is mainly because in this study we adopt piecewise linear/harmonic distance restraints on ith and i+3rd of EO beads so that we could use a larger time step of 40 fs and still maintain numerical stability (as discussed earlier). Such a potential function is not as restrictive as the dihedral potential functions when the CG beads are close to each other, i.e., when dihedrals are small. Lastly, note that if we use in our implicit water model the same dihedral potentials as used in the standard Martini FF of Lee et al.10 (along with a time step of 10 fs to maintain stability), the probability distributions of dihedrals also become nearly identical. The radius of gyration (Rg) of a single PEO chain as a function of its molecular weight (Mw) is presented in Figure 2 for the implicit water model and compared with that from the Lee et al.10 Martini FF with explicit water. Both show a powerlaw dependence of Rg on Mw: Rg ∼ Mwν with a least-squares exponent ν = 0.55 for the latter and ν = 0.61 for the implicit model, both of which are close to the experimental value of 0.583.15 For the implicit water model, we find that the nonbonded interaction parameter ε describing the interaction between two EO beads strongly influences the exponent ν. For example, using ε = 0.5 kJ/mol in the implicit water model, we could obtain nearly identical Rg and distribution of end-to-end distances as obtained using the standard Martini FF. However, ε = 0.2 kJ/mol is a better choice to eliminate some problems associated with stronger EO interactions, such as unrealistic clustering of short polyoxyethylene alkyl ether CnEm surfactants and penetration of EO beads into the hydrophobic cores of C12E5 micelles. Table 2 gives Rg for PEO77 at different concentrations. At low concentration, Rg is 19.6 Å, a value very close to 19.1 ± 1 Å measured by small-angle neutron scattering using Guinier analysis.16 Increasing the concentration of PEO77 results in only a very minor decrease in Rg, a similar trend observed by Lee et al.10 using the standard Martini FF with explicit water. In fact, with 72 PEO77 molecules and a concentration of 148 mg/ 7712

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules

Figure 2. Plot of log Rg versus log Mw for PEO chains containing 10− 600 EO beads using the implicit water model. The results by Lee et al.10 using the coarse-grained Martini FF with explicit water are also shown for comparison. Fittings of these data using straight lines and their corresponding slopes (s) are shown.

Table 2. Average Radii of Gyration for PEO77 for Different Concentrations Calculated Using the Implicit Water Modela implicit water model

Lee et al.

concn (mg/cm3)

no. of PEO77

2.1 16.4 49.3 74.0 115.1 148.0 148.0

1 8 24 36 56 72 72

Rg (Å) 19.6 19.6 19.5 19.4 19.0 18.8 18.6

± ± ± ± ± ± ±

3.6 1.2 0.7 0.6 0.5 0.4 0.1

a

A datum obtained for a simulation of 72 chains of PEO77 from Lee et al.10 using the Martini FF with explicit water is also shown for comparison.

Figure 1. Probability distribution of (a) bond lengths, (b) angles, and (c) dihedrals of PEO37 from our implicit water model and from the results of Lee et al.10

cm3, the Rg of 18.8 ± 0.4 Å obtained using the implicit water model is almost the same as Rg of 18.6 ± 0.1 Å calculated by Lee et al.10 using the Martini FF with explicit water. These results from both Martini and our implicit water model indicate that PEO chains are coiled polymers with water acting as a good solvent and that both intra- and interchain interactions are weak. 3.2. Parameterization of EO−C1 Interactions. To parametrize the EO−C1 nonbonded interactions, we calculate the potential of mean force (PMF) for pulling a PEO8 out of a dodecane slab. The PMF as a function of distance of the center of mass of a PEO from the surface of the dodecane slab from the implicit water model is compared with those from the SDK and Martini FFs in Figure 3. The zero on the x-axis corresponds to the average surface position of the dodecane slab determined using the inflection point of the hydrocarbon density; negative x values are inside the slab and positive ones

Figure 3. Potential of mean force (PMF) as a function of distance of the center of mass of PEO8 from the surface of a dodecane slab calculated using different force fields. The ε values given are for the interactions between the C1 and EO beads.

outside. It is seen that with increasing distance along the x-axis, the PMF decreases, reaches a minimum when the PEO slightly penetrates the dodecane, and then increases to a plateau when the PEO chain has left the slab. The overall change of the PMF for the standard Martini FF is around ∼3.5kT, which is less than half that of the SDK FF, which is around ∼8.0kT. This discrepancy is mainly due to the weaker interaction between EO and C1 beads in the standard Martini FF. As shown previously by an atomistic simulation,38 PEO adsorbs strongly onto sodium dodecyl sulfate micelles, mainly due to a strong 7713

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules interaction between PEO and hydrocarbon. A recent study39 using GROMOS 53a6OXY+D reported an even higher PMF value in excess of 10kT for PEG8 pulled from a dodecane surface. These results, together with results from following adsorption study, suggest that the interaction between C1 and EO beads is weaker in the standard Martini FF, and the SDK FF is more appropriate than the Martini FF with explicit water. In Figure 3, we show the PMF calculated using the implicit water model with ε = 1.90 kJ/mol and ε = 1.80 kJ/mol for the LJ interactions between C1 and EO beads. ε = 1.80 kJ/mol gives slightly better overall agreement with that from the SDK FF and ε = 1.90 kJ/mol gives a slightly deeper PMF profile. In this study, we selected ε = 1.90 kJ/mol for two reasons: (1) Previous atomistic simulations38,39 found a C1 and EO interaction even stronger than that of the SDK FF as discussed above. (2) As shown below from adsorption studies, the amount of EO adsorbed onto a hydrocarbon slab or droplet using ε = 1.90 kJ/mol in the implicit water model is slightly smaller than that calculated using the SDK FF due to the neglect of water/PEO interactions in the former. Please note that a further increase of ε, for example, to 2.10 kJ/mol, results in unrealistic penetration of EO beads into micelle hydrophobic core of CnEm micelles or into the hydrocarbon slab, as discussed below. To further validate the EO and C1 parameters, we show in Table 3 the equilibrium adsorption of PEO chains onto either a

does the Martini FF, we compare the structural properties of C12E5 micelles of different aggregation numbers calculated using the implicit water model with those from the SDK FF in Table 4. For SDK, Rg is 17.0 ± 0.3 Å for a micelle with an Table 4. Structural Properties of C12E5 Micelles force fields

aggregation no.

implicit water

40 57 80 40 80

SDK

Rg (Å)

I1/I3

I2/I3

I1/I2

± ± ± ± ±

1.27 1.27 1.31 1.23 1.20

1.16 1.16 1.19 1.13 1.11

1.10 1.09 1.01 1.09 1.08

16.8 18.4 20.8 17.0 20.8

0.4 0.3 0.3 0.3 0.2

aggregation number N of 40 and 20.8 ± 0.2 Å for N = 80. Using the implicit water model, we obtain very similar results: 16.8 ± 0.4 Å for N = 40 and 20.8 ± 0.3 Å for N = 80. The ratios of the three moments of inertia I1, I2, and I3 evaluated along each of the principal axes of the radius of gyration tensor, I1/ I3:I2/I3:I1/I2, for different aggregation numbers are also very similar for SDK and the implicit water model as shown in Table 4. We note that results from both force fields are also close to the results obtained by Velinova et al.17 using the standard Martini FF. This means that the interaction of C1 with EO does not strongly influence the C12E5 micelle size and shape. This is understandable because the EO beads are crowded together in the C12E5 micelle corona, where their interactions with the hydrocarbon core can have little influence on the micelle structure. Figure 4 shows the radial distribution function, g(r), for hydrophobic C beads and EO beads as a function of distance

Table 3. Number of EO Beads of PEO40 and PEO20 Chains Adsorbed onto a Slab of 300 Hexadecane Molecules or onto a Droplet of 60 Hexadecane Molecules for the Standard Martini FF, SDK FF, and the Implicit Water Modela no. of adsorbed EO beads scenarios slab droplet

cases C300P40N1 C300P20N1 C60P40N1 C60P20N12

Martini FF 27.5 14.3 16.4 96.4

± ± ± ±

4.9 3.2 7.4 16.7

SDK FF 39.2 19.7 38.6 192.6

± ± ± ±

1.3 0.9 2.1 1.5

implicit 37.4 18.6 36.0 185.6

± ± ± ±

2.2 1.6 2.0 2.6

The simulations are designated by nomenclature “CxPyNz”, where x is the number hexadecane chains, y is the number of EO beads in the PEO chain, and z is the total number of PEO chains.

a

planar slab or a spherical droplet of hexadecane for the standard Martini FF, the SDK FF, and the implicit water model. As is seen, the number of EOs adsorbed onto either the hydrocarbon slab or onto the droplet for the standard Martini FF is consistently the smallest at both low and high PEO concentrations. For the droplet, for example, only 16.4 EO beads of PEO40 on average are adsorbed for the standard Martini FF compared to 38.6 EO beads adsorbed using the SDK FF and 36 using the implicit water model. For 12 PEO20 molecules, which are enough to nearly achieve saturation by EO beads, only 96.4 EO beads are adsorbed for the standard Martini while 192.6 EO beads are adsorbed for the SDK FF and 185.6 for the implicit water model. Note that the adsorption is slightly smaller for the implicit model than for the SDK FF because of the neglect of water/PEO interactions in the former. We could not further increase ε for EO and C1 beads without increasing the solubility of PEO into hydrocarbon, as noted earlier. 3.3. Applications of Implicit Water Models. 3.3.1. CnEm Micelles. Because the SDK FF shows a more realistic result for the EO interaction with hydrocarbon than

Figure 4. Radial distribution function, g(r), for C12E5 micelles with aggregation numbers N of 80 and 40. C beads and EO beads are shown for both cases. Lines represent results from the implicit water model and the data points are for SDK.

from micelle center of mass (COM) for C12E5 micelles. As seen, g(r) calculated using the implicit solvent model is almost the same as those obtained using SDK FF for micelles with n = 40 and 80. These functions indicate that the C12E5 micelle hydrophobic core is composed of exclusively C beads and hydrophilic EO beads forms the micelle corona. There is little penetration of EO beads into the micelle core even when there is strong adsorption of EO onto hydrocarbons as discussed above. We also simulate the self-assembly of surfactant with a much longer EO block, i.e., C12E100, and in Figure 5a give the micelle size distribution for three surfactant concentrations 7714

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules

3.3.3. Self-Assembly and Phase Diagram of Telechelic Polymers. Telechelic PEO chains that are capped at both ends by alkyl groups form flower-like micelles in solution where under dilute conditions both hydrophobic terminal groups on any one chain are in the same aggregate. At high concentrations, the two hydrophobic groups on a chain can bridge between two different micelles. The self-assembly of telechelic polymers have been studied by self-consistent field theory42,43 and other numerical methods,44 but molecular dynamics simulations of flower-like micelle is rare;45 to the best of our knowledge, there have been no previous simulations of the self-assembly of long chain telechelic polymers. Using the implicit solvent model, we here simulate the selfassembly of C12E100C12 and C12E200C12, both of which are model HEUR molecules. We calculate the cluster size distribution by defining a “cluster” as an aggregate of terminal hydrophobes (excluding the attached PEO chains). The typical cluster size distribution at a given concentration is similar to that of C12E100 shown above; i.e., there is a large number of free hydrophobes, a few small clusters of hydrophobes, and a large number of clusters of size around the most probable value. The latter can also be described by a Gaussian function. The average cluster sizes as functions of concentration for both C12E100C12 and C12E200C12 are compared with C12E100 in Figure 5b. For C12E200C12, the average cluster size is similar to that of C12E100, even though the former contain bridged flower-like micelles while the latter form only regular micelles. The average cluster size increases from 13 to 22.5 with increasing of polymer concentration, a trend similar to that seen in experiments: using fluorescent probe techniques, Elliot et al.46 found that the mean aggregation number for hydrophobically modified ethoxylated urethanes (HEUR) with C12 and EO195 coupled through 4,4-methylene bis(dicyclohexyl)isocynate units (in the presence of 0.1 mM SDS) increases with HEUR concentration and reaches a plateau value of around 27 at high concentrations. However, the concentration required to reach the plateau value is only around 2 mM (∼13.4 g/L), almost 10 times smaller than the lowest value used in our implicit water model. Using self-consistent field modeling, an aggregation number of around 20 was determined for C12E200C12 by Sprakel et al.42 At high concentration, around 200 g/L, the average cluster size for C12E100C12 is around ∼23, which is almost the same as that for C12E100 or C12E200C12 (see Figure 5b). This indicates that the average cluster size is barely affected by length of EO block or by the functionality of terminal ends when the EO length is long or at least around this EO length. Both observations are in agreement with previous studies.42,44 For example, Sprakel et al.42 observed no large differences between the aggregation numbers of the monofunctional chains and the telechelic chains and the aggregation number decreases only slightly with increase of EO length for surfactants and telechelic polymers with C12 alkyls and long EO blocks. Meng et al.44 also pointed out that the aggregation number of telechelic polymers is less sensitive to PEO molecular weight when the PEO block is long and screens the hydrophobes effectively from solvent. As concentration decreases, the cluster size decreases only slightly for C12E100C12, so that at lowest concentration we considered, around 30 g/L, the cluster size for C12E100C12 is still around 20, much larger than that for C12E100 or C12E200C12. The main reason for the large cluster size in C12E100C12 is the occurrence of a dense phase of interconnected flower-like micelles. As discussed below,

Figure 5. (a) Micelle size distributions for different concentrations (controlled by using different box sizes L) for 100 C12E100 using the implicit water model. Dashed lines correspond to the fits by Gaussian functions. The y-axis is truncated at 1.0, cutting of the large peak at singlet surfactants, to make clearer the distribution of larger aggregates. (b) Most probable cluster size as a function of concentration for C12E100, C12E100C12, and C12E200C12. Circle indicates phase separation.

above the critical micelle concentration. Each distribution features a relative large number of isolated free molecules (with peak cut off in the figure), a few small aggregates of size 2 or 3, and a relative large number of micelles distributed around the most probable size. The distribution can be fitted with a Gaussian function, as shown by the dashed lines in Figure 5a, from which we obtain the average micelle size and standard deviation. An increase in surfactant concentration results in a decrease in the number of free molecules and small clusters and an increase in the average micelle size. As shown in Figure 5b, the average micelle size increases gradually with surfactant concentration until it reaches a plateau value of around 23. The distributions of micelle sizes in Figure 5a are nearly the same for the two highest concentrations. Renou et al.40 determined the aggregation number of PEO micelles formed by hydrophobically end-capped poly(ethylene oxide) for chains with 100 EO monomers and between 12 and 22 carbons in the alkyl group. Using light scattering, they showed that the micelle aggregation number increases approximately linearly with alkyl end group and for C12E100 the aggregation number is approximately 18 by extrapolation of their data. Laflèche et al.41 estimated that the aggregation number of a mixture of 80% monofunctionalized and 20% difunctionalized PEO (5 kg/mol per C12) is roughly 20. In addition, Sprakel et al.42 also obtained an aggregation number around 20 for C12EO100 using self-consistent field theory calculations. So the average micelle size calculated using the implicit water model in this study at high concentration agrees well with the values previously obtained in both experimental and theoretical studies. 7715

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules telechelic polymers with short EO lengths tend to phase separate into an aggregate phase, where the local concentration is much higher than the average concentration. As a result, large clusters of hydrophobes can persist down to low overall concentration. Some previous experimental studies47,37 have found that HEUR molecules with short EO blocks tend to phase separate over a range of concentrations. Using our implicit water model, we study four different HEURs, i.e., C12E45C12, C12E136C12, C12E227C12, and C12E300C12, to compare with available experimental results. We use 150 chains and vary the box sizes to get different HEUR concentrations. For C12E45C12, we observe in Figure 6 that at high concentration

Figure 7. Same as in Figure 6, except for C12E136C12.

is much smaller than for C12E45C12 mainly due to the long EO block and the smaller sizes of hydrophobe clusters. In contrast to C12E45C12 in Figure 6 or C12E136C12 in Figure 7, no obvious phase separation is seen for C12E227C12 and C12E300C12, even at a low polymer concentration ∼0.005 g/g, as seen in Figure 8. This is mainly due to the formation of

Figure 6. Simulation snapshots for C12E45C12 at various concentrations. From (a) to (e) different box sizes with different volume fractions are shown. In each snapshot, the original simulation box is outlined, and the rest of the images are periodic copies. EO beads are shown in red, and C1 beads inside the original simulation box are shown in blue. (f) is a snapshot of only C1 beads at t = 1.6 μs for a case with box size of 13 nm and concentration 27.0%. The C1 beads in different clusters are shown in different colors. The colors of the respective C1 beads are retained as the C1 beads change from one cluster to another leading to the image (g) at later time t = 2.0 μs.

the hydrophobic tails form interconnected clusters. The hydrophobes in these cluster can also exchange among different cluster as shown in Figure 6f,g, where the hydrophobes in respective clusters at t = 1.6 μs are redistributed to different micelles at t = 2.0 μs. With decreasing surfactant concentration, regions with no C12E45C12 appear and increase in size. At low concentration, it seems clear that a concentrated and a dilute phase have appeared, although in the small boxes we use for our simulations, it is difficult to determine the volume fraction of each phase. If we assume, however, that the concentration of the C12E45C12 in its separated phase corresponds to the boundary of the two-phase region on a phase diagram, phase separation for C12E45C12 can be estimated to occur roughly between the box sizes of 13 and 14 nm at a concentration of 0.206−0.257 g/g. Please note for direct comparison with experimental results we calculate the HEUR density in unit of g/g by assuming bulk density of 1 g/cm3. For C12E136C12, a similar trend is observed, as seen in Figure 7. Decreasing the C12E136C12 concentration leads to separation into two different phases at box sizes between 27 and 30 nm or a concentration between 0.058 and 0.079 g/g. The critical concentration for phase separation of C12E136C12

Figure 8. Same as Figure 6, except for (a, b) C12E227C12 and (c, d) C12E300C12.

small hydrophobe clusters for telechelic polymer with long EOs as discussed above and facile dissociation and association of hydrophobes with these clusters. We find that the critical concentrations where phase separation occurs for both C12E45C12 and C12E136C12 determined from the implicit water model are in good agreement with experimental results for the same chain lengths by François et al.37 as shown in Figure 9. A similar decrease of the critical concentration upon increasing the length of EO interblock and no phase separation for HEURs with 227 EOs and longer have also been observed experimentally at 300 K.37 7716

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules

size distribution of hydrophobically end-capped poly(ethylene oxide): C12E100C12 and C12E200C12. We found that similar to C12E100, the most probable cluster size for C12E200C12 also increases with polymer concentration and reaches a hydrophobe aggregation number of around 22.5 at high concentration. Aggregation numbers for both the singlehydrophobe and two-hydrophobes end-capped molecules are similar to results from previous experimental and theoretical studies. We pointed out that the most probable cluster size for C12E100C12 decreases at a slower rate with decreasing concentration than that for C12E200C12 due to phase separation at low concentration for C12E100C12. We estimated using implicit water model the phase transition concentrations for C12E45C12 and C12E136C12 to be around 0.206−0.257 g/g and 0.058−0.079 g/g, respectively, and showed that longer chains C12E227C12 and C12E300C12 do not phase separate even at low concentrations, at 300 K, in agreement with previous experimental studies. In summary, we parametrized poly(ethylene oxide) using a coarse-grained model with implicit solvent. These newly developed parameters allowed us to investigate and explain some previous studies, such as the self-assembly (and phase behaviors) of long chain surfactants and HEURs. These parameters opened doors to investigate many other phenomena that involve long chain HEURs, such as the adsorption isotherm of HEURs on to latex particle surfaces, the formation of bridges, loops, and micelles between surfaces, or the competitive adsorption in the presence of other surfactants, and so on.

Figure 9. Comparison of the phase diagram for C12E45C12 and C12E136C12 between simulation results using implicit solvent CG model and experimental results by François et al.37 Simulation results are obtained at 27 °C or 300 K.

4. CONCLUSIONS We have developed parameters of poly(ethylene oxide) (PEO) for coarse-grained molecular dynamics simulations with implicit solvent. We showed that using these parameters, we could match the probability distributions of angles, bond lengths, and dihedrals of PEO chains from the Martini force field (FF) with explicit water. The radii of gyration (Rg) calculated using the implicit water model are close to those from the Martini FF with explicit water, and the Rg for both FFs exhibits an exponential dependence on Mw: Rg ∼ Mwν. The Martini FF yields an exponent of ν = 0.55 while the implicit water model gives ν = 0.61. We also showed that similar to results from the Martini FF, the Rg of the PEO chains from the implicit model exhibit very little dependence on PEO concentration, indicating weak interchain interactions. To model ethylene oxide chains capped at one or both ends by alkanes, we tuned the interaction parameters between the EO bead used to represent the ethylene oxide monomer and the C1 bead used to represent a four-carbon alkane chain. Using umbrella sampling to obtain potential of mean force (PMF), the tuning was done by matching PMFs from the implicit solvent model to those obtained from explicit-solvent models. We showed that with proper choice of the interaction parameters the PMF for pulling a PEO chain from a dodecane slab calculated using implicit water is close to that calculated using the SDK CG force field. We determined that further increase in the EO−C1 interaction parameter would result in unreasonable PEO penetration into a hydrocarbon slab. We verified that the adsorption of EO beads onto hydrocarbon slab or onto a small hydrocarbon droplet calculated using the implicit water model agrees with that from the SDK FF. We studied CnEm surfactant micelles using our coarsegrained implicit solvent model and obtained nearly the same sizes and monomer density profiles for C12E5 micelles as are obtained from the SDK FF, for two different aggregation numbers. In addition, the implicit solvent model allowed us to simulate directly the self-assembly and equilibration of poly(oxyethylene) alkyl nonionic surfactant with long EO blocks. Specifically, we calculated the cluster size distribution of C12E100 and found that at high concentration the most probable cluster aggregation number reaches a plateau value of 23 that is very close to that from previous experimental and theoretical studies. We also calculated the hydrophobe cluster

■

AUTHOR INFORMATION

Corresponding Author

*E-mail [email protected] (R.G.L.). Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS We are grateful to The Dow Chemical Company for funding this research and to Valeriy V. Ginzburg, Antony K. van Dyk, Tirtha Chatterjee, Alan Nakatani, and Alex Kalos of The Dow Chemical Company for valuable discussion and inspiration that guided the research. Computations were carried out on supercomputers provided by the Extreme Science and Engineering Discovery Environment (XSEDE).

■

REFERENCES

(1) Harris, J. M.; Chess, R. B. Nat. Rev. Drug Discovery 2003, 2, 214− 221. (2) Xue, Y.; O’Mara, M. L.; Surawski, P. P. T.; Trau, M.; Mark, A. E. Langmuir 2011, 27, 296−303. (3) Dong, R.; Hao, J. Chem. Rev. 2010, 110, 4978−5022. (4) Taylor, K. C.; Nasr-El-Din, H. A. J. Pet. Sci. Eng. 1998, 19, 265− 280. (5) Kästner, U. Colloids Surf., A 2001, 183−185, 805−821. (6) Israelachvili, J. Proc. Natl. Acad. Sci. U. S. A. 1997, 94, 8378−8379. (7) Dormidontova, E. E. Macromolecules 2002, 35, 987−1001. (8) Chen, C.; Depa, P.; Sakai, V. G.; Maranas, J. K.; Lynn, J. W.; Peral, I.; Copley, J. R. D. J. Chem. Phys. 2006, 124, 234901. (9) Fuchs, P. F. J.; Hansen, H. S.; Hünenberger, P. H.; Horta, B. A. C. J. Chem. Theory Comput. 2012, 8, 3943−3963. (10) Lee, H.; de Vries, A. H.; Marrink, S.-J. J.; Pastor, R. W. J. Phys. Chem. B 2009, 113, 13186−13194. (11) Rossi, G.; Fuchs, P. F. J.; Barnoud, J.; Monticelli, L. J. Phys. Chem. B 2012, 116, 14353−14362. 7717

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718

Article

Macromolecules (12) Shinoda, W.; Devane, R.; Klein, M. Mol. Simul. 2007, 33, 27−36. (13) Prasitnok, K.; Wilson, M. R. Phys. Chem. Chem. Phys. 2013, 15, 17093−17104. (14) Wang, Q.; Keffer, D. J.; Nicholson, D. M. J. Chem. Phys. 2011, 135, 214903. (15) Devanand, K.; Selser, J. C. Macromolecules 1991, 24, 5943− 5947. (16) Thiyagarajan, P.; Chaiko, D.; Hjelm, R. P. J. Macromolecules 1995, 28, 7730−7736. (17) Velinova, M.; Sengupta, D.; Tadjer, A. V.; Marrink, S.-J. Langmuir 2011, 27, 14071−14077. (18) Levine, B. G.; LeBard, D. N.; DeVane, R.; Shinoda, W.; Kohlmeyer, A.; Klein, M. L. J. Chem. Theory Comput. 2011, 7, 4135− 4145. (19) Lee, H.; Pastor, R. W. J. Phys. Chem. B 2011, 115, 7830−7837. (20) Shinoda, W.; DeVane, R.; Klein, M. L. Soft Matter 2008, 4, 2454−2462. (21) Shinoda, W.; Discher, D. E.; Klein, M. L.; Loverde, S. M. Soft Matter 2013, 9, 11549−11556. (22) Wang, S.; Larson, R. G. Langmuir 2015, 31, 1262−1271. (23) Chen, T.; Hynninen, A. P.; Prud’Homme, R. K.; Kevrekidis, I. G.; Panagiotopoulos, A. Z. J. Phys. Chem. B 2008, 112, 16357−16366. (24) Bedrov, D.; Ayyagari, C.; Smith, G. D. J. Chem. Theory Comput. 2006, 2, 598−606. (25) Fischer, J.; Paschek, D.; Geiger, A.; Sadowski, G. J. Phys. Chem. B 2008, 112, 13561−13571. (26) Juneja, A.; Numata, J.; Nilsson, L.; Knapp, E. W. J. Chem. Theory Comput. 2010, 6, 1871−1883. (27) Jusufi, A.; Sanders, S.; Klein, M. L.; Panagiotopoulos, A. Z. J. Phys. Chem. B 2011, 115, 990−1001. (28) Arnarez, C.; Uusitalo, J.; Masman, M.; Ingólfsson, H. I.; de Jong, D.; Melo, M.; Periole, X.; de Vries, A.; Marrink, S. J. Chem. Theory Comput. 2015, 11, 260−275. (29) Marrink, S. J.; Risselada, H. J.; Yefimov, S.; Tieleman, D. P.; de Vries, A. H. J. Phys. Chem. B 2007, 111, 7812−7824. (30) López, C. A.; Rzepiela, A. J.; de Vries, A. H.; Dijkhuizen, L.; Hünenberger, P. H.; Marrink, S. J. J. Chem. Theory Comput. 2009, 5, 3195−3210. (31) Monticelli, L.; Kandasamy, S. J. Chem. Theory Comput. 2008, 4, 819−834. (32) Rossi, G.; Monticelli, L.; Puisto, S. R.; Vattulainen, I.; AlaNissila, T. Soft Matter 2011, 7, 698−708. (33) Van Der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. C. J. Comput. Chem. 2005, 26, 1701−1718. (34) Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. J. Chem. Phys. 1984, 81, 3684−3690. (35) Torrie, G.; Valleau, J. J. Comput. Phys. 1977, 23, 187−199. (36) Hub, J. S.; de Groot, B. L.; van der Spoel, D. J. Chem. Theory Comput. 2010, 6, 3713−3720. (37) François, J.; Beaudoin, E.; Borisov, O. Langmuir 2003, 19, 10011−10018. (38) Shang, B. Z.; Wang, Z.; Larson, R. G. J. Phys. Chem. B 2008, 112, 2888−2900. (39) Huston, K. J.; Larson, R. G. Langmuir 2015, 31, 7503−7511. (40) Renou, F.; Nicolai, T.; Nicol, E.; Benyahia, L. Langmuir 2009, 25, 515−521. (41) Laflèche, F.; Nicolai, T.; Durand, D.; Gnanou, Y.; Taton, D. Macromolecules 2003, 36, 1341−1348. (42) Sprakel, J.; Besseling, N. A. M.; Leermakers, F. A. M.; Cohen Stuart, M. A. J. Phys. Chem. B 2007, 111, 2903−2909. (43) Sprakel, J.; Besseling, N. A. M.; Cohen Stuart, M. A.; Leermakers, F. A. M. Eur. Phys. J. E: Soft Matter Biol. Phys. 2008, 25, 163−173. (44) Meng, X. X.; Russel, W. B. Macromolecules 2005, 38, 593−600. (45) Kim, M.; Choi, Y. W.; Sim, J. H.; Choo, J.; Sohn, D. J. Phys. Chem. B 2004, 108, 8269−8277. (46) Elliott, P. T.; Xing, L. L.; Wetzel, W. H.; Glass, J. E. Macromolecules 2003, 36, 8449−8460.

(47) Pham, Q. T.; Russel, W. B.; Thibeault, J. C.; Lau, W. Macromolecules 1999, 32, 2996−3005.

7718

DOI: 10.1021/acs.macromol.5b01587 Macromolecules 2015, 48, 7709−7718