12032
Langmuir 2007, 23, 12032-12041
Structure and Interactions in Micellar Solutions: Molecular Simulations of Pluronic L64 Aqueous Solutions Dmitry Bedrov,* Grant D. Smith, and Jinyong Yoon Department of Materials Science and Engineering, UniVersity of Utah, 122 S. Central Campus Dr., Rm. 304, Salt Lake City, Utah 84112 ReceiVed March 13, 2007. In Final Form: July 30, 2007 Coarse-grained, implicit solvent molecular simulations have been conducted to investigate the structure and interactions of L64 Pluronic micelles in aqueous solutions. Simulations of an L64 solution beginning with monodisperse micelles (aggregation number Nagg ) 40 chains) resulted in a narrow Gaussian distribution of Nagg centered around 40. While not fully equilibrated, this distribution supports the supposition that L64 micelles with Nagg ) 40 are representative of the conditions considered and model employed. Detailed analysis of intramicellar monomers distribution and micelle shapes revealed that L64 micelles have a scalene ellipsoidal shape. Additional simulations of solutions containing 125 micelles constrained to have Nagg ) 40 at polymer volume fractions of 0.024 and 0.110 were performed to study micelle-micelle structure factor, single micelle form factor, and total scattering intensity. The ability of various models utilized in analysis of scattering profiles in micellar solutions to describe the structure of the model L64 solutions was investigated. Investigation of the potential of mean force between two micelles reveals that the interactions between micelles are repulsive but on a length scale smaller than the mean micelle diameter, indicating that the micellar shape fluctuations are important in determining intermicellar interactions.
I. Introduction Poly(ethylene oxide) (PEO) and PEO-based copolymers, including PEO-poly(propylene oxide)-PEO (PEO-PPO-PEO) triblock copolymers (also known as Pluronics), are utilized in a wide variety of applications in aqueous environments, including protein crystallization,1,2 modification of surfaces for biocompatibility,3,4 control of particle aggregation in solutions,5,6 and drug delivery.7 At low polymer concentration, Pluronic solutions exist as isolated polymer coils, or unimers.8,9 As the polymer concentration increases at constant temperature, or the temperature increases at constant concentration, the critical micelle concentration or critical micelle temperature is reached and a regime of equilibrium between unimers and spherical micelles is entered.10-16 Small-angle neutron scattering (SANS) studies reveal that Pluronic micelles consist of 20-100 molecules, are typically 4-20 nm in diameter, and have the expected corecorona architecture.17-21 Further reduction in solvent/polymer * Corresponding author. E-mail:
[email protected]. (1) McPherson, A. Methods Enzymol. 1985, 114, 112. (2) Cudney, B.; Patel, S. Acta Crystallogr. 1994, D50, 414. (3) Harris, J. M. In Poly(Ethylene Glycol) Chemistry: Biotechnical and Biomedical Applications; Harris, J. M., Ed.; Plenum Press, New York, 1992, 1. (4) Andrade, J. D.; Hlady, V.; Jeon, S. I. Polym. Mater.: Sci. Eng. 1993, 69, 60. (5) Sakai, T.; Alexandridis, P. Nanotechnology 2005, 16, S344. (6) Alexandridis, P. Curr. Opin. Colloid Interface Sci. 1997, 2, 478. (7) Adams, M. L.; Lavasanifar, A.; Kwon, G. S. J. Pharm. Sci. 2003, 92, 1343. (8) Mortensen K.; J. S. Pedersen. Macromolecules 1993, 26, 805-813. (9) Mortensen, K. In Amphiphilic Block Copolymers; Alexandridis P., Lindman, B., Eds.; Elsevier, New York, 2000. (10) Alexandridis, P.; Hatton, T. A. Colloid Surf. AsPhysicochem. Eng. Aspects 1995, 96, 1-46. (11) Alexandridis, P.; Holzwarth, J. F.; Hatton, T. A. Macromolecules 1994, 27, 2414-2425. (12) Almgren, M.; Brown, W.; Hvidt, S. Colloid Polym. Sci. 1995, 273, 2-15. (13) Wanka, G.; Hoffmann, H.; Ulbricht, W. Macromolecules 1994, 27, 41454159. (14) Hecht, E.; Hoffmann, H. Colloid Surf. AsPhysicochem. Eng. Aspects 1995, 96, 181-197. (15) Chen, W.-R.; Chen, S.-H.; Mallamace, F. Phys. ReV. E 2002, 66, 021403. (16) Pedersen, J. S.; Gerstenberg, M. C. Colloid Surf. AsPhysicochem. Eng. Aspects 1996, 213, 175-187. (17) Liao, C.; Choi, S.-M.; Mallamace, F.; Chen, S.-H. J. Appl. Crystallogr. 2000, 33, 677-681.
Figure 1. Distribution of the probability to find a chain in the micelle of size Nagg obtained from limited CGIS simulations without biasing potential at φ ) 0.110. Also shown are Gaussian fit of the distribution for Nagg > 30 and the theoretical expectation for the equilibrium Gaussian distribution for solution with average Nagg ) 40.
affinity (increase in temperature) or increase in concentration can lead to formation of micellar gels or crystals. Estimates of the size, structure, water content, and aggregation number for Pluronic micelles vary widely. For example, for one of the most characterized Pluronics, L-64 (EO13/PO30/EO13), estimates for the aggregation number at 40 °C vary from 25 to 70.12 Predictions of the water content in the PPO-rich core have ranged from 0% (dry core) up to 60% water.22 Such variations (18) Mortensen, K. J. Phys. Condens. Matter 1996, 8, A103-A124. (19) Pedersen, J. S.; Gerstenberg, M. C. Macromolecules 1996, 29, 13631365. (20) Mortensen, K. Colloid Surf. AsPhysicochem. Eng. Aspects 2001, 183, 277-292. (21) Zhou, D. L.; Alexandridis, P.; Khan, A. J. Colloid Interface Sci. 1996, 183, 339-350. (22) Goldmints, I.; Yu, G.; Booth, C.; Smith, K. A.; Hatton, T. A. Langmuir 1999, 15, 1651-1656.
10.1021/la700742z CCC: $37.00 © 2007 American Chemical Society Published on Web 10/23/2007
Molecular Simulations of Pluronic L64 Micelles
in interpretation of experimental results can be attributed to the accuracy of micellar models that are used for fitting the experimental data. In order to fit a scattering intensity I(q) obtained from SANS measurements, one has to make assumptions regarding the single micelle form factor [F(q)] and the nature of micelle-micelle interactions that define the structure factor [S(q)] of the solution. Models employed for this purpose range from simple hard sphere approximations both for the S(q) and F(q) to complex hydrated core-corona models for F(q) and attractive interactions between Pluronic micelles in aqueous solutions.8-25 While some combination of F(q) and S(q) models poorly represent experimental data and therefore have been discarded, there are a number of those combinations that adequately fit the data, relying on different assumptions and predicting different micellar structure. Molecular simulations of these systems can provide much needed insight into structure of individual micelles as well as micelle-micelle interactions. However, taking into account the slow kinetics of formation26,27 and relatively large dimensions of Pluronic micelles, it is impractical to use conventional atomistic molecular dynamic (MD) or Monte Carlo (MC) simulations to model these systems. For example, simulation of a single hydrated L64 micelle consisting of 40 triblock chains would involve about 150 000 atoms, while simulation of a solution containing multiple micelles would be unrealistic even with utilization of massively parallel architectures. To date the modeling of Pluronics solutions has been limited either to self-consistent field lattice models28-30 or mean-field density functional theory approaches (MesoDyn).31-33 In these studies, an ideal Gaussian chain representation of Pluronic chains and mean-field approximations were employed, while thermodynamic parameters (interaction between species) were either determined semiempirically31,32 or based on incomplete physical models.28-30 While these mean field approches allow computationally expedient prediction of solution morhology and phase behavior, their ability to incorporate important atomistic scale phenomena (i.e., hydrogen bonding, changes in hydration structure, hydrophobic interactions, local conformations, solvent clustering, etc.) operative in these systems is very limited. Recently, we have developed a coarse-grained, implicit solvent (CGIS) model for Pluronic chains in aqueous solution34 using a multiscale modeling hierachy in which computationally expensive degrees of freedom were systematically eliminated while their contribution/influence on thermodynamics, structure, and polymer conformations were implicitly retained in the coarser grained model. The CGIS model for Pluronics has been parametrized to reproduce structural and conformational correlations of PEO and PPO chains in water obtained from extensive atomistic explicit solvent (AES) MD simulations and currently represents the most detailed and (23) Chen, S.-H.; Chen, W.-R.; Mallamace, F. Science 2003, 300, 619-622. (24) Castelletto, V.; Hamley, I. W. Curr. Opin. Colloid Interface Sci. 2002, 7, 167-172. (25) Pedersen, J. S.; Gerstenberg, M. C. Curr. Opin. Colloid Interface Sci. 2002, 7, 158-166. (26) Goldmints, I.; Holzwarth, J. F.; Smith, K. A.; Hatton, T. A. Langmuir 1997, 13, 6130. (27) Hositza, M. J.; Bohne, C.; Alexandridis, P.; Hatton, T. A.; Holzwarth, J. F. Langmuir 1999, 15, 322. (28) Linse, P. Macromolecules 1994, 27, 6404. (29) Hurter, P. M.; Scheutjens, J. M. H. M.; Hatton, T. A. Macromolecules 1993, 26, 5592. (30) Svensson, M.; Linse, P. Macromolecules 1998, 31, 1427. (31) Lam, Y. M.; Goldbeck-Wood, G. Polymer 2003, 44, 3593. (32) Fraaije, J. G. E. M.; Zvelindovsky, A. V.; Sevink, G. J. A. Mol. Simul. 2004, 30, 225. (33) Guo, S. L.; Hou, T. J.; Xu, X. J. J. Phys. Chem. B 2002, 106, 11397. (34) Bedrov, D.; Ayyagari, C.; Smith, G. D. J. Chem. Theory Comput. 2006, 2, 598-606.
Langmuir, Vol. 23, No. 24, 2007 12033
chemically realistic coarse-grained representation of Pluronic chains in aqueous solution. The treatment of water degrees of freedom implicitly allows us to conduct simulations of solutions containing thousands of Pluronics chains forming hundreds of micelles and therefore to address important issues of intra- and intermicellar structural correlations and their relationship to structural correlations directly measureable in experiments. The ability to precisely control solution composition and apply constrains, such as the micelle aggregation number, during MD simulation as well as the ability to decompose structural correlations obtained from the entire systems (e.g., solution scattering intensity) into different contributions makes this approach very suitable and attractive for understanding important correlations in Pluronic micellar solutions. In this paper, we applied this methodology to aqueous solutions of L64 Pluronics. We have determined the micellemicelle interaction, the intramicellar structure, and the structure of micellar solutions as a function of concentration. The scattering intensities obtained from analysis of our simulations have been fitted using different combinations of models and constrains used for representation of single micelle form factor and micellemicelle structure factor. The sensitivity of model parameters to constrains and assumption applied in the fitting of scattering data has been investigated. II. System Description and Simulation Methodology A. Coarse-Grained, Implicit Solvent (CGIS) Model. Details of the parametrization of the CGIS model can be found in ref 34. Here we briefly summarize the essence of the model. In the CGIS model, water is treated as an implicit viscous solvent, while polymer chains are represented as bead-spring polymers where each bead represents a EO (-CH2OCH2-) or PO (-OCH(CH3)CH2-) monomer as a single spherical force center with the center located at the centerof-mass of the corresponding monomer in the atomistic model. The harmonic bonds (1-2 interactions) and bends (1-3 interactions) that connect ether monomer centers-of-mass have been introduced, and the strength of stretching and bending constants have been adjusted to reproduce the 1-2 and 1-3 intramolecular distributions obtained from AES MD simulations of PEO and PPO aqueous solutions. To reproduce intramolecular 1-4 correlations (beads separated by three bonds) we have introduced a specific nonbonded potential that is different from other, the nonlocal (1-5 and greater) intramolecular and intermolecular, nonbonded interactions. Using the inverted Boltzmann approach,35 numerical potentials for 1-4 and other intra- and intermolecular nonbonded interactions were parametrized such that all important structural correlations, characterized by pair distribution functions, obtained in the AES simulations have been reproduced by simulations using our CGIS model. The resulting potential functions for intra- and intermolecular ether-ether interactions were represented numerically and have a nontrivial shape with oscillatory modulations reflecting the waterinduced interactions that are manifested in the AES model. We point out that, while we found that potential functions in our CGIS are independent of solution concentration in a relatively wide range of compositions, they are only valid at 298 K and atmospheric pressure, conditions at which this model has been parametrized. B. System Description and Simulation Protocols. In this work Pluronic L64 (EO13/PO30/EO13) represented by 56 beads (30 PO and 26 EO) has been simulated using the molecular simulation package Lucretius.36 For simulations of micellar solutions, 5000 Pluronic chains were used. Systems with polymer volume fractions φ ) 0.024 and 0.110 have been investigated at 298 K. At this temperature, experiments predict a critical micelle concentration of around φcmc ) 0.10 with average aggregation number (Nagg) of approximately 40 chains per micelle.37 Taking this into account, as well as the fact (35) Reith, D.; Putz, M.; Muller-Plathe, F. J. Comput. Chem. 2003, 24, 1624. (36) http://lucretius.mse.utah.edu/.
12034 Langmuir, Vol. 23, No. 24, 2007
BedroV et al.
that aggregation numbers reported for L64 micelles in water range between 25 and 70 chains,12 we believe that Nagg ) 40 can be considered a reasonable representation of an average L64 micelle in aqueous solution. Initially, the 5000 chains have been divided into 125 micelles each consisting of 40 chains. Each micelle was created by applying an external force on PO monomers and forcing those units toward the micelle center-of-mass. Then the micelles were randomly distributed in a cubic simulation cell (L ) 1000.0 Å for φ ) 0.024 and L ) 587.13 Å for φ ) 0.110) with periodic boundary conditions. Subsequently, Langevin dynamics simulations using the CGIS model were conducted at NVT conditions over several million time steps, which allows multiple L64 chain relaxations inside the micelle.34 During these simulations, no external forces constraining micelle aggregation number were applied. We have also carried out additional simulations of L64 solutions in which all micelles were constrained to Nagg) 40. To implement constraining of Nagg in a way that causes the least perturbation of micellar structure, we analyzed the radial distribution of the centerof-mass of PPO segments (defined as all 30 PO monomers in a Pluronic chain) relative to the micelle core center-of-mass (defined as PPO segments of all chains belonging to the same micelle) from the unbiased simulations described above. We found that for micelles with Nagg) 40 the center-of-mass of PPO segments are always confined within a sphere of radius 40 Å. Therefore, an additional biasing potential that keeps the center-of-mass of each PPO segment inside a sphere of r0 ) 45 Å from the micelle core center-of-mass was applied. Here, the force on a PPO segment center-of-mass due to the biasing potential is calculated and then is redistributed evenly among all 30 PO monomers comprising the PPO segment. Each monomer of a PPO segment is therefore subjected to a biasing potential of the form UPO ) 0, rPPO com < r0 2 PPO UPO ) 0.5kbias(rPPO com - r0) , rcom > r0
(1)
where UPO is the biasing energy acting on each PO monomer, kbias is a force constant (kbias ) 0.5 kcal/mol/Å2), and rPPO com is the distance of the PPO segment center-of-mass from the center-of-mass of the micelle core. Analysis of density profiles of EO and PO monomers around the micelle core center-of-mass from simulations with and without this biasing potential revealed no difference, indicating that the applied biasing potential, while preserving the initial aggregation number (Nagg ) 40), does not perturb micellar shape and intramicellar correlations. We also found that during simulation at any given time step less than 0.1% of the chains has been affected by biasing potential, and the energy contribution due to biasing potential was negligible, indicating that there is no significant driving force for polymer chains to leave the micelle. In order to improve the sampling of micelle-micelle correlations we have coupled the Langevin dynamics simulations with Monte Carlo (MC) displacements of micelles. Here, after every 50 000 integration steps of Langevin dynamics simulation, 1000 MC moves to translate micelles (randomly selected for each move) are attempted. For each MC move, three random numbers between -∆rmax and ∆rmax are generated for translation of a randomly selected micelle along x, y, and z directions, in accordance with metropolis rule. ∆rmax has been adjusted to allow about 40% acceptance rate at each concentration investigated (∆rmax ) 50 Å for φ ) 0.024, and ∆rmax ) 2 Å for φ ) 0.110). The small micelle displacements allowed for MC moves in the more concentrated solution are due to intermicellar interactions, which preclude larger micelle translation. In this concentration regime, the MC moves can only improve the sampling of structural correlations with the nearest neighbor micelles. In order to obtain a good sampling on longer length scale, we have created several statistically independent systems representing this concentration (φ ) 0.110) by taking uncorrelated snapshots from simulations at low concentrations and slowly reducing the cell volume to desired dimensions. (37) Wu, G.; Chu, B.; Schneider, D. K. J. Phys. Chem. 1995, 99, 5094.
In addition to simulations of solutions with multiple micelles of constrained size, we have conducted simulations of two micelles with Nagg) 40 to obtained a potential of mean force (PMF) between those micelles as a function of their separation at infinite dilution (no other micelles or Pluronic chains were present in the simulation cell). In these simulations, an umbrella sampling approach38 with an additional harmonic biasing potential was applied to keep the separation between the centers-of-mass of the micelle cores in a certain window. Thirty umbrella windows were used to sample corecore separations with a 2 Å step in the range between 70 and 130 Å. A self-consistent multiple histogram approach39 was used to analyze the data and to determine the PMF. Finally, we have conducted Langevin dynamics simulations where each micelle was represented as a single spherical force center. In this coarse-grained, implicit solvent, implicit polymer (CGISIP) model the interactions between two force centers (micelles) were described by the PMF obtained from CGIS simulations of two micelles described above. Each system contained 125 spheres in a simulation cell with the same dimensions as used in the CGIS model.
III. Results for Unconstrained Solutions A. Distribution of Micelle Sizes. The distribution of micelle sizes obtained from simulations at φ ) 0.110 is shown in Figure 1. While initially monodisperse with Nagg ) 40, the final distribution is Gaussian with a mean of 40. Also shown in Figure 1 is the distribution expected on the basis of an ideal thermodynamic theory for the formation of spherical micelles40 for a mean size of 40. The ideal distribution is also Gaussian but is significantly broader than that found from simulations. This difference is a strong indication that we cannot access an equilibrium distribution on time scales accessible to our simulations. In common with the behavior observed for φ ) 0.110, which is above the measured cmc, simulations of the more dilute solution (φ ) 0.024), which is well below the measured cmc, resulted in a Gaussian distribution of micelles sizes centered around Nagg ) 40. In contrast, however, we observed a large and continuously increasing concentration of unimers in the low concentration solution, indicating that the majority of the micelles will eventually dissolve given sufficient simulation time. In addition to the two solutions with an initial Nagg ) 40, we performed simulations on a solution at φ ) 0.110 with Nagg ) 80. These micelles rapidly broke up into smaller micelles with a size distribution ranging from 30 to 50. B. Single Micelle Structure. It is typically assumed in analysis of scattering data that Pluronic micelles are spherical. We therefore begin the analysis of the structure of a single micelle by examining the radial number density profiles for EO and PO monomers as a function of separation from the micelle core center-of-mass. Density profiles obtained from simulations at φ ) 0.110 have been averaged for each micelle size separately and are shown in Figure 2 for Nagg ) 40. The density profiles reported in Figure 2 were determined by calculating the number of monomers of a particular type found in a spherical shell of radius r and thickness dr as schematically illustrated in the inset of the Figure 2. Examination of Figure 2 reveals a significant overlap between density profiles of EO and PO monomers, indicating an (apparently) relatively wide boundary between the micelle core and corona. Similar density profiles were observed for all micelles with Nagg > 33. Note that in Figure 2 the density profiles have been averaged assuming that micelles are spherically symmetric. However, the analysis of the principal moments of the micelle core gyration tensor (R2g,i, i ) 1, 2, 3; R2g,1 > R2g,2 > R2g,3), for the (38) Frenkel, D.; Smit, B. Understanding Molecular Simulation: From Algorithms to Applications; Academic Press: San Diego, 2002. (39) Ferrenberg, A. M.; Swendsen, R. H. Phys. ReV. Lett. 1989, 63, 1195. (40) Israelachvili, J. Intermolecular & Surface Forces; Academic Press: San Diego, 1992.
Molecular Simulations of Pluronic L64 Micelles
Langmuir, Vol. 23, No. 24, 2007 12035 Table 1. Components of the Principal Moment of Gyration Tensor for a Single Micelle Obtained for Different Size Micelles at O ) 0.110 from Simulations with and without Biasing Potential
Figure 2. Number radial density profiles of EO and PO units on linear (a) and logarithmic (b) scales as a function of separation from the micelle core center-of-mass as obtained from CGIS simulations of micellar solution at concentrations φ ) 0.024 and 0.110 (from simulations with and without biasing potential). The inset illustrates the definition of the radial density profiles; see the text for details. Arrows show micellar core, corona, and hard sphere radii as obtained from fitting I(q) using the fit 1 model.
Figure 3. A representative snapshot of an L64 micelle with Nagg ) 40 obtained from unbiased CGIS simulations. Dark-colored beads are PO, light colored EO.
φ ) 0.110 solution, where the core is defined by PPO units, reveals that simulations predict that instantaneously the core of the L64 micelles is aspherical, as can be clearly seen in Figure 3. The average principal moments for the L64 micellar core are given in Table 1 for the range of micelle sizes observed in the simulations. As we will discuss below, deviation of L64 micelles from spherical shape has important consequences on micellemicelle interactions.
IV. Results for Solutions with Nagg ) 40 Figure 1 reveals that due to slow kinetics of formation/breakup of micelles it is not possible for us to determine the equilibrium distribution of aggregation numbers on the time scale of our simulations. However, the results discussed above indicate that
Nagg
〈R2g,1〉 (Å2)
〈R2g,2〉 (Å2)
〈R2g,3〉 (Å2)
36 37 38 39 40 40 (Nagg constrained) 41 42 43
378.2 362.5 429.7 457.2 480.6 447.3 539.3 521.7 432.6
200.0 209.6 202.0 203.0 207.4 209.6 209.3 224.1 264.3
134.8 146.3 138.1 138.9 139.1 139.4 139.2 136.6 132.6
Nagg ) 40, while not necessarily the mean micelle size, is a representative size of L64 micelles. Hence, for the sake of clarity in analysis and interpretation, we performed extensive simulations and detailed analysis of single micelle structure and intermicellar interactions in L64 micelle solutions with the micelle size constrained to Nagg ) 40 in the manner described in section II. A. Single Micelle Structure. The ensemble average (over 125 micelles) radial number density profiles for EO and PO monomers as a function of separation from the micelle core center-of-mass obtained from simulations at two compositions are also shown in Figure 2. We find that changes in concentration in the composition range investigated do not result in any noticeable micelle deformation, as illustrated by almost identical density profiles for the two compositions. To address the issue of whether the biasing potential used to maintain Nagg ) 40 perturbs the micelle structure we have compared density profiles obtained from simulations using the biasing potential (eq 1) that constrains micelles to Nagg ) 40 with density profiles obtained for micelles with Nagg ) 40 from simulations without the biasing potential at φ ) 0.110. Figure 2 clearly shows that the density profiles obtained from simulations with and without biasing potentials are nearly identical with the exception of small deviations for low-probability configurations at large distances (Figure 2b). The unbiased simulations show higher probability at larger distances due to Pluronic chains that are in the process of leaving the micelle, which is not allowed in simulations with constrained aggregation number. Table 1 also indicates that R2g,i values for micelles with Nagg ) 40 using the biasing potential are nearly identical to those for unconstrained micelles. Therefore, we conclude that the biasing potential does not perturb the structure and shape of the micelles. B. Potential of Mean Force (PMF) Between Two Micelles. Examination of the density profiles shown in Figure 2 indicates that, as the separation between two micelles (defined by the center-of-mass of their cores) becomes less than 120 Å, an overlap between micellar coronas can be anticipated, and therefore, one can expect to observe interactions between two micelles. Whether those interactions should be attractive or repulsive is not obvious. The excluded volume interactions between interpenetrating coronas of two micelles should result in repulsive interactions. However, we have shown previously34 that effective EO-EO nonbonded interactions are slightly attractive under these thermodynamic conditions; therefore, some attractive interactions between the EO-rich coronas of the two micelles can be anticipated as well. Examination of Figure 4, where the PMF between two micelles is shown as a function of separation between the core centers-of-mass of two micelles, reveals that there is no significant interaction between the two micelles at separations R larger than 100 Å. For R < 100 Å, the PMF monotonically increases indicating repulsive interactions between the micelles. We note that energetic contribution to the PMF from the biasing potential
12036 Langmuir, Vol. 23, No. 24, 2007
BedroV et al.
Figure 4. Potential of mean force between two L64 micelles as a function of their core-core center-of-mass separations R obtained from CGIS simulations.
Figure 6. Planar number density profiles for the pair of micelles along directions parallel (l) and perpendicular (k) to the vector connecting the core centers-of-mass of two micelles. See top panel and discussion in the text for definition.
Figure 5. Representative snapshots of two micelles with core centersof-mass separations: (a) R ) 130 Å and (b) R ) 80 Å.
that constrains the micelle aggregation number is less than 0.1 kcal/mol, even at R ) 70 Å, indicating that changes in shape and orientation between two micelles as a function of their separation are not significantly influenced by the biasing potential. Furthermore, the structure of the micelle pair was found to remain stable and unperturbed for significant simulation time when the biasing potential was turned off. Eventually, the micelles begin to dissolve without the biasing function, which is expected at the very low concentrations of the PMF simulations. Why potentially significant overlap between micellar coronas is not accompanied with changes in intermicellar interactions can be understood if the nonspherical shape of L64 micelles discussed above is taken into account. When micelles are separated by more than 120 Å, they barely interact with each other and their relative orientations are uncorrelated, as illustrated in Figure 5a. At these separations, micelles can have their smallest principal moment of gyration oriented at any angle to the vector connecting the core center-of-masses of two micelles. However, when the separation between two micelles is reduced, the micelles try to
orient such that their smallest principal gyration moments are aligned parallel to the vector connecting two micelles, as illustrated in Figure 5b. This allows minimization of the overlap between the micelles, as shown in Figure 6. In this figure, the planar number density profiles for each micelle along the direction of the vector connecting the two core centers-of-mass (the l direction defined in Figure 5a) are shown for a micelle-micelle separation of 70 Å. The density profiles shown in Figure 6 were calculated by counting the number of monomers of each type (EO, PO) in a planar slab of thickness dl and at separation l from the core center-of-mass of micelle-1 located at the origin of the coordinate system as schematically illustrated in the figure. Even for separations as small as R ) 70 Å, the overlap between coronas is moderate and there is no overlap between micellar cores. Figure 6 also shows the number density profile for micelle-1 in a direction perpendicular to the vector connecting the two micelles (the k direction in Figure 5a) for comparison. EO and PO planar density profiles in this direction are similar to or somewhat more extended than those shown in Figure 3 for the spherically symmetric density profiles. The asymmetry in the density profiles in the l and k directions observed in Figure 6 indicates that, when two micelles approach each other, the orientation of their principal moments of gyration tensor is correlated and that the fluctuations in micelles shape are no longer isotropic. This is further illustrated in Figure 7, where the probability distribution for the 〈R2g,l〉1/2 and 〈R2g,k〉1/2 (components of Rg2 in the l and k directions) are shown for one of the micelles in the presence of another micelle at distance R. For R ) 130 Å, both components (parallel and perpendicular to the vector connecting the micelles) explore the same range of values for the radius of gyration, indicating that shape fluctuations in each direction are equivalent. However, the same figure shows that for R ) 70 Å the radius of gyration in the direction along the vector connecting two micelles is much smaller than the radius of gyration value in the perpendicular direction. C. Structure of L64 Solutions. Structural correlation in micellar solutions obtained from extensive simulations using our combined MD-MC approach on systems containing 125
Molecular Simulations of Pluronic L64 Micelles
Figure 7. Distribution of the micelle core radius of gyration in the l and k directions (see discussion of Figure 5 for definition) in the presence of another micelle at distance R.
Figure 8. Core-core center-of-mass radial distribution function.
micelles with Nagg ) 40 have been determined. Figure 8 shows the micelle core-core center-of-mass pair distribution functions gcore-core(R) for φ ) 0.024. For the larger concentration (φ ) 0.110), a pronounced first peak at separation around 100 Å can be clearly observed, consistent with the fact that the PMF between two micelles (Figure 4) becomes repulsive for separations less than 100 Å. At larger separations the gcore-core(R) shows a shallow minimum and then approaches to 1.0 - 1.0/Nmic ) 0.992, where Nmic is the number of micelles in the system (Nmic ) 125), therefore indicating homogeneous structure beyond the first coordination shell. For the more dilute concentration (φ ) 0.024), gcore-core(R) is only slightly larger than one at small R and then featurelessly approaches 1.0 - 1.0/Nmic at larger R. For both concentrations, gcore-core(R) values show behavior similar to what can be expected for systems of nonaggregating particles (e.g., hard spheres). In addition to CGIS simulations, where each polymer chain was represented explicitly, we have conducted CGISIP simulations, where each micelle was represented as a single force center and their pairwise interactions were described using the micellemicelle PMF. The gCGISIP(R) values (equivalent to the gcore-core(R) obtained from CGIS simulations) obtained from these simulations at equivalent thermodynamic conditions and system configuration (number of micelles and concentration) are also shown in Figure 8. The close agreement between pair distribution functions obtained from CGIS and CGISIP simulations indicates that, in this range of concentrations, micellar solutions can be adequately represented as a collection of single force centers with repulsive
Langmuir, Vol. 23, No. 24, 2007 12037
Figure 9. Micelle-micelle structure factor S(q) as a function of wave vector q obtained from analysis of CGIS simulations using eq 2. Also shown are fits using the HS model (eqs 3 and 4).
interactions given by the PMF shown in Figure 4. Note that the same spherically symmetric PMF was used for simulations of both concentrations, indicating that in this composition range the CGISIP micelle-micelle potential is concentration independent. D. Scattering Properties of L64 Solutions. In this section, we analyze the ability of simple models that are often used for interpretation of experimental scattering data to describe the interand intramicellar correlations obtained from simulations of our micellar systems. The total scattering intensity [I(q)] of the solution is usually the only data accessible from experiment. In order to extract information regarding inter- and intramicellar correlations from experimental I(q), the scattering intensity is fitted assuming that I(q) ) F(q) S(q), where F(q) is the single micelle form factor and S(q) is the structure factor defined by intermicellar correlations. For each property [F(q) and S(q)] a model, usually with several fitting parameters, is assumed. The advantage of analyzing simulations compared to experiments is that F(q) and S(q) can be directly calculated from simulation trajectories without any model assumptions. In order for the observations based upon the scattering properties of our model L64 system to be relevant to real L64 micellar solutions, the Nagg ) 40 micelles comprising our model solution must be reasonably representative of micelles formed in L64 solutions, which we have shown above to be the case. The lack of polydispersity in our model solutions, far from being a disadvantage, is an advantage in determining the applicability of the simple models that are used for interpreting scattering data, since these models usually assume monodisperse solutions. If a model fails to accurately describe scattering in a well-defined monodisperse solution (i.e., our model solutions), it is clearly not applicable to more complex, polydisperse solutions. Below we describe our calculations of the micelle-micelle structure factor [S(q)], single micelle form factor [F(q)], and total solution scattering intensity [I(q)], as well as fitting of these properties with various models reported in the literature. Structure Factor, S(q). S(q) for the solutions investigated has been calculated as
S(q) ) 1 +
∫0∞ (gcore-core(R) - 1) sinqRqR 4πR2 dR
(2)
where gcore-core(R) is the core-core center-of-mass pair distribution function (Figure 8). S(q) is shown in Figure 9. In fitting experimental data, S(q) is often either assumed to be 1.0
12038 Langmuir, Vol. 23, No. 24, 2007
BedroV et al.
(noninteracting micelles) or described using a hard-sphere (HS) approximation. The former usually fails to allow adequate fitting of I(q) data; therefore, the latter is usually employed. S(q) for the HS model can be represented as8
SHS(q) ) 1/(1 + 24φmicP(2qRHS)/(2qRHS))
(3)
Here RHS is the hard-sphere radius, φmic is the micellar volume fraction, and P is given as
P(x) ) R(sin x - x cos x)/x2 + β(2x sin x + (2 - x2) cos x - 2)/x3 + γ(-x4 cos x + 4((3x2 - 6) cos x + (x3 - 6x) sin x + 6))/x5 where
R ) (1 + 2φmic)2/(1 - φmic)4 β ) -6φmic(1 + φmic/2)2/(1 - φmic)4 γ ) (φmic/2)(1 + 2φmic)2/(1 - φmic)4
(4)
In our simulations, where the micelle aggregation number Nagg, the number of micelles, and the volume of the system V are fixed and well defined, no additional assumptions to determine the micellar volume fraction φmic are necessary. The φmic is defined as
φmic ) Nmic4πRHS3/3V
(5)
and therefore only one adjustable parameter, RHS, is available to optimize the fitting of S(q) data using the HS model. Figure 9 shows the HS model fits of S(q) obtained from our simulations. RHS ) 49.6 and 47.1 Å obtained for φ ) 0.024 and φ ) 0.110, respectively, provide a very good representation of simulation data, indicating that the HS model approximation can quite adequately describe intermicellar correlations of L64 micellar solutions in this composition range. Form Factor, F(q). The variety of models used to describe the single micelle form factor is very broad, ranging from uniform sphere to quite sophisticated core-corona models, where chain conformations and polymer density profiles are taken into account.41 Representation of a micelle as a sphere with uniform scattering density, while conceptually simple and requiring only one fitting parameter, is a poor model for Pluronic micelles with highly hydrated coronas and poorly hydrated cores. Moreover, the uniform sphere model cannot provide information about water distribution inside the micelle. Hence, the most commonly applied models for representing the Pluronic micelle form factor are the so-called core-corona models. In the simplest implementation of such a model, the Pluronic micelle is considered to consist of a poorly hydrated, spherical, PPO-rich core with radius Rcore and a well-hydrated, uniform, spherical PEO-rich corona with radius Rcorona, with no overlap between the core and corona regions.42,43 In this case, the form factor can be expressed as42,43
F(q) )
[
(
)
4πR3core 3(sin(x1) - x1 cos(x1)) (F1 - F2) + 3 x3
(
1
)]
3(sin(x2) - x2 cos(x2)) 4πR3corona (F2 - Fw) 3 x3 2
2
(6)
(41) Pedersen, J. S.; Svaneborg, C.; Almdal, K.; Hamley, I. W.; Young, R. N. Macromolecules 2003, 36, 416. (42) Yang, L.; Alexandridis, P.; Steytler, D. C. Langmuir 2000, 16, 8555.
with
x1 ) qRcore, x2 ) qRcorona F1 ) R1FPO + (1 - R1)Fs, F2 ) R2FPO + (1 - R2)Fs R1 ) (3NaggVPPO)/(4πR3core), R2 ) (3NaggVPEO)/(4π(R3corona - R3core))
(7)
FPO ) bPO/VPO, FEO ) bEO/VEO, Fw ) bwpVw where VEO, VPO, and Vw are the volumes of the EO monomer, PO monomer, and D2O molecule, respectively, and bEO, bPO, and bw are their coherent scattering lengths. VPEO and VPPO are the total volumes of EO and PO segments in the chain given by VEONEO and VPONPO, where NEO and NPO are the numbers of EO and PO monomers per chain (NEO ) 26, NPO ) 30). The volumes VEO ) 62.3 Å3, VPO ) 105.7 Å3, and Vw ) 29.8 Å3 were taken from our previous study,34 while bEO ) 4.139 × 10-5 Å, bPO ) 3.307 × 10-5 Å, and bw ) 19.15 × 10-5 Å were obtained by summation of atomic scattering lengths of all atoms comprising the unit, e.g., bEO ) 2bC + bO + 4bH. There are three unknown parameters, Rcore, Rcorona, and Nagg, that need to be fit in this model. Note that these parameters not only define the polymer distribution in the micelle but also define the water content in the core and corona through R1 and R2, which are the volume fractions of PPO and PEO in the core and the corona, respectively. To see how well the core-corona model presented above can describe results from our simulations of L64 micellar solutions, we first calculate the micelle form factor from the simulation trajectories. Note that, since water is treated implicitly in our simulations, there are no scattering centers associated with water molecules. Therefore, in order to predict the form factor for a realistic L64 Pluronic micelle in D2O, we have to take into account correctly the contrast in scattering between D2O and ether monomers. For each micelle the form factor has been calculated as
F(q) )
∑i ∑j bi/bj/
sin(qrij) qrij
(8)
where i and j go from 1 to the total number of monomers in the micelle Nbeads ) Nagg(NEO + NPO), rij is the distance between a pair of beads, and bi/ (i ) EO, PO) represents the excess (relative to the equivalent volume of water) coherent scattering length of EO and PO monomers, bi/ ) bi - bwVi/Vw. The form factor calculated using eq 8 is shown in Figure 10 for the φ ) 0.024 solution. F(q) obtained from the higher concentration (not shown) was nearly identical to the one obtained at φ ) 0.024, confirming that the intramicellar structure for the micelle with Nagg ) 40 does not depend on solution concentration in this composition range. The fit of F(q) by the core-corona model (eqs 6 and 7) is also shown in Figure 10. The fit describes the data nicely for small q values (q < 0.1). However, for higher q values, oscillatory behavior (clearly seen in the inset panel) predicted by the model is not observed either in our data or experiments.37 Parameters for the core-corona model obtained from the F(q) fit are given in Table 2. We remind the reader that while in our simulations all micelles were constrained to have Nagg ) 40, in the fitting of F(q) we allowed this parameter to be fitted without any constrains. Note that the best fit for the core-corona model resulted in Nagg ) 40.02. The radii for the core, Rcore ) 32.3 Å, and corona, Rcorona ) 59.5 Å (for φ ) 0.024), obtained from F(q) (43) Goldmints, I.; Gottberg, F. K.; Smith, K. A.; Hatton, T. A. Langmuir 1997, 13, 3659.
Molecular Simulations of Pluronic L64 Micelles
Langmuir, Vol. 23, No. 24, 2007 12039
Figure 10. Single micelle form factor F(q) as a function of wave vector q obtained from analysis of CGIS simulations at φ ) 0.024 using eq 8. The inset shows the same plot with a log scale for F(q). Also shown is the fit using the core-corona model (eqs 6 and 7). Table 2. Model Parameters Obtained from Fitting fitted property
Nagg
Rcore (Å)
Rcorona (Å)
S(q) (φ ) 0.024) F(q) (φ ) 0.024) I(q), fit 1 (φ ) 0.024) I(q), fit 2 (φ ) 0.024) I(q), fit 3 (φ ) 0.024) I(q), fit 1 (φ ) 0.110)
40.02 40.68 38.09 32.48 41.70
32.23 32.28 34.57 33.94 32.21
59.50 61.0 47.16 40.16 64.70
RHS (Å)
R1
R2
0.904 0.916 0.698 0.629 0.945
0.087 0.081 0.232 0.489 0.068
49.60 50.84 47.16 58.90 48.35
fit appear to be consistent with EO and PO density profiles shown in Figure 2. Note that the RHS obtained from the fitting of S(q) using the HS model lies in between the core and corona radii. Scattering Intensity, I(q). The total scattering intensity, I(q), a property that is directly measured in SANS experiments, can be expressed in terms of F(q) and S(q). Several formulations for defining the I(q) as a function of F(q) and S(q), which take into account polydispersity and shape anisotropy of the micelles, can be found in the literature.44-46 However, in most cases of experimental data analysis described in the literature, a simplified expression I(q) ) F(q) S(q) is assumed. Therefore, in this work we purposely decided to define our I(q) as F(q) S(q) (shown in Figure 11) in order to eliminate additional discrepancy between results of fitting models used in experimental analysis (described below) and the target I(q). Analysis of experimental I(q) data typically involves simultaneous fitting of parameters for the S(q) and F(q) models, often involving additional constraints on the parameters. In our fitting of I(q), we have tried several constraints previously used in the literature. In all fits we used the HS model for S(q) (eqs 3 and 4) and the core-corona model for F(q) (eqs 6 and 7) and adjusted four parameters, RHS, Rcore, Rcorona, and Nagg, with the following additional constraints/assumptions:
fit 1: φmic ) Nchains4πR3HS/3VNagg fit 2: φmic ) Nchains4πR3HS/3VNagg, RHS ) Rcorona fit 3: φmic ) Nchains4πR3corona/3VNagg The fit 2 approach, in which it is assumed that the corona radius is equal to the micelle HS radius, is often used to reduce the number of fitting parameters and has been applied to Pluronic (44) Pedersen, J. S. J. Chem. Phys. 2001, 114, 2839-2846. (45) Pedersen, J. S.; Hamley, I. W.; Ryu, C. Y.; Lodge, T. P. Macromolecules 2000, 33, 542-550.
Figure 11. Total scattering intensity I(q) ) F(q) S(q) as a function of wave vector q obtained from analysis of CGIS simulations: (a) φ ) 0.024 and (b) φ ) 0.110. Also shown are fits using various models.
solutions at similar concentrations to those investigated by us.43 The fit 3 approach has been used in the fitting of SANS data for L64 Pluronic solutions at φ ) 0.025 as a function of temperature.42 In this approach, while the RHS is fitted independently from Rcore and Rcorona, the micellar volume is defined by Rcorona. The best fit of experimental data on L64 solutions using the fit 3 approach resulted in RHS significantly larger than Rcorona (RHS ) 90 Å, Rcorona ) 56 Å, Rcore ) 34 Å),42 which seems to us unphysical. The fits to I(q) using all three approaches for φ ) 0.024 are shown in Figure 11a, while model parameters obtained from these fits are given in Table 2. The best fit is obtained using the fit 1 approach; however, the other two are reasonable as well. The fit 1 approach also results in values of RHS, Rcore, Rcorona, and Nagg that are the closest to the values obtained from fitting the F(q) and S(q) separately, indicating that in this approach a simultaneous fit of S(q) and F(q) does not affect the model parameters. The Nagg for fit 1 is also the closest to 40, which was a constraint in our simulations. Fit 2 slightly underpredicts Nagg and forces Rcorona to be significantly reduced compared to values obtained from fit 1 and F(q) fits. The latter indicates that in fitting I(q) with this approach the S(q) component dominates the fit and forces the RHS and Rcorona (which are constrained to be equal) to values expected for the RHS. In fit 3, Nagg differs 20% from the correct value, Rcorona is even smaller than in fit 2, and RHS is much larger than the Rcorona, qualitatively similar to what has been obtained from fitting SANS data for L64 solutions using the same approach.42 The parameters obtained from fit 3 (46) Frenkel. D.; Vos, R. J.; de Kruif, C. G.; Vrij, A. J. Chem. Phys. 1986, 84, 4625-4630.
12040 Langmuir, Vol. 23, No. 24, 2007
BedroV et al.
Figure 13. Local volume fractions of EO, PO, EO + PO monomers and water as a function of separation from the micelle core centerof-mass calculated at φ ) 0.024 using density profiles from Figure 2. Also shown are the volume fractions assumed/obtained using fit 1-3 models. Table 3. Water Volume Fraction in the Core and Corona for O ) 0.024 core
corona
core φ ) ∫R0coreφw(r) dr φ )1 - R1 φ ) ∫RRcorona φw(r) dr φ )1 - R2
fit 1 fit 2 fit 3
Figure 12. (a) F(q) and (b) S(q) obtained using parameters obtained from fits 1-3 of I(q). Also shown are the data obtained from analysis of simulation trajectories and their direct fits with core-corona [for F(q)] and HS [for S(q)] models.
differ significantly from parameters obtained from individual fits of S(q) and F(q), indicating that constraining the micellar volume to be defined by Rcorona is a poor assumption that can result in unreasonable parameters with unphysical trends. Interestingly, Rcore is only weakly dependent on the model used. Another assessment of the adequacy of the I(q) various fitting protocols can be made by examining predictions of the corecorona [for F(q)] and HS [for S(q)] models using parameters obtained from I(q) fitting. In Figure 12 we compare F(q) and S(q) for φ ) 0.024 predicted by the models with parameters obtained from fits 1-3, directly calculated from our simulations and from a separate fitting of F(q) and S(q) with the corecorona and HS models, respectively. While fits 1-3 shown in Figure 11 give an adequate description of the I(q), the constraints employed in some of the fits significantly affect the quality of F(q) and S(q) descriptions, as can be seen from Figure 12. Parameters obtained from fit 1 provide the best prediction of F(q) and S(q) contributions. Fit 2 and particularly fit 3 show significant deviations from the actual F(q) and S(q), indicating that in these approaches F(q) and S(q) contributions to the I(q) are partitioned incorrectly. Taking into account the comparison of different fit models discussed above, we conclude that fit 1 is the most reliable approach, and have applied this approach for fitting the I(q) data at φ ) 0.110. Figure 11b shows that the fit obtained using the fit 1 model provides an adequate description of I(q) at this concentration, although not as good as at the lower concentration. Since F(q) is nearly identical for the two compositions and can be accurately described by the core-corona model, this degrada-
0.35 0.41 0.40
0.08 0.30 0.37
0.91 0.80 0.73
0.92 0.77 0.51
tion in the quality of the fit can be attributed to a poorer description of intermicellar correlations [S(q)] using the HS model at the higher concentration. The latter can be also seen in Figure 9, where the intermicellar S(q) is fitted by the HS model for both concentrations. Because of the nonspherical micellar shape and the ability of micelles to deform upon interaction with other micelles, it can be anticipated that, as the solution concentration increases, the HS model might become inadequate for description of intermicellar correlations. We will address this issue in more detail in a future publication. Water Distribution. One of the main advantages of the corecorona model described above for representation of the single micelle form factor is the ability to predict the water content inside the core and corona regions. Therefore, it is instructive to compare predictions of the different fitting protocols (fits 1-3) for the water content inside the core and corona regions with those directly obtained from simulations. In Figure 13 we show the mean volume fraction of EO, PO, and EO + PO monomers as a function of separation from the micelle core center-of-mass obtained from our CGIS simulations of the φ ) 0.024 solution. These volume fractions were obtained using the radial density profiles reported in Figure 2. Also shown in this figure is the volume fraction of water estimated as φw ) (Vtotal - nEOVEO - nPOVPO)/Vtotal under the assumption that water occupies all the remaining volume not occupied by EO or PO units in our CGIS simulations. Using values for Rcore and Rcorona obtained from fits 1-3 as core and corona boundaries, a volume fraction of water in the core and corona, φcore and φcorona, can be calculated by integration of φw(r), as shown in Figure 13. These volume fractions are given in Table 3. Fitting models fits 1-3 assume a homogeneous distribution of water inside each region (core and corona), and the corresponding water content can be obtained as 1 - R1 and 1 - R2 (R1 and R2 are reported in Table 1) for the core and corona, respectively. Water volume fractions predicted by the fitting models are also shown in Figure 13 and given in Table 3.
Molecular Simulations of Pluronic L64 Micelles
Langmuir, Vol. 23, No. 24, 2007 12041
approach. For example, fitting of PO volume fraction profiles with hyperbolic tangent φ(r) ) tanh(2r/w) allowed interfacial width (w) of 12.4 Å for the cone analysis, which is noticeably narrower than the value of 18.6 Å obtained for spherical averaging. However, even in the analysis where the micelle shape anisotropy is taken into account, the interfacial region of EO and PO monomers overlap is still relatively broad, and a sharp rise in water content from the value at the core center coincides with the onset of increasing volume fraction of EO monomers. Finally, similar analysis of the water volume fraction in the corona φcorona region indicates that predictions from fits 2 and 3 significantly underestimate the corona hydration. The reduced values of Rcorona [compared to the F(q) fit] obtained for fits 2 and 3 do not reflect the real extension of PEO segments, which leads to significantly increased volume fraction of EO monomers in the corona.
V. Conclusions
Figure 14. Local volume fractions of EO, PO and water in the core-corona interfacial region obtained from spherical and cone averaging (see the text and inset for definition) at φ ) 0.024. Lines show hyperbolic tangent fit for PO profiles.
Upon comparing the volume fraction profiles, it is clear that fit 1 gives the best prediction of the water content (φw ≈ 0.08) for the region well inside the PPO-rich core (r < 15 Å), while fits 2 and 3 significantly overestimate the water content in this region. However, φw(r) obtained from our simulations begins to increase sharply for r > 15 Å and reaches a value of around 0.60 at r ) Rcore. This results in a volume fraction inside the core φcore (calculated using a corresponding Rcore for each model and given in Table 2) of around 0.40, indicating that fits 2 and 3 give a better estimate of the volume fraction of water inside of a sphere with Rcore radius. Such inconsistence can be explained by realizing that (a) L64 micelle shape is not spherical and (b) the interface between core and corona is not sharp, as it is assumed in all fitting models. Figures 2 and 13 indicate that there is significant overlap between core (PO) and corona (EO) regions. However, spherically averaged density and volume fraction profiles shown on these figures do not take into account that the instantaneous shape of each micelle is ellipsoidal. Spherical averaging broadens the core-corona interface (or overlap region), even if the interface would be a perfect step function. We therefore enhanced our analysis of density profiles by considering narrow cones (θ ) 20°) oriented along the eigenvectors of the instantaneous principal moments of the core gyration tensor, as schematically illustrated in the inset of Figure 14. This allows us to avoid having both core and corona regions inside the same shell due to mismatch in the micellar shape and averaging. To allow averaging of density profiles over all cone orientations (three eigenvectors for each micelle), the coordinate origin (r ) 0) was set at the point where φPO ) φEO for each core, which we define as a core-corona interface. In Figure 14, we compare volume fraction profiles of EO, PO, and water near the interface obtained from spherical averaging and when the shape anisotropy is taken into account. The figure clearly shows that locally the core-corona interface is narrower than is predicted using the spherically symmetric
Extensive simulations using a coarse-grained, implicit solvent model of L64 micellar solutions with aggregation number constrained to 40 have been conducted in the composition range below 0.110 to investigate intra- and intermicellar structural correlations, intermicellar interactions, and the ability of simple models to describe the scattering data of these solutions. We find that L64 micelles with Nagg ) 40 instantaneously are not spherical. Anisotropy in shape allows micelles to minimize their interaction at shorter separations by aligning their smallest gyration moments parallel to the vector connecting two interacting micelles. The potential of the mean force for two micelles in dilute solution shows an onset of repulsive interactions for separations less than 100 Å. There are no signs of attractive interactions in the potential of mean force between L64 micelles at any separation. The micelle-micelle structure factor and the single micelle form factors obtained from simulations of L64 solutions were fitted separately using the hard-sphere and core-corona models, respectively. We find that the hard-sphere model provides an adequate description of the solution structure factor, while the core-corona fit of the micelle form factor results in an accurate description of the single micelle structure. However, when parameters of two models are fitted simultaneously to reproduce the total scattering intensity of solution, a strong sensitivity to constraints/assumptions used in the fit and incorrect partitioning between form factor and micelle-micelle structure factor are observed. We find that the most reasonable/accurate model parameters are obtained when the hard-sphere, core, and corona radii and the micelle aggregation number are fitted independently (with out constraints) with the assumption that the micelle volume is defined by the micelle hard-sphere radius. Application of additional constraints equalizing the hard-sphere and corona radii or the assumption that micellar volume is defined by the corona radius lead to inaccurate prediction of aggregation number, qualitatively incorrect micellar dimensions, and overestimation/ underestimation of water content in the micellar core/corona. Finally, in the composition range investigated, the structure of micellar solution (characterized by intermicellar distribution) can be adequately predicted using the coarse-grained, implicit solvent and implicit polymer simulations, where micelles are represented as a single force center with pairwise interaction described by the potential of mean force between two micelles. Acknowledgment. Authors would like to acknowledge support of the National Science Foundation through grants CRCCHE0304807 and MRSEC-DMR0213918. LA700742Z
