Simulations Study of Single-Component and Mixed n-Alkyl-PEG Micelles.

Here, we study one-component and mixed n-alkyl-poly(ethylene glycol) (C mE n) micelles with varying poly(ethylene glycol) (PEG) chain lengths n using coarse-grained molecular simulations. These nonionic alkyl-PEG surfactants and their aggregates are widely used in bio and chemical technology. As expected, the simulations show that increasing the PEG chain length decreases the alkyl-PEG micelle core diameter and the aggregation number but also enhances PEG chain penetration to the core region and spreads the micelle corona. Both the core and corona density are heavily dependent on the PEG chain length and decrease with increasing PEG length. Furthermore, we find that the alkyl-PEG surfactants exhibit two distinct micellization modes: surfactants with short PEG chains as their hydrophilic heads aggregate with the PEG heads relatively extended. Their aggregation number and the PEG corona density are dictated by the core carbon density. For longer PEG chains, the PEG sterics, that is, the volume occupied by the PEG head group, becomes the critical factor limiting the aggregation. Finally, simulations of binary mixtures of alkyl-PEGs of two different PEG chain lengths show that even in the absence of core-freezing, the surfactants prefer the aggregate size of their single-component solutions with the segregation propelled via enthalpic contributions. The findings, especially as they provide a handle on the density and the density profile of the aggregates, raise attention to effective packing shape as a design factor of micellar systems, for example, drug transport, solubilization, or partitioning.

Introduction

Nonionic surfactants, such as n-alkyl-poly(ethylene glycols) CE, where m is the number of carbon atoms in the alkyl chain and n is the number of oxyethylene units in the alkyl ethoxylate surfactant, have a wide array of applications in detergents, cosmetics, and in industrial processes, especially in agriculture, textile, paper, and oil industries.[1] Also, their macromolecular equivalents, block copolymers containing poly(ethylene oxide), have received significant attention, see refs,[2−5] especially because of both their relatively simple chemical structure and biocompatibility that permit their use as drug carriers.[6,7]

Linear alkyl ethoxylates, CE, are especially interesting poly(ethylene glycol) (PEG)-based nonionic surfactants because their aggregation characteristics, in particular micellization, are easily tunable via modification of the alkyl and ethoxylate chain lengths. A further handle on alkyl ethoxylate micellar solutions is provided by the lower critical solution transition (LCST) response that these surfactants exhibit.[8,9] The phase separation above LCST is mainly driven by the increased entropy of solvent water molecules coupled with the weak attractive van der Waals interaction between polymer chains. In general, decreasing the tail-to-head ratio increases the cloud point temperature, decreases micelle size, and delays the onset of any temperature-dependent micellar growth, if present. For reviews, see refs (8) and (9).

Linear alkyl ethoxylates with longer alkyl chains typically form micellar, relatively monodisperse, spherical aggregates in aqueous solution.[10] This response is propelled by the low critical micellization concentration (CMC) and large interfacial tensions of such systems.[10] However, the aggregation and aggregates can be tuned via, for example, concentration, surfactant species, temperature, and additives. Specifically, the effects of surfactant concentration,[5,11] temperature,[4,5,12] and surfactant structure[13,14] on the aggregate size and morphology have been previously studied through various scattering experiments. The aggregation number has been reported to increase linearly with the alkyl chain length and decrease exponentially with the PEG chain length.[13] On the other hand, the effect of the surfactant concentration on the aggregate size is limited to high surfactant concentrations.[5] The core-corona interface of alkyl ethoxylate micelles is known to be generally very broad.[5,13] Furthermore, the aggregate structure of alkyl ethoxylates with longer PEG chains is greatly affected by the size asymmetry of the hydrophobic and hydrophilic blocks.[5,13,14]

Alkyl-PEG micelles are also thermoresponsive. The temperature dependency in the alkyl ethoxylate response arises in part from the solvent quality of the PEG block changing with temperature: water is a good solvent for PEG at room temperature but becomes a poorer solvent when the temperature approaches the theta temperature (around 100 °C).[5] The change in solvent quality changes the PEG chain extension in the solvent which influences the intermicellar interactions. Nonionic micelles of surfactants with a large head group compared to the hydrophobic part, such as C12E8 and C8E5, are considered to present negligible aggregate growth with increasing temperature, but an increase of attraction between the micelles has been argued as the temperature approaches the cloud point, typically based on self-diffusion and scattering experiments.[15−18] However, both Glatter et al.[19] and Bernheim-Groswasser et al.[20] have argued also for thermotropic aggregate growth. On the other hand, C12E8 aggregates are thought to undergo only minor changes when approaching the cloud point.[21] Additionally, direct electron microscopy imaging of pure C12E5 and C16E6 systems has revealed that spherical and rodlike micelles, as well as, more complex branched and loop-like structures can exist in alkyl ethoxylate systems.[20,22] This survey on the existing literature shows that alkyl ethoxylates present a wide variety of aggregation responses.

Computer simulations provide complementary means to characterize the aggregates and the aggregation control factors. Existing atomistic detail computational studies of CE surfactants and related molecules are limited to rather short PEG chains that can be described with reasonable computational cost; see refs (23) and (24). However, most experimental works and applications consider much longer PEG chain lengths. Reaching longer PEG chain lengths, some studies involving coarse-grained (CG) modeling with explicit[25−27] and implicit[28] water, as well as, with continuum representations[29] exist. While the CG models can reproduce the extension of a pure PEG polymer as a function of molecular weight very well,[25] and the micellar aggregation of small alkyl-PEG polymers,[26] the aggregation response of longer alkyl-PEG polymers and mixed micellar systems have received much less attention because of the higher computational demand. Simulations of pegylated lipid micelles[30−33] have, however, demonstrated that modeling can provide useful insight into the aggregation of related molecules in bulk water and at interfaces.

Here, we target the gap in alkyl-PEG simulational studies by examining medium chain length alkyl-PEG surfactants via C18E10, C18E20, and C18E50 surfactant aggregation response in single-component surfactant solutions and mixed surfactant systems containing C18E10 and C18E50 surfactants. We characterize the PEG chain length dependencies of the aggregation response, as well as, mixed micelle structure and composition presenting, to our knowledge, the first chemically specific, computational block-copolymer micellization study in which the polymeric head groups are long enough to exhibit both hydrophobic block-limited and hydrophilic block-limited micellization regimes. This enables identifying two distinct modes of micellization that result in different micellar density profiles, as well as, that the surfactant species segregation, which has been experimentally reported for mixed alkyl-PEG micelles, bears an enthalpic component that contributes toward surfactant species segregation also in the absence of core-freezing.

Computational Methods

The MARTINI CG[34,35] model within the GROMACS 4.6.5 simulation package[36,37] was employed for the alkyl-PEG micelle simulations of this work. Specifically, the MARTINI PEG model developed by Rossi et al.[26] combined together with the standard MARTINI v2.1 lipid alkyl chain beads and standard MARTINI water beads were used. The MARTINI CG model was chosen for this work, as in comparison to fully atomistic models, the CG model allows a larger system and longer time-scale to be studied with the same computational effort. These features are important for the feasibility of examining the CE micelle equilibrium structures because the studied systems need to be large enough to contain multiple micelles and the aggregation times are exponential. The speed-up is mainly due to the CG approach smoothing out the interaction potential which enables the use of a longer time-step. Additionally, the model has been successfully used to describe alkyl-PEGs previously.[26,38]

The MARTINI CG model is based on a four-to-one mapping.[34,35] This means four heavy atoms and their associated hydrogens are represented by a single CG bead in the model. This mapping is extended to the solvent so that four water molecules are effectively mapped by one CG solvent bead. The model has four different types of interaction sites (charged, polar, nonpolar, and apolar), and each site can be categorized into subtypes. The subtypes are based on, for example, polarity or hydrogen-bonding capabilities of the described group of atoms. This enables retaining the chemical nature of the molecules despite the CG representation. The MARTINI model is parametrized to reproduce the partition behavior and energies between polar and nonpolar phases.[34,35] Therefore, the model is especially suitable for modeling processes driven by, for example, hydrophobicity which includes micellization.

Micelle formation in systems containing single surfactant species was investigated by examining the behavior of 500 C18E molecules in an aqueous solution (Figure ). Three different surfactants with a fixed alkyl chain length but varying the PEG chain length were studied: C18E10, C18E20, and C18E50. Mixed surfactant systems were modeled by a random binary mixture of C18E10 and C18E50 molecules at a 1:1 mixing ratio (250 molecules of each species in all simulations). The C18E molecules were first placed randomly into the periodic simulations box and then solvated by water. The simulation systems were constructed so that the C18E w/v concentration was 0.1 g/mL for all the systems after equilibration. The concentration is well-above the observed CMCs of C18E type surfactants in the examined temperature range.[39,40]

Figure 1

Examined C18E surfactants in MARTINI representation and snapshots corresponding to a t = 0 μs initial and t = 10 μs final configuration in a simulation of 500 C18E20 molecules solvated in water. The temperature is T = 296 K and water is omitted in the visualization for clarity.

The pressure of the system was set to 1.0 bar using the Parrinello–Rahman[41] pressure control with τp = 12 ps. The temperature was controlled by the stochastic velocity rescaling thermostat developed by Bussi et al.[42] with τT = 1.0 ps. The Lennard-Jones interactions used a cutoff of 1.1 nm with a shift to zero between 0.9 and 1.1 nm. A time-step of 20 fs was used. This constitutes a standard MARTINI simulations protocol to be used with the employed PEG and alkyl chain beads.

The single surfactant species simulations were run at temperatures of 275, 296, 310, and 320 K. Unless noted otherwise, the results presented in this work are based on simulations at 296 K, while the other temperatures were examined to capture the behavior of the model outside the original parametrization temperature. For the binary surfactant mixtures, temperatures of 265, 270, 275, 296, and 310 K were examined. This range of temperatures was examined because, on the one hand, increasing the temperature enhances dynamics (makes equilibration faster), and on the other hand, the entropy contribution to free energy decreases with the decreasing temperature (enthalpy contributions to, for example, ordering can be expected to show at low temperatures). As is standard for the MARTINI model, for temperatures at and below 275 K, 10% of the CG water particles W were replaced by CG antifreeze water particle WF to prevent freezing of the CG solvent. It is worth noting that the surfactant alkyl tails freeze at low temperatures in experiments,[12,43] but the MARTINI model cannot be expected to reproduce this response as it is a coarse grained model parametrized to function around biological temperatures.

Additionally, the radii of gyration of single-PEG chains of 15 monomers in length measured by an atomistic and the MARTINI CG model were compared. The atomistic detail description was given by the Gromos53a6 united atom force-field[23] within GROMACS 4.6.5 simulation package.[36,37] In the temperature testing, the simulation temperatures in the 15 monomer PEG chain simulation were 270, 296, 310, and 350 K. The CG simulations are set up as described above, and in the atomistic simulations, the cubicle contains 6349 water molecules which corresponds to a system size of (5.8 × 5.8 × 5.8) nm3.

For the atomistic simulations, the Berendsen weak pressure coupling[44] with a reference pressure of 1.0 bar and τp = 1 ps is used. The temperature is controlled by the stochastic velocity rescaling thermostat developed by Bussi et al.[42] with τT = 1.0 ps. The van der Waals interactions were truncated at 1.4 nm. A time step of 2 fs was employed. The atomistic simulations were run for a total of 500 ns. Out of these total simulation durations, for the single-chain simulations, the first 0.5 μs and 300 ns were considered as relaxation time for the CG and atomistic simulations, respectively, and omitted in the analysis in the temperature response testing.

Prior to starting the simulation runs, all systems were energy-minimized using the steepest descent algorithm. The CG micelle simulations were run for a total of 10 μs out of which the first 5 μs was considered as relaxation time and omitted in the analysis. Because of the smoothing of the interaction potential, the CG simulation time corresponds to a significantly longer actual time. For MARTINI, a factor of 4 based on, for example, water diffusion rates is typically used in converting the simulation time to real time.[34,35] The simulation times reported in this work are not converted by this factor.

Two surfactants were classified to be in the same micelle if any two of their hydrophobic tail beads (alkyl beads) were closer than 0.6 nm to each other. This classification provides a robust means of identifying aggregates based on their hydrophobic cores and enables calculating the aggregation numbers for single micelles, as well as, a size distribution of the aggregates in the system. The principal component analysis (PCA) is done with the MATLAB pca function based on the raw data. VMD was used for visualizing the molecular configurations and micelles.[45] Solvent-accessible area Asas was calculated by the algorithm developed by Eisenhaber et al.[46] available in the Gromacs simulation package tools. A probe radius rprobe = 0.263 nm was used as this corresponds closely to the radius of a MARTINI water bead.

The consistency of our modeling with previously published data was checked via the PEG chain length response of the CG model on 5–120 monomer PEG chains at 296 K. The simulation protocol followed the abovementioned description, but the simulation time was 1 μs out of which the first 0.5 μs was disregarded in the analysis. For each examined chain length, the simulation box edge length was set at dbox = 2(Rg + 2) nm = 2(N0.51 + 2) nm, and the resulting box was filled with water beads. A PEG chain length-dependent simulation box is introduced to avoid interactions of the PEG chain with its periodic images. The box size is based on the estimated radius of gyration Rg = N0.51 for a random coil of N PEG monomers in a theta solvent.

Results and Discussion

We begun by checking the consistency of our modeling with previously published data on the same model. The comparison is presented in the Supporting Information. We find the radius of gyration Rg as the function of the number of PEG monomers N scales as expected based on ref (26) and report the slope of log(Rg) versus log(N) to be 0.67 with 95% confidence interval (CI) [0.66,0.68] based on linear regression; see Figure S1 in the Supporting Information. The slope is consistent with the findings of refs (25) and (26) that use the same model. However, Lee et al.[25] observed additionally a slight transition in the PEG chain extension with increased molecular weight. They concluded that low-molecular weight PEG chains exhibit an ideal chain-like behavior in aqueous solutions, while longer PEG chains (Mw > 1630 g/mol, which corresponds to approximately 36 PEG monomers) exhibit more random coil-like behavior as the radius of gyration approaches the empirically derived relation proposed by Devanand and Selser.[25,47] In our examination, and in that of Rossi et al.,[26] no clearly visible change in behavior as a function of the PEG chain length is observed. Altogether, the results of the consistency check are in agreement with the prior works.

Next, we address the micellization of C18E surfactants and the formed micelles with different PEG block sizes. Figure shows the time evolution of the average aggregation number, and Figure shows the final simulation configurations as snapshots for C18E10, C18E20 and C18E50 surfactants at T = 296 K. The data show that initial aggregation takes place fast, during the first 0.5 μs. After this, the aggregation number continues to rise gradually. The aggregates of C18E10, which have the shortest PEG chain, reach a constant aggregation number in the simulations, but for the surfactants with the longer PEG chains, even the 10 μs simulation time is not sufficient to fully equilibrate the aggregate size. However, the aggregates stabilize in size, and very little structural evolution takes place even for the alkyl surfactants with the 50 monomer long PEG chains.

Figure 2

Time evolution of the average aggregation number Nagg of C18E10, C18E20, and C18E50 surfactant micelles at T = 296 K. At 10 μs simulation duration, the average aggregation numbers are 33.3 ± 8.1, 18.5 ± 3.3, and 12.5 ± 3.1, where the error estimate is based on the standard deviation of the aggregate sizes in the system.

Figure 3

VMD snapshots of 500 (a) C18E10, (b) C18E20, and (c) C18E50 surfactants forming micelles. The snapshots correspond to the final, 10 μs simulation configurations. PEG head group color is same for surfactants in the same micelle. The shading is to capture depth in the visualization, that is, the micelles are clearly separate from each other also in the seemingly crowded C18E50 system. The temperature T = 296 K and the concentration of C18E is 0.1 g/mL in all systems. Water is omitted in the visualizations for clarity.

Figure shows that increasing the PEG block size while keeping the alkyl part constant decreases the aggregation numbers. The final mean aggregate sizes for the C18E10, C18E20, and C18E50 surfactants are 33.3 ± 8.1, 18.5 ± 3.3, and 12.5 ± 3.1 molecules. The error estimate is based on the standard deviation of the aggregate sizes in the system. The observed response is well-known based on prior experimental[7] and simulational[30] studies, and there are several thermodynamic model-based scaling laws for block copolymers that capture the response effectively; see refs (48) and (49). Indeed, this decrease of Nagg with the increasing PEG block size is captured well by an exponential theory, and our prior work on pegylated lipids demonstrates that theories for starlike micelles provide a good match.[30] In general, the observed response can be explained by the spatial hindrance and repulsion between the elongated hydrophilic blocks, and the effect is especially pronounced for surfactants containing a short hydrophilic block.[30,50]

The micelle size distribution produced by the MARTINI model here in this work is systematically biased toward smaller aggregates than what would be expected based on experiments; see refs (5)(51), and (52). For example, Zinn et al.[51] reported the aggregation number for C18E91 micelles as 27, which is significantly larger than what our simulations predict for C18E50. Sommer et al.[5] reported an aggregation number of 30 for Brij 700 (≈C18E100) surfactants. The smaller aggregation numbers in simulations can in part be explained by the n-alkyl poly(ethylene glycol) samples being typically prepared in experiments at a much higher temperature compared to the actual measuring temperature followed by a typical equilibration time span of hours before the actual measurement to ensure full equilibration of the system. Although the aggregate size has stabilized within the 10 μs simulation duration here, the distribution still very likely represents the lower end of the micelle size distribution especially for the longer PEG chains; see refs.[53−55] The employed model is also known to underestimate the free energy difference between isolated and aggregated surfactant species[26] which also favors the formation of smaller aggregates.

The snapshots of Figure show that the surfactants with a short PEG chain micellize stronger than the surfactants with longer PEG chains. The resulting decrease of Nagg with the increasing PEG chain length results from the increase of repulsion between the longer hydrophilic PEG chains. Consequently, the surfactants with longer hydrophilic PEG chains form a larger number of aggregates with the same 500 surfactant molecule number density in the box volume. The longer chains also intertwine significantly, but the micelles remain as clearly identifiable aggregates that diffuse in the simulation box as separate aggregate units.

The observations about the effect of the PEG head length here are fully consistent with the scattering experiments of Zinn et al.[13] in terms of the effect of the PEG block on the micelle formation and aggregation number. Specifically, Zinn et al. observed a clear drop in the aggregation number as the molecular mass of the PEG block increased for C27PEG amphiphilic block copolymers until a molecular mass of 10 kg/mol of the PEG block which corresponds to approximately 225 PEG monomers. For larger PEG blocks, no further significant change in the aggregation number was observed.[13] They attributed the decrease of Nagg to a change in the entropic and enthalpic contributions as the PEG block became larger.[13] A larger PEG block also causes larger spatial hindrance for the head groups and also between the micelle corona.[13,30,31]

Figure presents the radial number density and the cumulative radial particle count of the C18E micelle core and corona regions. The data are an average over all the aggregates in the system during the last 100 ns of the 10 μs simulation time. The density figure shows two partially overlapping peaks for each micelle: the short distance peak is the micelle core and the longer distance peak is the corona. Their sum is the total density. Structural examination of the simulated micelles reveals that the PEG chains partially penetrate into the alkyl core which induces a fuzzy core–corona interface. Furthermore, the core is much more extended and tightly bound for the C18E10 micelles than for the C18E20 and C18E50, which have longer PEG chains. Proportionally, also the PEG corona extends relatively longer for the short PEGs than for the longer PEGs. However, the average end-to-end distance of the C18 alkyl chains in the micelles is 15.1 ± 0.7 Å for all the examined micelles. This means that the alkyl chains of the micelle cores do not change significantly in length even though the aggregation number and volume occupied by them changes. The chains are in a liquid-disordered phase as a full extension of the alkyl chains corresponds to a length 23.0 Å.[56] Likewise, PEG segments can be found at the same distance from the micelle core center independent of the PEG chain length or Nagg, meaning that the innermost alkyl chains to which these PEG chains are attached adopt similar positioning in all micelles.

Figure 4

(a) Radial MARTINI bead number density distribution, (b) cumulative number count (CN) of the MARTINI beads, and (c) corresponding representative C18E micelle snapshots, as well as, cartoons visualizing the occurring core-limited and PEG sterics-limited aggregation of the micelles. Temperature is 296 K. The total number density is obtained as the sum of the core and corona contributions. The data are an average over all micelles in the system and calculated over the last 100 ns of simulation time as a function of the distance from the micelle center of mass r. The normalization of the radial number density distribution is so that the volume integration over r gives the number of beads of the species in the system.

In total, these observations reflect that the mean aggregation numbers of the micelles are larger for the short PEGs than for the longer PEGs (see Figure ); when the aggregation number is high, alkyl chains pack tightly to form a larger core. Both the alkyl chain density and the PEG density are high and localized to a well-defined region. Furthermore, at high aggregation numbers, the PEG density near the tightly packed alkyl core is high. However, this tight packing does not show in the radius of gyration Rg values presented in the Supporting Information: the PEG chains in the self-assembled micelles behave very similar to the individual chains and remain well-hydrated. Furthermore, the PEG heads in the micelles maintain a conformational ensemble that is very close to that of a free, fully hydrated PEG chain even when the aggregation starts to be limited by the PEG sterics and a local densification of the PEG enforced by the PEG packing occurs. The PEG heads maintain their conformational ensemble because it is energetically more favorable to change the micelle packing than the PEG conformation ensemble and the degree of hydration.

In Figure , the radial density function corresponding to the PEG corona of micelles with the longest PEG chains shows a bend at around r = 20 Å. The significance of the data is that the micelles undergo a transition with the increase of the PEG head group size from an aggregation mode where the micellar aggregation is limited by the alkyl chain packing in the core (short PEG chain heads) to an aggregation mode where the aggregation is limited by the PEG head group sterics (long PEG chain heads). To understand the origins of the transition, it is useful to consider the packing adopted by a free, solvated PEG chain: the Rg of a free chain scales as N0.67; see the Supporting Information. For short PEG chains, this scaling with a number of monomers in the PEG chain N means that the aggregation number of the micelle is determined by alkyl chain packing in the core: the radial volume element determined by the alkyl tails in the micelle core remains large enough for a short PEG chain head at all packings. However, as the number of PEG monomers N increases, the effective volume occupied by the PEG chain increases. Eventually, the long PEG chain heads of the surfactants overlap sufficiently in the micelles so that further aggregation is no longer energetically favorable: the aggregation becomes PEG sterics-limited. Figure c shows snapshots of the micelles and cartoon representations showing the corresponding core-limited and PEG sterics-limited aggregation response.

A PEG sterics-limited aggregation number leads to a smaller core and the available volume next to the core: the PEG chain segments near the core are heavily constrained in space which leads to the local radial difference in the PEG density around r = 20 Å for the C18E50 system. The changes in density can also be considered in terms of the PEG becoming its own solvent in the aggregate. This response is enhanced by a larger Nagg. Furthermore, the shift from a tighter, symmetrical distribution in the PEG chains to a density distribution with an elongated tail in Figure is consistent with results reported for pegylated lipids with significantly longer PEG blocks.[33]

Figure also shows that the size of the hydrophobic core decreases with the increasing PEG block length. The core radius of C18E10 micelles is close to 2 nm, but for C18E50, it is around 17–18 Å. For comparison, the mean end-to-end distance of the C18 alkyl chain in our simulations is 15.1 ± 0.7 Å. The decrease of the micelle core size with the increasing PEG chain length results from the increase of steric repulsion between the longer hydrophilic PEG chains which leads to a smaller aggregation number.

Even though an underestimation of the aggregation numbers by the model could contribute to these findings, the observed decrease of the core size is in full agreement with the SAXS data by Zinn et al. for C28E micelles.[43] Their data show a similar response to the increasing PEG block length.[43] Furthermore, the observed micelle core size dimensions, as well as the mean alkyl chain length in our simulations, are in close match with the core radius of 15 Å reported for the C18PEO5 (≈C18E110) micelles by Zinn et al.[13] The core size observed here is in agreement also with the 15–17 Å core radius reported for Brij 700 (≈C18E100) by Sommer et al.[5]

Figure presents the solvent-accessible area Asas of single-PEG chains. The analysis does not capture the absolute solvent-accessible area because the model is a CG model, but the trends with respect to the PEG chain length can be expected to be visible. At first sight, the average solvent accessible-area appears to grow linearly, but actually for short PEG blocks, the response deviates from linear; see the inset. The deviation can be understood via considering the ends: each PEG chain end contributes to Asas by an amount differing from the constant addition corresponding to adding PEG monomers to the bulk chain away from the ends. This is most visible in the data of short PEG chains.

Figure 5

Solvent-accessible area Asas for coarse-grained MARTINI PEG chains of varying lengths N at T = 296 K. For comparison, the data corresponding to n-alkyl linked PEG chains C18E in micelles are also presented. The inset shows the same data plotted as Asas/N vs N to show the nonlinear response of short PEG blocks. Error bars correspond to the standard deviation.

The Asas values calculated for the PEG chains in the micelle corona and for single free PEG chains are very similar. This means that the PEG chains remain very well-hydrated, that is, accessible to the solvent, also in the micelles. This persistence of the PEG chain conformation and hydration between free PEG chains and those linked to alkyl chains in the micelles is also visible in the similarity of the Rg values of the isolated PEG chains and those attached to alkyl chains and forming micelle head groups discussed earlier. However, for Asas, the alkyl tail linker and crowding in the micelle at one end decrease the PEG block size independent contribution to Asas which shows as the slightly smaller Asas values of the alkyl chain-linked PEG chains in the micelles.

The elongated tail end of the density distribution for longer PEG chains, namely C18E50, shown in Figure , is similar to the transition in corona density observed by Sangwai and Sureshkumar,[57] as surfactant micelles undergo shape transition from spherical to rodlike. Here, the C18E micelles remain close to spherical at all examined PEG chain lengths. To show this, Figure presents the principal moments of inertia for the micelle cores and the cross sections of ellipsoids corresponding to the square root of the inertia moment values. To show the maximum deviation from the sphere, the presented cross sections correspond to the first and the third principal axis of rotation for the ellipsoid, that is, the xz-plane cross section. The data correspond to an average of the micelle cores over the last 0.5 μs of the simulation time, and the values are normalized so that for a perfect sphere, each component has a value of 1. Although the principal moments of inertia indicate some deviation from a perfectly spherical form, the ellipsoid cross sections at the xz-plane show that this deviation remains very small for all studied systems. The principal moments of inertia are also relatively independent of the aggregation number (see the Supporting Information for data), but for the surfactants with longer PEG chains, more core-shaped fluctuations occur because the micelle aggregation numbers are smaller. This is in agreement with the PEG head group size forcing the micelle to assume a looser packing because of steric hindrance of the PEG chains. Notably, for shorter CE surfactants than examined here, growth in the aggregation number has been linked with transition in the micelle shape from spherical to rodlike.[58]

Figure 6

(a) Normalized principal moments of inertia for the C18E micelle carbon cores for PEG chains of varying lengths and (b) cross sections of the ellipsoids corresponding to the square root of the normalized principal moments of inertia. The cross section corresponds to the largest, that is x, and the smallest, that is z, moments; the y moment falls between the x and z in magnitude, so the xz-cross section reveals the maximum deviation from the spherical shape. The principal moments in Cartesian coordinates are I, where i = x,y,z, and the total moment of inertia is Itot. The normalization is set as for a sphere. The error bars represent the standard deviation.

Figure presents the PEG chain length dependency of the C18E micelle corona thickness. The plot shows a log–log plot of the mean end-to-end distance d of the PEG blocks over the last 100 ns of the 10 μs simulations versus PEG chain length. The PEG chains are slightly more extended in the micelle coronas than as free, unconstrained PEG chains in aqueous solution. At a denser PEG chain packing in the head group region, significant extension of the chain would be expected (see ref (30)), but in the micelles here, the head group density is low enough to keep the extension modest.

Figure 7

Log–log presentation of the end-to-end distance d of the PEG chains as a function of the PEG chain length N at T = 296 K. The data for both the PEG chains in the alkyl-PEG micelles and as individual, free chains are presented. The error bars represent ± log-transformed error based on the standard deviation. The slope of a line fitted to the data for individual, unconstrained PEGs, is 0.6960, with R2 = 0.9971.

Mixed Micelles

The aggregation numbers in Figure and the analysis of the micelle structure for alkyl-PEGs with different PEG chain lengths show that the PEG length influences the resulting micelle significantly. To examine how this PEG chain length dependency in the aggregate structure translates to the response of a binary mixture of alkyl-PEGs of different PEG chain lengths, we examined a 1:1 mixture of C18E10 and C18E50 surfactants. Here, we examine a range of temperatures between 265 and 310 K. The temperature is lowered to reduce the entropy of mixing, but on the other hand, a higher simulation temperature ensures faster dynamics and a more complete mixing.

The work employs a CG model where the balance of entropic and enthalpic effects is shifted significantly by the coarse-graining.[35] In general, CG models tend to produce a skewed response outside their parametrization range because of this. The MARTINI CG model for CE surfactants developed by Rossi et al.[26] is constructed to accurately reproduce the phase behavior of C12E2, C12E4, and C12E6 under isothermal conditions with the goal parametrization temperature being 300 K. The T = 296 K employed as the main simulation temperature in this study is within the parametrization range, but the range of temperatures in which the mixed micelle system is examined is not.

Therefore, before proceeding to examining the mixed micellar systems at varying temperatures, Figure presents a brief summary of the basic temperature response of the employed model. Panel (a) shows that the radius of gyration of a short, free PEG chain in water remains approximately constant both in a more rigorous detail atomistic PEG model and in the MARTINI CG model as the temperature varied. The CG model systematically underestimates the Rg, but the difference is approximately one CG bead in size. In total, the MARTINI model appears to capture the hydrodynamic radius of PEG chains in water, which is in agreement with the atomistic model. We note that the comparison mapping here is done for relatively short PEG chains; longer chains could experience differences in the persistence length between the models. Panel (b) shows that the final micelle aggregation numbers also are relatively independent of temperature for the C18E20 and C18E50 surfactant systems within the examined temperature range using this model. However, the surfactant with the shortest PEG chain, C18E10, appears to form micelles that have a slightly higher aggregation number at lower temperatures (T = 275 K). This curious change in the aggregation behavior could be associated with the presence of the antifreeze particles in the system which changes the water dynamics, as well as its character as the PEG solvent; we emphasize that the change is not due to core-freezing as demonstrated by panel (c) which shows that the n-alkyl chains in the alkyl-PEGs do not exhibit any temperature dependency in the examined temperature range in the CG model.

Figure 8

(a) Radius of gyration Rg of 15 monomer long PEG chains by the atomistic and CG MARTINI models. (b) Average aggregation number Nagg of C18E10, C18E20, and C18E50 surfactant micelles at different temperatures. c) Average radius of gyration Rg and end-to-end distance d for the C18E10 carbon cores and n-alkyl chains. The data show that the alkyl chains do not experience freezing. The data corresponding to the micelles are calculated at 10 μs simulation duration.

The lack of temperature dependency in the n-alkyl data in Figure c means that the alkyl core remains in a liquid-disordered phase throughout the temperature range, that is, core crystallization does not occur at the lowest examined temperatures in our simulations. Experimental works point clearly toward core crystallization taking place; see refs (12)(43), and (59). The T = 296 K employed as the main simulation temperature in this study is within the CG model parametrization range, but it is not really surprising that in temperatures well below, the behavior is less than stellar. According to Nawaz and Carbone,[60] the poor temperature transferability of the MARTINI model is mainly due to the lack of accurate representation of nonbonded interactions, such as hydrogen bonds. This means that the examination of the effect temperature change from the parametrization temperature here is meaningful only as a means to reduce the entropy of mixing (lower temperature) or introduce faster dynamics and a more complete mixing (higher temperature); the model cannot be relied on reproducing temperature-dependent structural changes.

Figure presents the ratio of C18E10 and C18E50 surfactants in individual micelles of different aggregate sizes at 265, 296, and 310 K temperatures along with an ellipse describing the two principal components of the data distribution, as obtained from PCA: the plotted ellipse is based on the eigenvectors of the covariance matrix scaled by the square root of the corresponding eigenvalue. The 296 K data set is presented in match with the preceding single alkyl-PEG component micellization studies. The data sets corresponding to a whole range of temperatures between 265 K and up to 310 K are shown in the Supporting Information. On the basis of Figure , a nonuniform distribution of the alkyl-PEGs into micelles of different sizes depending on their PEG chain length can be expected. Indeed, the data in Figure show signs of this but also of significant scatter in both the alkyl-PEG component distribution in the micelles and in the aggregate size. This is quite expected, as even for the micelles that contain only one type of alkyl-PEG, the micelle size distribution is relatively wide. The systems contain 500 surfactants, which indicates that tens of aggregates form. Despite the scatter of the data points, the PCA analysis reveals already for the 296 K micelles a weak bias of the alkyl-PEGs with the shorter PEG chains to preferentially reside in the larger aggregates. This bias is reflected by the tilt of the major axis of the PCA ellipse drawn on the figure, whereas the ellipse width describes the uncertainty (scatter of the data). The CI of the PCA analysis is 95%. At low temperatures, where entropy of mixing contribution is reduced, this bias of alkyl-PEGs with the shorter PEG segment preferentially being in the larger aggregates shows more clearly; see Figure .

Figure 9

Ratio of C18E50 and C18E10 surfactants in micelles of binary surfactant systems simulated at 265, 296, and 310 K for 10 μs. The ellipses correspond to 95% confidence based on PCA of the data.

Even though the simulations indicate that a bias for the surfactants to reside in micelles matching their native aggregation shape and size, there is no clear phase separation of the two surfactants even at very low temperatures. This is contrary to what has been reported by Plazzotta at al.[12] for aqueous mixtures of C18E20 and C18E100. The PEG chain lengths of the surfactants in ref (12) are exactly twice those of this study (same length ratio). The segregation observed by SAXS, NMR, and differential calorimetry experiments occurred below 7–8 °C, while mixed micelles were preferred at higher temperatures, and the separation was reported to stem from the difference in core-freezing temperatures of the two block copolymer components by Plazzotta et al.[12] This interpretation is supported by the findings of Zinn et al.[43] who have concluded that n-alkyl-PEG micelle cores exhibit a melting point depression that is predominantly controlled by the core size; as the PEG blocks differ in size, the core sizes of the micelles are different.

The foremost reason for the lack of a clear species separation in the simulations is the inability of the MARTINI model to accurately replicate transition temperatures: the micelle core does not freeze in our simulations. Additionally, the simulations here suffer from small size and limited simulation time, that is, surfactant exchange between the micelles does not occur to the full extent. The initial aggregation occurs by neighboring surfactants forming premicelles and then micelles. As the initial configuration is a random mixture of the surfactants, in principle, the resulting distribution in the micelles is even. Thus, any bias toward species separation visible here is underestimated; the observed segregation propensity means that the surfactants have preferred micelle sizes and distribute unevenly to micelles of different sizes even in the absence of core crystallization.

Conclusions

Here, we examined the micellization of medium chain length alkyl-PEG surfactants of different PEG chain lengths in single-surfactant component and mixed-surfactant systems via molecular modeling using a CG MARTINI force-field. We mapped the PEG chain length dependency of the micelle structure, dimensions, and the trends in density changes resulting from different PEG chain lengths. The findings show how the hydrophobic core and the hydrophilic corona size, as well as, density can be tuned via the polymeric head group size. Furthermore, the simulations show indications of the micelles transitioning from the alkyl chain packing-limited aggregation response to the polymeric head group packing-limited aggregation response with the increasing PEG head group length. To our knowledge, the micelles of alkyl-PEG micelles of comparable or longer PEG chains have not been priorly characterized to such a detail. The dependencies and differences in the micelle structure uncovered here as the function of the PEG chain length, especially in terms of micelle hydrophobicity–hydrophilicity profile and density, could have significance, for example, in small molecule partitioning, as well as, usage of these block-copolymer type molecules in solubilization; see refs (7) and (31). Even though the study focused on alkyl-PEG surfactants, the findings result from the polymeric nature of the PEG head group, namely, the scaling of the effective volume occupied by the head group, and thus are generalized to a large family of block-copolymer type surfactants.