598
J. Chem. Theory Comput. 2006, 2, 598-606
Multiscale Modeling of Poly(ethylene oxide)-Poly(propylene oxide)-Poly(ethylene oxide) Triblock Copolymer Micelles in Aqueous Solution Dmitry Bedrov,* Chakravarthy Ayyagari, and Grant D. Smith Departments of Materials Science & Engineering and Chemical Engineering, UniVersity of Utah, 122 South Central Campus DriVe, Room 304, Salt Lake City, Utah 84112 Received December 28, 2005
Abstract: We present a multiscale modeling approach for simulation of poly(ethylene oxide)poly(propylene oxide)-poly(ethylene oxide) triblock copolymer micelles in aqueous solution. We rely on systematic elimination of computationally expensive degrees of freedom yet retain implicitly their influence on the remaining degrees freedom in a coarser-grained model. Quantum chemistry (QC) calculations, atomistic explicit solvent (AES) molecular dynamics (MD) simulations, and coarse-grained implicit solvent (CGIS) simulations have been applied to investigate physical properties of these important self-assembling triblock copolymers. High-level QC calculations have been used to parametrize classical atomistic force fields that implicitly take into account and reproduce the important energetic and structural features due to correlations of electronic degrees of freedom. AES MD simulations utilizing the QC-based force fields have been used to provide structural and conformational properties of polymers in aqueous solution which were subsequently used for parametrization of the CGIS model using the Inverted Boltzmann method. The CGIS simulations were then employed to investigate structural properties of two PEO-PPO-PEO micelles (EO13-PO30-EO13 and EO99-PO65-EO99 also known as Pluronic L64 and F127, respectively) in aqueous solution.
I. Introduction Poly(ethylene oxide) (PEO) and PEO-based copolymers, including PEO-poly(propylene oxide)-PEO or PEOPPO-PPO, 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 room and physiological temperatures PEO of any molecular weight is soluble in water at any concentration, while PPO has very limited solubility in water for short oligomers only. In dilute solutions at room-temperature PEO-PPO-PEO copolymers can usually dissolve in water. As the temperature or concentration of polymer increases PEO-PPO-PEO triblocks form isolated micelles (usually spherical in shape but * Corresponding author e-mail:
[email protected].
10.1021/ct050334k CCC: $33.50
with other shapes possible depending on the triblock architecture) with PPO units in the core and PEO chains forming the corona. Upon further increase in polymer concentration or temperature, the supramicellar structures such as disordered gels and crystals as well as cylindrical and lamellar morphologies can form. The formation of micelles and supramicellar structures is quite sensitive to solvent quality. Slight changes in temperature,8 small addition of cosolvent,9 or change in pH10 can significantly alter the phase behavior of PEO-PPO-PEO aqueous solutions by shifting the delicate balance between hydrophilic and hydrophobic interactions in the system. A typical Pluronic chain contains about 50-300 ether monomers, while a Pluronic micelle formed in aqueous solution consists of up to 100 chains and has a diameter of several nanometers.11 Dynamic light scattering measurements of micellization kinetics for Pluronic chains in water upon © 2006 American Chemical Society
Published on Web 03/28/2006
Multiscale Modeling Approach
increasing solution temperature have shown that formation of micelles is a complex multistage process.12,13 For example, for L64 Pluronic (EO13-PO30-EO13) the micellization was found to consist of three processes: (1) The fast process associated with unimers joining/leaving incipient micelles with characteristic relaxation time on the order of microseconds. (2) The relaxation of micelle size distribution, a process during which the initially formed micelles adjust their size in accordance with equilibrium micelle size distribution corresponding to imposed thermodynamic conditions. During this process, the initially formed micelles can collide and merge into one larger micelle, or a large micelle can split into smaller micelles, or a micelle can dissolve into unimers which can then join other micelles. While the nature of operative mechanism for this process is a topic of ongoing debate, the characteristic relaxation time for this process for L64 Pluronic is on the order of milliseconds. (3) A process observed only for some Pluronics at temperatures much higher than the critical micellization temperature and has been associated with clustering of micelles. This process has operative relaxation times on the order of seconds. It was also found that all these kinetic processes slow dramatically with increasing molecular weight of the Pluronic chain.12 Taking into account the slow kinetics of formation and relatively large dimensions of Pluronic micelles it is impractical to expect conventional, brute force atomistic molecular simulations (molecular dynamics (MD) or Monte Carlo (MC)) to be able to model these systems. For example, a simulation of a single L64 micelle consisting of 40 triblock chains and hydrated by water would involve about 150 000 atoms. For such a large system accessing the operative time scales for the first (fast) micellization process is unrealistic even with utilization of massively parallel architectures. To date the modeling of Pluronics solutions has been limitted either to self-consistent field lattice models14-16 or mean field density functional theory approaches (MesoDyn).17-19 In both cases, an ideal Gaussian chain representation of Pluronic chains and mean-field approximations were employed, while thermodynamic parameters (interaction between species) were either determined semiempirically17,18 or based on incomplete physical models.14-16 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 limitted. In this paper we discuss an alternative approach to study Pluronic micelles in aqueous solution using a multiscale modeling hierachy in which computationally expensive degrees of freedom are systematically eliminated, while their contribution/influence on thermodynamics, structure, and conformations is implicitly retained in the coarser-grained model. We demonstrate for L64 and F127 Pluronics that our multiscale modeling approach can be used to efficiently study Pluronic micelles both at atomistic and coarse-grained levels.
J. Chem. Theory Comput., Vol. 2, No. 3, 2006 599
II. Development of the Multiscale Model for Pluronic Chains in Aqueous Solution The strategy of our multiscale modeling relies on systematic elimination of computationally expensive degrees of freedom while retaining implicitly their influence on the remaining degrees freedom in the coarser-grained model. For each level of coarsening we employ the corresponding most accurate and effective method/simulation technique available to investigate physical properties of the system at that level. Then, using the obtained information we parametrize a coarser model that incorporates all essential physics/ phenomena observed at the finer level. Here, we briefly outline the general strategy of our multiscale modeling approach: (1) We perform high-level quantum chemistry (QC) calculations on model compounds and their clusters to determine relative conformational and binding energies. Due to computational expense of high-level QC calculations this investigation is limited to short (one-three repeat units) oligomers and clusters with relatively few molecules (e.g., an ether oligomer with two water molecules). Optimal geometries, relative conformational energies, and cluster binding energies as a function of molecules separation and orientation obtained from QC calculations are subsequently utilized in parametrization of fully atomistic force fields for classical MD simulations where all subatomic scale degrees of freedom and correlations are averaged out and represented implicitly through effective nonbonded van der Waals interactions and through simplified intramolecular potentials that describe bond and bend vibrations and dihedral rotations. (2) Using parametrized QC-based force fields we conduct extensive atomistic, explicit solvent (AES) MD simulations of low molecular weight polymeric systems (e.g., PEO and PPO oligomers in aqueous solutions). These simulations serve two purposes: they (a) allow us to validate the developed force fields by direct comparison of thermophysical, structural, and dynamical properties obtained from MD simulations with data available from experiments and (b) provide us with a vast amount of molecular scale structural, conformational, and dynamic properties/correlations that help to identify important interactions/correlations which can be used to parametrize coarse-grained models. (3) Using conformational and structural properties obtained from AES MD simulations we parametrize the coarsegrained, implicit solvent (CGIS) model in which polymer chain segments or monomers are represented as single force centers (beads) and solvent is treated implicitly. Langevin dynamics or Lattice-Boltzmann simulations are then conducted using the CGIS representation of polymer chains. These simulations are about 3 orders of magnitude computationally less expensive than simulations with the AES model, therefore allowing us (a) to simulate realistic size polymers (e.g. Pluronic chains) and (b) to significantly extend the accessible time scales (e.g., extend the simulation time window up to the range where the first process of micellization can be resolved for some Pluronics). (4) Finally, equilibrium configurations of polymeric systems obtained from CGIS simulations can be reverse mapped back to the AES model, and AES MD simulations
600 J. Chem. Theory Comput., Vol. 2, No. 3, 2006
can be conducted in order to investigate atomistic scale structural and dynamic properties/correlations in the system. Below we discuss in detail each step of our multiscale modeling strategy as applied to Pluronics in aqueous solution. A. Quantum Chemistry Calculations and Force Field Fitting. To parametrize the atomistic force field for PEO and PPO in aqueous solution previously we have conducted quantum chemistry calculations for 1,2-dimethoxyethane (DME) and 1,2-dimethoxypropane (DMP) compounds in the gas phase and studied the conformational energetics and rotational energy barriers in these oligomers. We also have calculated complexation energy of these ethers with water as a function of their separation for two different paths.20 Details of these calculations can be found in refs 20-23. Molecular geometries, conformational and complexation energies, and electrostatic fields surrounding the model molecules determined from QC were then utilized to parametrize intra- and intermolecular interactions for the classical MD simulations in which each atom is represented as a single point force center. Local intramolecular interactions were described by bonds, bends, and dihedrals, while nonbonded intra- and intermolecular interactions were described by Lennard-Jones or Buckingham potentials and Coulomb interactions. The parameters for bonds, and bends, were fitted to reproduce optimal geometries of compounds obtained from QC calculations. Partial atomic charges were fitted to reproduce the dipole moment of and an electrostatic field around the molecule obtained from QC calculations. Nonbonded van der Waals interactions were then fitted to reproduce QC binding energies for the molecular clusters (e.g. paths between DME or DMP and water). Finally, the dihedral potentials for ethers were parametrized to reproduce relative conformational energies and rotational barriers obtained from QC calculations for model compounds such as DME and DMP. The force field parameters and details regarding force field parametrization in general and specifically for ether/water systems can be found in our previous work.20-23 In this study, the parameters for the Lennard-Jones interactions between PPO atoms and water were the same as used for PEO and water20 for identical atom types (e.g. ether oxygen, ethylene carbons, ethylene hydrogen, etc.), while the distribution of atomic partial charges for each ether varied due to difference in their chemical structure. B. Atomistic, Explicit Solvent MD Simulations. The force fields described above have been utilized extensively in MD simulations of PEO and PPO melts24 as well as PEO in aqueous solution.25,26 Specifically for the purpose of parametrization of the CGIS model for Pluronic chains in aqueous solution additional AES MD simulations of PEO and PPO oligomers in water were performed at atmospheric pressure and 298 K: • PEO/water solutions with PEO Mw ) 530 (12 repeat units) at 0.17 and 0.52 weight fraction of PEO • PPO/water solutions with PPO Mw ) 349 (6 repeat units) at 0.17 and 0.52 weight fraction of PPO • PPO/PEO/water solutions with one PPO Mw ) 349 (6 repeat units) chain and several PEO Mw ) 215 (6 repeat units) chains at 0.52 weight fraction of ether molecules.
Bedrov et al.
Figure 1. Intermolecular monomer-monomer pair distribution functions for ethers in aqueous solution at 298 K and 0.52 weight fraction of ethers as obtained from AES MD (symbols) and CGIS (lines) simulations.
All AES MD simulations were performed using the simulation package Lucretius.27 Each system contained several ether chains and up to 1000 water molecules (depending on concentration) represented using the TIP4P water model.28 Bond lengths have been constrained using the SHAKE algorithm.29 A reversible multiple time step integrator30 was employed with a 0.75 fs time step for bends and torsions, a 1.5 fs time step for nonbonded and the real part of electrostatic interactions within 6.0 Å, and a 3 fs time step for nonbonded interactions for separations between 6.0 Å and 10.0 Å as well as for the reciprocal part of the Ewald summation.31 Systems were equilibrated over 3 ns using an NPT ensemble, while production runs of 10 ns were performed in an NVT ensemble with V equal to the average value obtained from the NPT simulations. These AES MD simulations provided us with structural and conformational properties of ethers in aqueous solution. In particular we focused on intermolecular pair correlation functions between monomer (defined as -CH2-O-CH2for EO, and -O-CH(CH3)-CH2- for PO) centers-of-mass in solution. In Figure 1, we show EO-EO, PO-PO, and EO-PO intermolecular pair correlation functions, g(r), as obtained from AES MD simulations for solutions with 0.52 weight fraction of ethers. The difference between EO-EO and PO-PO g(r) at short separations reflects the larger size of the PO monomer. Particularly striking is the EO-PO g(r) which clearly shows unfavorable interaction of EO monomers with PO monomers. Our initial expectation was that the EO-PO g(r) would lie somewhere between EO-EO and PO-PO g(r), which is clearly not the case. This indicates that despite the amphiphilic nature of PEO, PEO prefers to minimize its interaction with the hydrophobic PPO chain in the presence of water. This result, while unexpected, appears to be consistent with density functional theory modeling of Pluronic solutions where it was found that a relatively large positive χ (effective Flory-Huggins parameter) is necessary to describe PEO-PPO interaction in water.17,18 The PEOPPO χ used in these studies was an order of magnitude larger than the one estimated from SANS measurements for melts
Multiscale Modeling Approach
J. Chem. Theory Comput., Vol. 2, No. 3, 2006 601 Table 1. Intramolecular Parameters for Ethers Bonds and Bends in the CGIS Model Ubond(R) ) 0.5*Kbond*(R0-R)2 bond type EO-EO PO-PO EO-PO
Kbond (kcal/mol/Å2)
R0 (Å)
30.0 20.0 25.0
3.20 3.52 3.36
Ubend(R) ) 0.5*Kbend*(θ0-θ)2
Figure 2. Intramolecular distributions of distances between neighboring monomers (1-2), monomers separated by one monomer (1-3), and monomers separated by two monomers (1-4) for PEO and PPO chains in water as obtained from AES (symbols) and CGIS (lines) simulations.
of PEO-PPO block copolymers32 and about two times larger than the PPO-PPO χ in water. In addition to intermolecular correlations we have analyzed intramolecular correlations for PEO and PPO chains in water. In Figure 2, we show the distribution of intramolecular distances between centers-of-mass of two adjacent (1-2 interaction), separated by one (1-3 interaction), and separated by two (1-4 interaction) EO or PO monomers. These distributions characterize local conformations of PEO and PPO chains in water and need to be accurately reproduced in the CGIS model. More detailed discussion of local structure and conformations in ether/water solutions as a function of temperature and concentration can be found in our previous work.25 C. Parametrization of the Coarse-Grained, Implicit Solvent 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 or PO monomer. Our initial guess was to represent EO (-CH2-O-CH2-) and PO (-O-CH(CH3)-CH2-) monomer as a single spherical force center with the center located at the center-of-mass of the corresponding monomer in the atomistic model. While this was the simplest and intuitively obvious geometric mapping between two models, we were prepared to use a more complicated mapping (e.g. representing a monomer with two or three force centers)
bend type
Kbend (kcal/mol)
θ0 (deg)
EO-EO-EO PO-PO-PO PO-EO-EO PO-PO-EO
10.0 5.0 7.5 7.5
128.0 110.5 119.25 119.25
between two models. However, as we illustrate below, the obtained agreement between CGIS and AES models did not require any additional and more complex structural mapping. The harmonic bonds and bends that connect ether monomer centers-of-mass have been introduced, and the strength of stretching and bending constants has been adjusted to reproduce the 1-2 and 1-3 intramolecular distributions obtained from AES MD simulations as illustrated in Figure 2. The force constants, equilibrium bond lengths, and bend angles in the CGIS model for PEO and PPO are given in Table 1. To reproduce 1-4 interactions one can either introduce a dihedral potential or a specific nonbonded potential that can be different from other nonbonded interactions. In this work, we chose the latter approach, and the resulting nonbonded potentials for intramolecular EO-EO and PO-PO 1-4 interactions are shown in Figure 3. We note that, in principle, the local intramolecular interactions (1-2, 1-3, and 1-4) can depend on nonbonded intra- and intermolecular interactions and vice versa, necessitating iterative fitting of intramolecular and nonbonded interactions. However, in this work, we found that the description of the nonlocal (1-5 and greater) nonbonded intramolecular and nonbonded intermolecular interactions had essentially no influence on the description of the local (1-2, 1-3, and 1-4) intramolecular interactions shown in Figure 2.
Figure 3. Converged 1-4 nonbonded interactions in the CGIS model.
602 J. Chem. Theory Comput., Vol. 2, No. 3, 2006
Bedrov et al.
Figure 4. The intermolecular EO-EO nonbonded interaction in the CGIS model obtained after several iterations using the Inverted Boltzmann method.
To parametrize nonlocal (1-5 and greater) nonbonded intramolecular and nonbonded intermolecular interactions we have employed the Inverted Boltzmann approach33 in which a numerical potential is determined iteratively based on deviation of properties obtained from short CGIS simulations using the current iteration of the potential, Ui(r), and the desired target properties obtained from AES MD simulations. The intermolecular g(r) shown in Figure 1 for EO-EO, POPO, and EO-PO monomers obtained from AES MD simulations were chosen as the target functions, gtarget(r), that we wish to reproduce from simulations using the CGIS model. These gtarget(r) reflect all important correlations between ethers in aqueous solution, and therefore the coarsegrained potentials that reproduce these target functions would implicitly incorporate important water-induced interactions between the ethers. The initial guess for all interactions was a Lennard-Jones type potential as illustrated in Figure 4 for the EO-EO interaction. Short CGIS simulations using Langevin Dynamics (LD) were performed for the system exactly equivalent to the one used in AES MD simulations (the same number of ether chains and volume of the simulation cell), and the corresponding gi(r) were determined. These simulations must be sufficiently long to allow accurate determination of gi(r). In our CGIS model all nonbonded interactions were represented by numerical functions which are not restricted to any particular analytical form. Based on deviations between gi(r) and gtarget(r) we predict the new numerical potential Ui+1(r) for the i+1 iteration of the CGIS model using the following expression Ui+1(r) ) Ui(r) + w(r)kBTln
(
gi(r)
gtarget(r)
)
(1)
where w(r) ) exp(-r) is a weighting function and kB is the Boltzmann constant. Then, we repeat the CGIS simulations to obtain the new gi+1(r) pair correlations and continue this procedure until gi(r) reproduces gtarget(r) as accurately as possible. In Figure 4, we illustrate this process by showing Ui(r) for EO-EO at several iterations, while in Figure 1 we show the agreement between the final (8th) iteration gi(r) obtained from CGIS simulations and gtarget(r) from AES MD
Figure 5. Converged nonbonded interactions in the CGIS model obtained using the Inverted Boltzmann method.
simulations. In the fitting nonbonded interactions we first fitted the EO-EO and PO-PO interactions to match gtarget(r) from AES MD simulations of solutions containing only PEO or PPO chains, respectively. Using these potentials for interactions between beads of the same type the cross-term EO-PO nonbonded interactions were fitted to reproduce EO-PO gtarget(r) obtained from the PEO/PPO/water system. In Figure 5 the final ether-ether nonbonded interactions obtained for the CGIS model are shown. We point out that all potentials have a nontrivial shape with oscillatory modulations due to structure of the hydrating water around the ethers that results in water-induced interactions which are manifested in gtarget(r) and therefore reflected in the CGIS potentials. We also note that the EO-EO monomermonomer interaction has a reasonably deep minimum at separations around 4.8 Å indicating an apparent short-range attraction between EO monomers in water. However, this minimum is followed by a positive maximum and then by a very shallow second minimum at larger separations which results in an overall positive second virial coefficient B ) 71 Å3 for EO monomers consistent with the solubility of PEO in water. The PO-PO potential shows much stronger attraction resulting in a negative second virial coefficient B ) -56 Å3 consistent with the poor solubility of PPO in water. Finally, the EO-PO monomer-monomer potential is repulsive for almost all separations yielding B ) 193 Å3. Finally, the monomer friction coefficient for CGIS LD simulations of Pluronic micelles was set to yield a viscosity of the effective solvent 100 times smaller than that of water thereby facilitating equilibration and sampling of the micelle structure. To establish an exact correspondence between time scales in the CGIS simulations and real time units a headto-head comparison of motion and relaxation in the micelles should be conducted for CGIS and AES simulations which are currently underway. In this paper, we report all dynamical properties as a function of the number of integration steps. All systems were simulated using integration time step ∆t ) (/mσ2)1/2 ) 0.003 where ) 0.189 kcal/mol is the lowest energy in the EO-EO nonbonded interaction shown in Figure 5, m is the mass of the EO monomer, and σ ) 4.213 Å is the size of the EO monomer defined as the smallest
Multiscale Modeling Approach
J. Chem. Theory Comput., Vol. 2, No. 3, 2006 603
Figure 6. The average time autocorrelation function for relaxation of angles formed by two intermolecular vectors connecting the micelle center-of-mass and PEO-PPO junction point for each chain as obtained from CGIS LD simulations of L64 and F127 micelles.
separation between two EO-EO monomers at which their nonbonded interaction is zero (see Figure 5).
III. Multiscale Modeling of Pluronic Micelles A. Coarse-Grained, Implicit Solvent Simulations of Pluronic Micelles. Following parametrization of the CGIS model the single micelle simulation for L64 (EO13-PO30EO13) and F127 (EO99-PO65-EO99) Pluronics has been set up. For each Pluronic micelle we used a number of chains within the range of aggregation numbers estimated from SANS measurements.34-36 Specifically, we used Nagg ) 40 and 30 for L64 and F127, respectively. Initially, randomly generated configurations of Pluronic chains were simulated with an external field that imposed an attractive interaction between the center of the simulation cell and all PO monomers. This allowed us to quickly assemble a micelle with all PO monomers within a spherical micelle core and EO segments extending into the solution making the micelle corona. Following this forced “micellization” we turned off the external field and allowed the micelles to equilibrate their structure during the CGIS simulation. Our primary objective was to obtain equilibrium configurations of a micelle with a fixed Nagg, and, therefore, it was necessary to confirm that the complete decorrelation/relaxation of the micelle structure occurs on time scales accessible to the CGIS LD simulations. To characterize relaxation of the micellar structure we have defined vectors that connect the micelle core center-of-mass and PEO-PPO junction points for each chain. We then calculated the time autocorrelation function (ACF) ACF(t) )
〈cosθ(t)cosθ(0)〉 - 〈cosθ〉2 〈cos2θ〉 - 〈cosθ〉2
(2)
where θ is the angle between two intermolecular vectors, and 〈 〉 represents averaging over all possible intermolecular vector pairs in the micelle. When the ACF, shown in Figure 6, decays to zero, the micelle configuration has completely decorrelated relative to the time zero orientation of Pluronic
Figure 7. Representative snapshots of L64 and F127 micelles in the CGIS model.
chains in the micelle. Figure 6 illustrates that even for the relatively long chains of the F127 the micelle relaxation time is on the order of 106 integration time steps, while simulations using the CGIS model can be routinely performed over 108109 integration steps. The CGIS simulations carried out over multiple characteristic micellar relaxation times (defined above) provide us with equilibrium configurations of Pluronic micelles independent of how the initial micelle configuration has been created for a given Nagg. Representative configurations of L64 and F127 Pluronic micelles obtained from CGIS simulations are shown in Figure 7. We note that over the entire simulation for both micelles (108 time steps) none of the Pluronic chains was able to leave the micelles despite the fact that the structure of the micelles has relaxed multiple times. In Figure 8 we show volume fractions of EO and PO monomers (φEO and φPO) as a function of distance from the micelle core center-of-mass averaged over the entire CGIS simulation for both micelles. To calculate these profiles we used the AES simulations of pure water and ether/water
604 J. Chem. Theory Comput., Vol. 2, No. 3, 2006
Bedrov et al.
Figure 9. A representative snapshot of the L64 micelle in water using the AES representation reversed mapped from the CGIS configuration. Hydrogen atoms are not shown.
Figure 8. Volume fraction profiles of EO and PO monomers and water for L64 and F127 micelles obtained from CGIS LD simulations.
solutions to determine specific volumes occupied by a bulk water molecule (VW)29.8 Å3) and by EO and PO monomers in aqueous solutions (VEO)62.3 Å3 and VPO)105.7 Å3). Then, using configurations obtained from CGIS simulations the local number density of EO and PO monomers (and hence their volume fractions) have been determined, while the local volume fraction of water, φw, has been calculated as φw ) 1 - (φEO+φPO) and is also shown in Figure 8. For both micelles we observe a noticeable fraction of water in the micellar core. Also shown in Figure 8 are values for the core (R1) and corona (R2) radii reported for L64 and F127 micelles of similar sizes by fitting the static structure factors from SANS measurements using a core-corona form factor model and hard-sphere approximation for the micelle interactions.35,36 The best way to make a direct comparison between experiments and our CGIS simulations would be to calculate the static structure factor for micellar solutions at finite concentrations and use the same models and fitting protocol as in interpretation of SANS measurements. These efforts are currently underway. In addition to volume fraction profiles the CGIS simulations can provide detailed information on configuration of polymer chains, fluctuations of micellar shape, and intermicellar interactions/correlations. B. Reverse Mapping onto Atomistic, Explicit Solvent Model. The representative configurations of Pluronic micelle obtained from CGIS simulations can be reverse mapped back to the AES model for further, more detailed investigation of
structural, conformational, and short-time dynamical properties. In the reverse mapping procedure we begin with placing the atomistic Pluronic chains along the contour of the coarsegrained chains. To do this, we first generate atomistic Pluronic chains in random configurations randomly placed in a simulation cell. Then, for the center-of-mass of each ether monomer in the atomistic model a corresponding monomer in the CGIS model with a defined position in space is assigned. An external field is then applied that forces the center-of-mass of the atomistic monomer to occupy the same position as the corresponding monomer in the CGIS model. A short atomistic, implicit solvent LD simulation in the presence of these external fields in which all intermolecular interactions are turned off (phantom chains) allows all atomistic chains to be placed along the contours of corresponding CGIS chains for a given micelle configuration. Since chains arrangement in the configuration from the CGIS simulation already takes into account the excluded volume interactions between ether monomers, the mapped atomistic configuration has only minor overlap between atoms. In the next step, the water molecules are inserted around and inside the micelle proportionally to the water density profile (shown in Figure 8). The unphysical overlap between water and ether atoms is then removed by a subsequent short (several picoseconds) AES MD simulation with upper bounds on nonbonded forces. Figure 9 shows a snapshot of the AES system of an L64 micelle reversed mapped from the configuration obtained from CGIS simulations. The AES system contains 30 000 water molecules. Following the reverse mapping procedure, we are currently conducting extensive AES MD simulations of single micelles hydrated by water. Because the AES system is extremely large the observation times will be limited to a few nanoseconds which is sufficient for redistribution of water molecules inside the micelles and accurate sampling of chain local conformations, characteristics of hydrogen bonding, and
Multiscale Modeling Approach
other atomistic scale properties in the core and corona of the Pluronic micelle. If the micellar structure obtained from AES MD simulations using the reverse-mapped configuration significantly changes/differs from predictions using the CGIS model, the latter will be further revised based on the new structural and conformational data obtained from AES MD simulations and the cycle will be repeated. C. Further Extension of Multiscale Modeling of Pluronic Micelles. The CGIS simulations of preformed individual micelles or ensemble of micelles can provide equilibrium sampling of structural and conformational properties of those micelles; however, they can hardly access the time scales for the second process of micellization which establishes the equilibrium micelle size distribution in the system. While at some thermodynamic conditions it is possible to observe the individual events of Pluronic chains leaving/ joining the micelle or micelle breaking up into two smaller ones, the CGIS simulations are not able to effectively access the milliseconds time scales which would allow the sufficient number of those events to adequately sample the equilibrium micelle size distribution. Therefore, the multiscale modeling protocol described above has to be further enhanced to effectively access larger length and longer time scales and to allow the prediction of the equilibrium distribution of micelle sizes as a function of Pluronic concentration and temperature. One of the approaches which we are currently investigating is the Parallel Tempering method,37 modification of which we have recently developed and successfully applied for equilibration of model self-assembling polymer networks.38 Other potential extensions of the multiscale modeling of Pluronic micelles include coupling with singlechain mean field approach,39,40 density functional theory,17,18 and self-consistent theory approaches.15 However, for all these methods it is important to establish the reliable method to incorporate implicitly all important atomistic and molecular scale physics and phenomena. Acknowledgment. The authors would like to acknowledge the support through NSF ITR, CHE0312226, and MRSEC (University of Colorado), DMR-9809555, grants. Note Added after ASAP Publication. This article was released ASAP on March 28, 2006, with the incorrect Received Date. The correct version was posted on April 12, 2006.
J. Chem. Theory Comput., Vol. 2, No. 3, 2006 605 (8) Almgren, M.; Brown, W.; Hvidt, S. Colloid Polym. Sci. 1995, 273, 2. (9) Pandit, N. K.; McIntyre, H. J. Pharm. DeV. Technol. 1997, 2, 181. (10) See, A. S.; Gast, A. P.; Butun, V.; Armes, S. P. Macromolecules 1999, 32, 4302. (11) Wanka, G.; Hoffman, H.; Ulbricht, W. Macromolecules 1994, 27, 4145. (12) Goldmints, I.; Holzwarth, J. F.; Smith, K. A.; Hatton, T. A. Langmuir 1997, 13, 6130. (13) Hositza, M. J.; Bohne, C.; Alexandridis, P.; Hatton, T. A.; Holzwarth, J. F. Langmuir 1999, 15, 322. (14) Linse, P. Macromolecules 1994, 27, 6404. (15) Hurter, P. M.; Scheutjens, J. M. H. M.; Hatton, T. A. Macromolecules 1993, 26, 5592. (16) Svensson, M.; Linse, P. Macromolecules 1998, 31, 1427. (17) Lam, Y. M.; Goldbeck-Wood, G. Polymer 2003, 44, 3593. (18) Fraaije, J. G. E. M.; Zvelindovsky, A. V.; Sevink, G. J. A. Mol. Sim. 2004, 30, 225. (19) Guo, S. L.; Hou, T. J.; Xu, X. J. J. Phys. Chem. B 2002, 106, 11397. (20) Smith, G. D.; Borodin, O.; Bedrov, D. J. Comput. Chem. 2002, 23, 1480. (21) Smith, G. D.; Borodin, O.; Bedrov, D. J. Phys. Chem. A 1998, 102, 10318. (22) Bedrov, D.; Pekny, M.; Smith, G. D. J. Phys. Chem. B 1998, 10, 996. (23) Smith, G. D.; Borodin, O. In Molecular Simulation Methods for Predicting Polymer Properties; Galiatsatos, V., Ed.; John Wiley & Sons: New York, 2005. (24) Londono J. D.; Annis, B. K.; Habenschuss, A.; Smith, G. D.; Borodin, O.; Tso, C.; Hsieh, E. T.; Soper, A. K. J. Chem. Phys. 1999, 110, 8786. Annis, B. K.; Borodin, O.; Smith, G. D.; Benmore, C. G.; Soper, A. K.; Londono J. D. J. Chem. Phys. 2001, 115, 10998. Borodin, O.; Douglas, R.; Smith, G. D.; Trouw, F.; Petrucci, S. J. Phys. Chem. B 2003, 107, 6813. (25) Smith, G. D.; Bedrov, D.; Borodin, O. J. Am. Chem. Soc. 2000, 122, 9548. Smith, G. D.; Bedrov, D.; Borodin, O. Phys. ReV. Lett. 2000, 85, 5583. Smith, G. D.; Bedrov, D. Macromolecules 2002, 35, 5712. Smith, G. D.; Bedrov D. J. Phys. Chem. B 2003, 107, 3095.
(1) McPherson, A. Methods Enzymol. 1985, 114, 112.
(26) Borodin, O.; Bedrov, D.; Smith, G. D. Macromolecules 2002, 35, 2410. Borodin, O.; Bedrov, D.; Smith, G. D. J. Phys. Chem. B 2002, 106, 5194. Borodin, O.; Trow, F.; Bedrov, D.; Smith, G. D. J. Phys. Chem. B 2002, 106, 5184.
(2) Cudney, B.; Patel, S. Acta Crystallogr. 1994, D50, 414.
(27) http://lucretius.mse.utah.edu.
(3) Harris, J. M. In Poly(Ethylene Glycol) Chemistry: Biotechnical and Biomedical Applications; Harris, J. M., Ed.; Plenum Press: New York, 1992, p 1.
(28) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. J. Chem. Phys. 1983, 79, 926.
References
(4) Andrade, J. D.; Hlady, V.; Jeon, S. I. Polym. Mater.: Sci. Eng. 1993, 69, 60.
(29) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. J. Comput. Phys. 1977, 23, 327.
(5) Sakai, T.; Alexandridis P. Nanotechnology 2005, 16, S344.
(30) Martyna, G. J.; Tuckerman, M.; Berne, B. J.; Martyna, G. J. J. Chem. Phys. 1992, 97, 1990.
(6) Alexandridis, P. Curr. Opin. Colloid Interface Sci. 1997, 2, 478.
(31) Allen, M. P.; Tildesley, T. J. Computer Simulations of Liquids; Oxford University Press: Oxford, 1987.
(7) Adams, M. L.; Lavasanifar, A.; Kwon, G. S. J. Pharm. Sci. 2003, 92, 1343.
(32) Hamley, I. W.; Castelletto, V.; Yang, Z.; Price, C.; Booth, C. Macromolecules 2001, 34, 4079.
606 J. Chem. Theory Comput., Vol. 2, No. 3, 2006 (33) Reith, D.; Putz, M.; Muller-Plathe, F. J. Comput. Chem. 2003, 24, 1624. (34) Wu, C.; Chu, B.; Schneider, D. K. J. Phys. Chem. 1995, 99, 5094. (35) Yang, L.; Alexandridis P.; Steytler, D. C.; Kositza, M. J.; Holzwarth, J. F. Langmuir 2000, 16, 8555. (36) Thurn, T.; Coulderc, S.; Sidhu, J.; Bloor, D. M.; Penfold, J.; Holzwarth, J. F.; Wyn-Jones, E. Langmuir 2002, 18, 9267.
Bedrov et al. (37) Yan, Q.; de Pablo, J. J. J. Chem. Phys. 1999, 111, 9509. (38) Ayyagari, C.; Bedrov, D.; Smith, G. D. J. Chem. Phys. 2005, 123, 124912. (39) Muller, M.; Smith, G. D. J. Polym Sci., Part B: Polym. Phys. 2005, 43, 934. (40) Szleifer, I. J. Chem. Phys. 1990, 92, 6940. CT050334K
