Rationalizing the Phase Behavior of Triblock Copolymers through Experiments and Molecular Simulations

: In this paper, we develop a new coarse-grained model, under the MARTINI framework, for Pluronic block copolymers that is able to describe the self-assembly mechanism and reproduce experimental micelle sizes and shapes. Previous MARTINI-type Pluronic models were unable to produce realistic micelles in aqueous solution, and thus our model represents a marked improvement over existing approaches. We then applied this model to understand the e ﬀ ects of polymer structure on the cloud point temperature measured experimentally for a series of Pluronics, including both normal and reverse copolymers. It was observed that high polyoxypropylene glycol content leads to dominant hydrophobic interactions and a lower cloud point temperature, while high hydrophilic polyoxyethylene glycol content shields the micelles against aggregation and hence leads to a higher cloud point temperature. As the concentration increases, the e ﬀ ect of polymer architecture (normal versus reverse) starts to dominate, with reverse Pluronics showing a lower cloud point temperature. This was shown to be due to the increased formation of cross-links between neighboring micelles in these systems, which promote micelle aggregation. Our results shed new light on these fascinating systems and open the door to increased control of their thermal responsive behavior.


INTRODUCTION
The extraction and purification of valuable compounds from biomass or fermentation media are complex affairs requiring a series of economically and environmentally demanding processing steps often based on potentially hazardous organic solvents. 1The growing legislative and popular pressure has driven the redevelopment of existing processes toward more sustainable alternatives.In line with this current evolution, there is a parallel need for the development of "greener" solvents with low toxicity and cost, which provide measurable advantages through their use.
−4 Pluronics are constituted by sequences of polyoxyethylene glycol (PEG) and polyoxypropylene glycol (PPG) units (often also called poly(ethylene oxide), PEO, and poly(propylene oxide), PPO) in which the hydrophilicity of the LPs can be adjusted by varying the PEG to PPG unit ratio. 5hese triblock copolymers are biodegradable, relatively cheap, biocompatible, and present very low volatility, making them excellent cosolvents for biological and other applications. 6In fact, there are several Pluronics already approved by the Food and Drug Administration (FDA) being applied in personal care formulations, food additives, and medical use. 7The Pluronic family is divided in two major groups, normal (L) and reverse (R), based on the relative positions of the PEG and PPG units.L-Pluronic is characterized by a PEG−PPG−PEG sequence, while R-Pluronic is composed of PPG−PEG−PPG. 5−10 The ability to tune the thermoresponsiveness of the system by varying the molecular weight of the LP as well as the PEG to PPG ratio and its respective sequence confers to these compounds a design flexibility not commonly encountered in simple LPs. 5 Despite significant advances, detailed knowledge of the molecular scale interactions driving the self-assembly of Pluronics is still lacking, and this is of great importance to design and expand the range of thermal responsive LP systems.In this context, computer simulations can address the current lack of experimental data and provide a platform for the inexpensive screening of potential systems.−27 Coarse-grained (CG), 3,19,23,27−37 mean-field coarse-grained (MF-CG), 38 or dissipative particle dynamics (DPD) 39−42 simulations were also deployed to investigate mesoscale structures of these polymers.
The time scale and size limits inherent to quantum density functional theory (DFT) calculations and atomistic MD simulations limit the application of such approaches to systems beyond the initial stages of micelle formation for a few molecules (DFT) or short simulation times (tens of nanoseconds) in AA MD simulations. 32,43Attempts to overcome this limitation by using preassembled structures only provide a qualitative resolution lacking in thermophysical significance as pointed out by Shelley and co-workers. 44In contrast, the CG approach overcomes these limitations by reducing the number of interaction centers involved in the system.In this manner, several atoms can be mapped as one interaction center, speeding up the dynamics and allowing the study of much larger systems as well as increasing the simulation time by at least one order of magnitude.Nevertheless, despite the ability of the CG approach to qualitatively map the entire phase behavior of triblock copolymers 36 or to form complex structures in diblock copolymers such as bilayers or vesicles, 39 they have been hindered by the lack of transferability of the models. 19The development of the MARTINI model 34,45 provides a generally applicable CG framework in which the molecules can be mapped based on an energy matrix of interactions with 18 different bead types. 45The CG models based on the MARTINI force field were able to address the complexity of the self-assembly process of many amphiphilic molecules 46−51 as well as copolymer solutions. 27,28,30,34,35,37It is therefore an ideal framework to study the behavior of triblock copolymers in combination with different solvents and ionic species, in which we are ultimately interested.
In this work, a simple and transferrable CG model for classical MD simulations is proposed and validated with experimental results to address the micelle formation and phase behavior of triblock Pluronic surfactants around the cloud point temperature.The main aim is to find a set of two CG MARTINI beads through which any Pluronic molecule can be built by assigning a single bead to each PEG and PPG monomer unit.An essential prerequisite of the model is that it is able to realistically describe the self-assembly process and the micelle shapes and sizes found experimentally in these systems.For this purpose, six archetypal Pluronics used in liquid−liquid extraction of biomolecules were studied bearing in mind their cloud points and thermoresponsiveness. Computational results were further rationalized by comparison with obtained experimental results, providing an explanation for the different behaviors observed in the studied Pluronic systems.This work provides a new platform for the study of current and future responsive systems of industrial relevance for the development of sustainable and integrated systems mainly focused on the micellar phase behavior of dilute Pluronic aqueous solutions.

METHODOLOGY
2.1.Experimental Section.2.1.1.Materials.The copolymers used in this work as nonionic polymeric surfactants were Pluronic L-31, L-35, L-61, 17R4, 31R1, and 10R5, all of them acquired from Sigma-Aldrich.The general structure of normal and reverse Pluronics, as well as details of specific systems studied here, is presented in Figure 1.
2.1.2.Measurement of the Copolymers Cloud Point.The cloud point determination of 1, 2.5, or 10 wt % aqueous solutions of copolymer was carried out through the visual determination of the onset of turbidity in a solution heated in a temperature-controlled water bath with an error of ±0.01 °C following a methodology described elsewhere. 9,10Each copolymer solution was both prepared in distilled water and in an aqueous solution of 0.18 M McIlvaine buffer (pH 7.0) up to a final volume of 10 mL.The binodal curves were established in McIlvaine buffer solutions since this is the medium of preference to work with pH-sensitive or labile biomolecules.For each system, three replicas were determined, and the average value and standard deviation (always <1 °C) were reported.
2.2.Molecular Dynamics Modeling.The CG MD simulations were carried out with the Gromacs 5.1.4package 52 using a leapfrog algorithm 53 to integrate the equations of motion with a time step of 10 fs.The potential energy function comprised bond stretching, angle bending, and dihedral torsion terms for bonded interactions.The Lennard-Jones (LJ) potential and the Coulomb term were considered for nonbonded interactions.A potential cutoff with the forceswitch modifier function was used in LJ with a cutoff radius of 1.2 nm where the energy decays smoothly to zero from 1.0 to 1.2 nm.The temperature was fixed with the velocity-rescaling thermostat 54 in the equilibration stage and the Nose-Hoover thermostat 55 in the production runs, to ensure adequate ensemble sampling.The pressure coupling was considered as isotropic, and the pressure was fixed at 1 bar using the Berendsen pressure coupling 56 in the equilibration stage and the Parrinello−Rahman barostat 57 in the production runs.Bond lengths were constrained by the LINCS algorithm. 58oulombic interaction energies beyond the cutoff were computed with the Potential-Shift-Verlet coulomb modifier.All simulated systems in this study were enclosed in cubic boxes with periodic boundary conditions with the molecules arranged randomly in the initial configurations.Simulation outputs were visualized using the VMD software package. 59All the simulations carried out in this work are summarized in Table 1, and the same protocol was followed: an energy minimization step using the steepest descent algorithm to prevent short-range contacts between atoms prior to two short equilibrium steps of 0.5 and 2 ns in the NVT and NpT ensembles, respectively.Unless otherwise stated, the duration of the NpT production run was 1 μs.It must be noticed that with the 4:1 or 3:1 mapping used in the CG MARTINI model, a correction factor of about four should be applied to yield a realistic time scale. 45Nevertheless, as we are mostly interested in equilibrium phenomena, this was not applied here; i.e., the time reported in our results is simply the number of steps multiplied by the nominal time step.
As previously mentioned, each PEG and PPG group was represented by a single interaction bead, with the mapping scheme depicted schematically at the top of Figure 1.We started off by assessing the performance of several existing MARTINI CG models for Pluronic triblock copolymer systems. 27,28,30,34,35,37Unfortunately, as discussed in detail in the Supporting Information (SI), none of these models was able to even qualitatively describe the self-assembly and micelle structure of real Pluronics.The Hezaveh et al. 27 model predicted unrealistic bundle-like aggregates with stretched, fully solvated, chains (Figure S1b), displaying no differences between reverse (10R5) and normal (L-35) Pluronics.The model of Nawaz et al. 37,60 was found to have PPG beads that were too hydrophilic, leading to no significant micelle formation (Figure S1d).Conversely, the model of Hakateyama et al. 35 had PEG beads that were too hydrophobic, leading to desolvation of PEG chains and incipient phase separation (Figure S1c).Since the ability to reproduce the self-assembly of Pluronics and the micelle shapes and sizes was central to this study, we developed a new model consisting of "small" MARTINI beadsa size of 0.43 nm with the self-interaction energy scaled by 75% and 3:1 mapping.We chose the intermediate apolar character of SC3 bead types to reproduce the PPG segment, whereas the lowest polar SP1 bead type was selected to represent the polar PEG segments, following the initial attempt of Hezaveh et al. 27 The bonded interactions between PPG and PEG beads were also taken from Hezaveh et al. 27 The water molecules were represented by a 4:1 mapping of P4 beads, following the recommended protocolregular size with a diameter of 0.47 nmand an addition of 10% of antifreeze BP4 beads to avoid the unrealistic tendency of standard CG water beads to freeze. 45It has been demonstrated that the addition of antifreeze beads maintains the essential physicochemical characteristics of water, reproducing, with a reasonable accuracy, the phase behavior of surfactants in aqueous solutions in a wide range of surfactant concentrations and temperatures. 51The new model yielded pure liquid and solvation properties of monomers and dimers of both PEG and PPG in reasonable agreement with experiment, with the exception of the density of PPG, which is somewhat overestimated (see Table S4).More importantly, the model was able to qualitatively reproduce the aggregation behavior and the structure of both normal and reverse Pluronic micelles (Figure S1a).Furthermore, we obtained a reasonable quantitative agreement between the simulated average micelle size and experimental DLS results for a variety of Pluronic systems (Table S3).Details about the model validation, including full tables of parameters (Tables S1 and S2) can be consulted in the Supporting Information.
The individual micelles present in each system obtained in the MD simulations were characterized in detail using a specially developed code 61 based on the Hoshen−Kopelman cluster counting algorithm. 62The criteria followed was that two individual Pluronics chains belong to the same cluster when the distance between the central PPG beads in the chain is lower than 1.2 nm.Tests using reasonable variations of this criteria (e.g., shorter contact distances; different tail beads) showed that the results were quite robust to this choice.The code allows the determination of the average micelle aggregation number as well as the radial density profiles from the center of mass of the micelle using the trajectory files All systems contain 300 Pluronic chains and W is the number of CG water beads (each corresponding to four real water molecules).The temperature was fixed to the experimental cloud point temperature values obtained in the laboratory.
of each MD simulation.The micelle radius was estimated as the distance between the micelle center of mass and the position of the maximum in the PEG density profile (i.e., the average position of the head groups).The diffusion coefficients of Pluronic chains included in the micelles were obtained with the mean square displacement gmx msd gromacs tool.
It is well-known that different configurations of PPG and PEG segments in Pluronics produce diverse micelle sizes in aqueous solutions.Bearing this issue in mind, after a careful micelle distribution analysis over all Pluronics, the MD results can be distributed in three micelle size families of <10, 10−15, and >20 chains per micelle.The majority of micelles, for all Pluronics, fall within the range of 8 ± 3 chains per micelle, as can be seen in Table S5 and Figure S2, and therefore this size range was taken as the most representative to compare the Pluronic micelle density profiles.The aggregation number (number of Pluronic chains per micelle), N a , was obtained with the above-mentioned cluster counting code.Furthermore, a normalized aggregation number N a *, defined as the number of PPG segments (N a * = N a × no.PPG units), was introduced for a fairer comparison of Pluronics with diverse molecular weights.

RESULTS AND DISCUSSION
3.1.Copolymer Self-Assembly.PEG/PPG block copolymers display a lower critical solution temperature (LCST) in water.As the temperature increases, water becomes a worse solvent for these compounds, leading to their self-aggregation and clouding phenomena as well as the formation of mesophases.Upon self-aggregation, the PPG units form the micelle hydrophobic core, while the PEG units form the micelle hydrophilic corona and establish hydrogen bonds with water.These block copolymers form mostly spherical micelles above their critical micelle concentration (CMC) and critical micelle temperature (CMT), which depend on the PPG and PEG molecular weight and their sequence in the block copolymer. 5,63In copolymers, temperature has an even greater effect on the aggregation than in conventional nonionic surfactants due to the impact of PPG and PEG dehydration.
The micellization of block copolymers is an entropy-driven and endothermic process, 64,65 favoring the aggregation with an increase in temperature through two complementary effects: (i) the increase of the hydrophobic effect of the PPG and (ii) the lower solubility of the PEG segments as a result of dehydration.−67 Copolymer micellization is further affected by the presence of cosolutes, such as salts and impurities, cosolvents, and cosurfactants.Depending on their nature, they can either favor or delay the self-assembly process. 5s mentioned in section 2.2, our model was able to describe the self-assembly process of Pluronic surfactants in molecularlevel detail.In Figure 2, we follow the micelle evolution as a function of simulation time, showing that equilibration proceeds in five stages.At the start of the simulation, it is possible to see the isolated Pluronic 10R5 and L-35 chains well distributed along the simulation box in Figure 2, panels i and vi, respectively.In a second stage, at 3 ns of simulation time, in Figure 2ii, the PPG segments of each 10R5 chain join together to form isolated loops, whereas in L-35 the PEG segments remain stretched while the PPG segment coils to avoid contact with water, thus forming relatively straight "knots" as shown in Figure 2vii.In a third stage, after 8 ns, two or three isolated Pluronic chains merged to form the initial micelles.In 10R5, Figure 2iii shows the initial cross-linking process between small micelles, connected by the PEG branches colored in orange.Conversely, in Figure 2viii, no such cross-links are able to form between the L-35 proto-micelles.In the fourth stage, at several tens of nanoseconds of simulation time, the small 10R5 and L-35 micelles (Figure 2, panels iv and ix, respectively) begin to fuse to form aggregates of ∼7 Pluronic chains each, i.e., close to the final micelle size, combined with some isolated Pluronics and other small micelles.The fifth stage shows the equilibrium of micelle distribution obtained in both systems after 100 ns, with flower-like micelles in 10R5, some of them physically connected as shown in Figure 2v, and star-like micelles in L-35, Figure 2x.Notwithstanding the diverse micelle size distribution obtained in Pluronics, the self-assembly process is generally similar to that observed in ionic surfactants. 61However, it is important to notice that the spread in micelle sizes is quite large, as evident in Figure 2v,x.This is in agreement with the large polydispersity experimentally observed in Pluronics. 65,68,69he average micelle sizes obtained in our simulations are in reasonable agreement with experimental DLS measurements (Table S5), particularly taking into account that the high degree of polydispersity leads to large uncertainties in those values.Our model is also able to qualitatively represent the differences in micelle structure between normal and reverse Pluronics.The reverse Pluronics form "flower-like" micelles, named after the tendency of the PEG segments to form loops resembling flower petals. 8Conversely, the PEG peripheral units in the normal Pluronics point outward from the micelle core, arranging in a "star-like" shape (see snapshots in Figure 1).This gives us confidence that our model is providing a realistic description of these complex systems, allowing us to interpret the experimental cloud point measurements.
3.2.Copolymer Cloud Points.As copolymers display a clouding phenomenon, the cloud point is an excellent property to characterize the thermoresponsiveness properties of aqueous solutions of Pluronics.Low molecular weight PPG homopolymers display a clouding behavior at very low temperatures, whereas PEG homopolymers require very high temperatures.By presenting both PPG and PEG units, triblock copolymers, such as Pluronics, phase separate at a wide range of temperatures depending mainly on the number and molecular weight of the PEG and PPG units.With an LCST behavior, Pluronic systems change from a monophasic to a biphasic system above the cloud point. 5,63The standard cloud point for each copolymer was measured both in distilled water and in 0.18 M of McIlvaine buffer pH 7.0, as shown in Figure 3, for a comparison of the behavior of the various copolymers studied here.The cloud point of copolymers in the presence of McIlvaine buffer, composed of the salting-out hydrogen phosphate and citrate anions, was determined to (i) demonstrate the effect of salt addition on the cloud point and (ii) to extend the potential of the studied system to biological applications.In addition to the mixture point of 1 wt % of copolymer, 10 wt % of each Pluronic was also studied in distilled water and McIlvaine buffer to ascertain the effect of polymer concentration on its cloud point.The cloud point of Pluronic L-61 was measured only for concentrations up to 2.5 wt % since it is not soluble at higher concentrations.
Regarding the cloud points found for 1 wt % of the reverse copolymers in water, namely, Pluronics 31R1, 17R4 and 10R5, the obtained values are in good agreement with those previously reported. 5,70The cloud points obtained for 1 wt % of the normal Pluronics L-31, L-35, and L-61 in water are higher (3 to 10 °C higher) than those previously reported. 5he main reason for this difference may be related with the copolymer composition and purity, since the molecular weight is an average parameter over samples with potentially different polydispersity, and the triblock copolymers often present contaminations with monomers and diblock copolymers. 4hen water is replaced by McIlvaine buffer, all the cloud points decrease between 4 and 17 °C, depending on the Pluronic used, as expected when a salting-out agent is added to the system. 63This cloud point difference is kept constant when the copolymer concentration is increased, indicating that independently of the studied solvent, namely, water or buffer, the system displays an identical behavior, varying merely in the temperature range.
3.3.Binodal Curves of Pluronic Copolymers in Aqueous Solutions.Two determinant factors seem to influence the temperature responsiveness of the studied LPs, namely, the PEG/PPG ratio, i.e., the relative hydrophobicity of the Pluronic, and micelle surface effects.The latter is largely dependent on the Pluronic nature, whether it is normal or reverse.To further investigate these issues, the binodal curves of six Pluronic copolymers were determined in McIlvaine buffer pH 7.0 and presented in Figure 4.The Pluronics were selected to evaluate the influence of different copolymer features over the cloud points, for example: (i) copolymers  with the same amount of PEG but different amounts of PPG (L-31 and L-61); (ii) copolymers with the same amount of PPG but different amounts of PEG (L-31 and L-35); (iii) copolymers with gradually increasing PEG/PPG ratio (e.g., 10R5, 17R4 and 31R1); (iv) a normal and a reverse copolymer with the same PEG content and molecular weight (10R5 and L-35).These effects will be analyzed in turn below, using MD simulations at each corresponding cloud point temperature to gain a microscopic insight into the governing phenomena for the thermomorphic transition of Pluronic copolymers.
3.3.1.Effect of PEG and PPG Content.As discussed in section 3.1, Pluronics tend to form nearly spherical micelles in aqueous solutions at dilute concentrations.These micelles are composed of a hydrophobic PPG core and a hydrophilic PEG shell, or corona.Therefore, an increase in the PPG content of the copolymer will tend to produce micelles with a larger hydrophobic core.The experimental measurements show that, all else being equal, a larger micelle core leads to a decrease of the cloud point temperature over the whole range of concentrations.Comparing L-31 with L-61, both have the same amount of PEG, but the latter has almost twice as much PPG.As can be seen in Figure 4, the curve for L-61 is displaced almost uniformly to lower temperatures by about 12 °C.This decrease is accompanied by an increase in the L-31 and L-61 micelle core size from 2.8 to 3.1 nm, while the micelle shape is maintained, as shown in the density profiles and snapshots of Figure 5.A larger micellar core will lead to stronger hydrophobic interactions between micelles, facilitating micelle aggregation and phase separationhence a lower cloud point temperature.
It must be noticed that the density of the PPG core is significantly higher than the experimental density of polypropylene glycol (approximately 1033 kg/m 3 ). 71,72To test whether this was a consequence of the aggregation process, we carried out simulations with 30 pure PPG or PEG chains, and the densities obtained were 1760 and 1807 kg•m −3 , respectively, similar to those found in the micelle density profiles.This means that the density overestimations are an inherent limitation of our CG model, probably arising from the combination of the "small" MARTINI beads with the bonded parameters of Hezaveh et al. 27 Although PEG models using the MARTINI framework exist and provide better agreement with experimental properties of pure homopolymers, 73,74  they have not yet been applied to simulate Pluronic self-assembly.Furthermore, to the best of our knowledge, no such model exists for PPG homopolymers.It is important to note, as discussed in detail in the Supporting Information, that the density overestimation does not affect the micelle formation process or the micelle structure and shape, leading to micelle size distributions in good agreement with experimental data.Therefore, this limitation does not affect the qualitative conclusions of our study.Nevertheless, we intend to address this shortcoming in future work.
It is interesting to note that the effect of PPG content is, at least qualitatively, rather independent of the polymer architecture (normal or reverse).For instance, 10R5 and 17R4 have almost the same PEG content (in this case, there is a small difference of 2 units), but the latter has much more PPG.Again, this leads to an increase in the micelle core (see profiles and snapshots in Figure S3) and a corresponding decrease in the cloud point temperature.However, in this case the decrease is much more pronounced (∼20−22 °C).As we will discuss in more detail below, this difference is most likely due to the much larger corona size (i.e., PEG content) in the 10R5 and 17R4 copolymers, relative to L-31 and L-61.
The PEG content has an important impact on the cloud point due to the increase in hydrophilicity of the micelle surface.Indeed, it has been proposed that the main driving force for the existence of a LCST in PEG-containing surfactants is a conformational change in the PEG chains with increasing temperature, leading to weaker interactions with water and consequently phase separation. 75The effect of PEG content can be clearly seen from our data by comparing Pluronics with similar amounts of PPG but different number of PEG units.For example, L-35 has the same PPG content as L-31 but more than 5 times as much PEG.This increase in PEG content causes a massive increase in the cloud point temperature, by about 30 °C at low concentration and more than 40 °C at the highest concentration reported in Figure 4. Again, this can be rationalized by considering the effect on intermicelle interactions.The shell, or corona, is what ensures that the micelles are stable against aggregation in aqueous solution, as these hydrophilic PEG moieties are soluble in water.It is clear from the density profiles and snapshots shown in Figure 6 that the PEG "arms" assume a rather extended configuration at the micelle surface, aiming to maximize interactions with water.Naturally, the larger the diameter of  the PEG shell, the higher the degree of shielding of the micelle core.Conversely, the much smaller shell of L-31 allows micelle cores to come into relatively close proximity, where hydrophobic interactions take over, and this consequently leads to an earlier (i.e., lower T) onset of phase separation.
A similar effect can be observed by comparing L-61 with 17R4 (both have similar PPG content, but the latter has a much larger PEG group), where an increase of about 18 °C in cloud point temperature can be seen.In fact, Figure S4 shows that the L-61 system yielded larger micelles as a consequence of shorter PEG "arms"presumably due to more PPG micelle core contacts enhancing the micelle fusion processes and aggregation.Here, however, it is important to note that we are comparing a normal with a reverse Pluronic, and hence more than one effect is at play in such a case.We will return to the effect of polymer architecture later in the paper.
3.3.2.Effect of the PEG/PPG Ratio.The joint effect of PPG and PEG content can be encapsulated in the PEG/PPG ratio for each copolymer (Figure 4).Bearing this in mind, the binodal curves of the three reverse Pluronics studied here (Figure 7I) follow the trend: 31R1 < 17R4 < 10R5.In 31R1, 17R4, and 10R5 the PEG/PPG ratio is 0.15, 0.86, and 1.37 respectively, suggesting that the cloud point temperature is roughly proportional to this ratio.A similar qualitative behavior was reported in previous experimental studies. 5,63,76t must be remarked that the L-61 and 31R1 systems, which have quite different molecular weights, have very similar cloud point curves.Once again, this is likely due to their similar PEG/PPG ratios (0.13 and 0.15, respectively), which suggests that they share a similar balance between attractive core−core interactions and the protective effect of the PEG shell.
Figure 7II shows the MD simulation snapshots for the three reverse Pluronics at 1 wt % near their respective cloud transition temperatures after 1 μs of simulation time.In all three snapshots, the initial stages of formation of micellar aggregates can be seen, which would eventually lead to phase separation.Unfortunately, even with CG models, we cannot access the necessary large sizes and long times to allow the phase separation to unfold.Nevertheless, it is interesting to see that all three systems seem to be at a similar stage with regard to intermicelle interactions, even though the simulations were run at significantly different temperatures.At these low concentrations, as discussed above, the relative size of the core and shell of the micelles seems to be the determining factor in driving the phase separation.We will analyze the systems at higher copolymer concentration in the next section.
Table 2 shows the PEG/PPG ratio, the diameter (⌀), the aggregation number (N a ), and the standardized aggregation  number N a * obtained in the MD simulations near the cloud point temperature.The results show different micelle sizes according to the distinct molecular weight, with a micelle diameter increase of ∼0.3 nm from 10R5 to 17R4 and ∼0.6 nm from 17R4 to 31R1.The average number of Pluronic chains per micelle, N a , is 5, 6, and 7 in 31R1, 17R4, and 10R5, respectively, but this does not properly reflect the micelle size due to the different PEG/PPG chain sizes in each case.Figure 8 illustrates the cloud point temperature trend with the normalized aggregation number N a * depicting once again a nonlinear behavior.The cloud point temperature seems to reach a plateau for 31R1 or L-61; thus N a * is a good reference to assess the cloud point temperature behavior of Pluronics, at least in dilute aqueous solutions.The micelle diameter increases from 10R5 to 31R1, as shown in the density profiles of Figure 9 and in Table 2.It must be stressed that the low PEG/PPG ratio and the relatively long PPG chains in 31R1 (high PPG core and small PEG micelle surface) the PPG micelle core more accessible to contact with the PPG core of neighboring micelles, enhancing the micelle fusion in the initial stages of micelle formation (Figure 7II).This propensity to aggregate is reflected in the lower cloud point temperature of the 31R1 and the formation of relatively big micelles (N a * = 260).In fact, the 31R1 core showed a lower PPG density compared with the other normal and reverse Pluronics.This could be due to steric restrictions, caused by the relatively small PEG chain in 31R1 having to arrange in a flower-like configuration, somehow affecting the PPG micelle core arrangement.The N a * shows that the PPG core size follows the trend 31R1 > L-61 > 17R4 > L-31 > L-35 > 10R5, which is in agreement with the obtained micelle density profiles.
3.3.3.Micellar Surface Effects on the Thermomorphic Behavior.Figure 10I shows the cloud points obtained in 10R5 and L-35 at different concentrations, the two Pluronics with the highest cloud points in this study.Both systems have the same PEG and PPG content, thus the same relative amphiphilic character (same N a *); however, the cloud point temperature strongly differs as soon as the Pluronic concentration increases above 1 wt %.This effect could be related to the different shapes of the micelle crown composed of PEG units, 77 as illustrated in Figure 1.The results in Figure 10 suggest that the structural variations of the micelle surface play an important role in the polymer behavior at higher concentrations.To gain a detailed microscopic understanding of this phenomenon, six MD simulations were carried out for 10R5 and L-35 copolymers at 1, 5, and 20 wt % concentrations, respectively, using the CG model developed in this work.The temperature was set at the experimental cloud point: 61 °C for 1 wt % in both, whereas above 1 wt %, 50 °C was set for 10R5 and 65 °C for L-35.
Figure 10II shows the results obtained after 1 μs of simulation time with the MD simulation snapshots displaying the final micellar organization.In Figure 10, II-i and iv, both systems at 1 wt % concentration depict a similar micelle configuration homogeneously distributed along the simulation box.However, when the concentration is raised above 1 wt %, intermicelle interactions become dominant, as illustrated in Figure 10II-ii, iii, v, and vi for 5 and 20 wt % Pluronic concentrations.It is clear that at 5 wt %, the reverse Pluronic system is much more compact than the normal Pluronic system, with micelles closely packed together.At even higher concentration, both systems are quite densely packed, but again the compactness is enhanced for 10R5.
In Figure 11, we take a closer look at the structure and evolution of the two systems at 5 wt %.A dynamic picture of this micelle formation can also be seen in the simulation movies MS1 and MS2 for 10R5 and L-35, respectively (see the Supporting Information).From an early stage, it was observed that the 10R5 system was able to form "cross-links" between different aggregates; i.e., the PPG moieties, being on the edges  of the chain, were able to attach to two separate micelles with the PEG moiety acting as a bridge.Conversely, the PPG moieties of the L-35 system, lying at the center of the chain, are only able to attach to a single aggregate.This physical cross-linking was observed to promote micelle fusion during the self-assembly process, by drawing neighboring micelles closer together and facilitating the establishment of core−core interactions.In the L-35 system, on the other hand, micelle fusion took place only through weak PEG−PEG interactions and was thus relatively slower.
To quantify this effect, we analyzed the systems shown in Figure 10II-i-vi using our cluster counting code.The aggregation number N a (N a , cf. Figure S5), the total number of micelles (N), the number of cross-linked micelles (N c-l ), and the averaged micelle diameter (⌀) for 1, 5, and 20 wt % Pluronic concentrations are presented in Table 3.In can be seen that the simulations for 10R5 and L-35 at 1 wt % both yielded an average aggregation number of N a = 7 and similarly sized micelles, as expected.Above this concentration, at 5 wt %, the micelle density profiles remain very similar between the two Pluronics (see Figure S6), but the cluster counting indicates the presence of several 10R5 micelles with 20 chains, whereas L-35 barely increases its micelle size to 12 chains.As explained above, this is due to more micelle−micelle fusion events enhanced by the physical cross-linked micelles.Finally, the N a was increased at 20 wt % in 10R5 and L-35, with 27 and 19 Pluronics per micelle, respectively.
Table 3 also reports the number of physically connected micelles evaluated with the cluster counting code.At 1 wt % of    10R5, there were less than 0.33% of cross-linked micelles.However, as the copolymer concentration increased, the percentage of physically cross-linked micelles in that system was 27% and 55% for concentrations of 5 and 20 wt %, respectively.In the L-35 system, no physically cross-linked micelles were observed, as expected, and only PEG−PEG contacts played a role in intermicelle interactions.The above results demonstrate the ease of 10R5 micelles to aggregate and, more importantly, highlight the role of physical cross-linking between 10R5 micelles, as shown in Figure 10II−v and in more detail in Figure 11a and Figure S7.
A dynamic property such as the coefficient of diffusion (D cd ) of the PPG chains can be used to study the mobility of the Pluronic micelles obtained in the above simulations (Table 3).As a reference, the D cd for isolated 10R5 or L-35 chains in water (representing the infinite dilution limit, with <0.5 wt %) was obtained, and the results were D cd = 0.77 × 10 −5 cm 2 •s −1 and D cd = 0.55× 10 −5 cm 2 •s −1 .As expected, the diffusion coefficient of the PPG moieties decreases when they are incorporated into micellar aggregatesthe corresponding values for 10R5 and L-35 at 1 wt % are D cd = 0.17 × 10 −5 cm 2 •s −1 and D cd = 0.16 × 10 −5 cm 2 •s −1 , respectively (Figure 10II-i and iv).Importantly, both in the infinite dilution limit and in the system with small micelles, the diffusion coefficients for the two Pluronic systems are very similar.However, an increase in the Pluronic concentration to 5 wt % yielded markedly different behavior, with D cd decreasing much more significantly for 10R5 than for L-35.The 10R5 decrease in D cd is caused by the formation of linkages between micelles, as described above (Figure 11), which slow down the dynamics of the system.This behavior was accentuated at 20 wt %, where the D cd = 0.0017 cm 2 •s −1 in 10R5 was over 1 order of magnitude lower than for the L-35 system, D cd = 0.0301 cm 2 • s −1 .
Taken together, our results show that the 10R5 system requires less energy to coalesce and separate into two macroscopic phases due to the presence of physically crosslinked micelles.Consequently, this is reflected in its lower cloud point temperature compared to L-35 (Figure 10I), in which only relatively weak PEG−PEG interactions are possible.

CONCLUSIONS
This work provides a detailed study on the effect of Pluronic copolymer structures on their binodal curves by molecular dynamics simulations and experiments.The main goal of this work was the development of a CG framework for MD simulations of Pluronic aqueous solutions.Experimental results such as the micelle size distribution and micelle shape were taken as a reference to validate the CG framework presented in this work.Different CG models of Pluronics available in the literature were also tested.However, our approach was the only one able to reproduce the experimental micelle distribution (aggregation number and micelle diameter) obtained in our laboratory as well as the expected micelle shapes, with PEG branches forming flower petals loops on the micelle surface for reverse Pluronics and star-shaped architectures for normal Pluronics.Nevertheless, our model overestimates the density of the homopolymer segments, a shortcoming that we intend to address in future work.
Experiments and computer simulations demonstrated that the PPG has an important impact in the cloud point temperature; thus, the higher the PPG content, the lower the cloud point temperature.When the PPG is at a higher proportion compared with PEG, the balance easily shifts toward hydrophobic core−core interactions, promoting aggregation and decreasing the cloud point temperature.In systems with larger PEG coronas, these are able to shield the micelles against aggregation, leading to higher cloud point temperatures.A particularly interesting effect was observed when comparing a normal and a reverse Pluronic with the same hydrophobic/hydrophilic strength.Experiments depicted a very different thermal response at concentrations above 1 wt %, with the reverse Pluronic showing a significantly lower cloud point temperature (by about 10 °C).Our simulation results showed that the micelle architecture plays an important role in controlling the cloud point temperature.In particular, we described in detail the formation of physically cross-linked micelles since the very early stages of micelle formation in the reverse system.This markedly affects the mobility of micelles, corroborated by a decrease in the PPG diffusion coefficient observed in our simulations.By acting as bridges between micelles, these cross-links make it easier for micelles to aggregate and hence lower the cloud point temperature.The formation of cross-links increases with concentration and was only observed above 1 wt %, thus explaining the changes in the binodal curves observed experimentally above this concentration.
The CG simulations performed here enable us to rationalize and control the thermal response of Pluronic aqueous solutions by manipulating the PEG and PPG content, as well as the polymer architecture.This work provides, for the first time, an intuitive computer simulation framework to study the phase behavior of Pluronic aqueous solutions, opening the door for designing tailor-made thermally controlled solvents by computer simulations.

* S Supporting Information
The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acs.jpcc.9b04099.
Detailed description of the model validation, comparing against previous coarse-grained models for MD simulations.Dynamic light scattering experiments carried out in our laboratory to obtain reference micelle sizes for the six Pluronics studied in this work.Additional figures are also shown to facilitate the comprehension of the manuscript (PDF) Simulation movies MS1 and MS2 in MP4 format (ZIP)

■ AUTHOR INFORMATION
Corresponding Authors

Figure 1 .
Figure 1.Pluronic systems used in this study, three normal L-31, L-61, and L-35 and three reverse 31R1, 17R4, and 10R5.The micelle pictures were taken from the MD simulations carried out in this work.The color code is as follows: the apolar PPG is in black, and the polar PEG is in orange.The % PEG is given by the molecular weight of each coarse-grain bead, 58 and 44 UMAS, for PPG and PEG, respectively.

Figure 2 .
Figure 2. 10R5 and L-35 initial self-assembly stages of micelle formation obtained in the MD simulations.Both Pluronics were run at 5 wt % of Pluronic concentration at 298 K.The simulation time in nanoseconds of every snapshot is specified.Water molecules were removed for clarity, and the color code is the same as in Figure 1.

Figure 3 .
Figure 3. Cloud point of aqueous solutions of 1 wt % (circles) and 10 wt % (squares) of each copolymer either in distilled water (blue) or in 0.18 M of McIlvaine buffer pH 7.0 (orange).For L-61, the highest concentration is 2.5 wt % (triangles) and not 10 wt % due to solubility issues.

Figure 5 .
Figure 5. Micelle density profile comparison of L-31 (solid) and L-61 (dashed) at 1 wt % of Pluronic concentration at the cloud point temperatures 35 and 20 °C, respectively.The color code is PPG in black, PEG in orange, and water in blue.

Figure 6 .
Figure 6.Micelle density profile comparison of L-31 (solid) and L-35 (dashed) at 1 wt % of Pluronic concentration at the cloud point temperatures 35 and 61 °C, respectively.The color code is the same as in Figure 5.

Figure 7 .
Figure 7. (I) Experimental binodal curves taken from Figure 3 for 31R1 ( ⧫ ), 17R4 (▲) and 10R5 (•).Red dashed line corresponds to 1 wt % of Pluronic concentration used in the MD simulations.(II) MD simulation snapshots after 1 μs of simulation time at 21, 38, and 61 °C for 31R1, 17R4, and 10R5, respectively.Water molecules in the MD simulation snapshots were removed for clarity.Color code is the same as in Figure 1.

Figure 8 .
Figure 8. Cloud point temperature trend at 1 wt % of Pluronic concentration versus the normalized aggregation number N a * (N a × number of PPG segments) of ( □ ) reverse and (Δ) normal Pluronics.

Figure 10 .
Figure 10.(I) Experimental binodal curves for L-35 (Δ) and 10R5 (•) at different concentrations.(II) Six MD simulations were carried out at 1 (i) and (iv), 5 (ii) and (v), and 20 wt % (iii and vi) concentrations at the cloud point temperature conditions for L-35 (top) and 10R5 (bottom) after 1 μs of simulation time.The diffusion coefficients cm 2 •s -1 are reported below each snapshot.Water molecules in the MD simulation snapshots were removed for clarity.Color code is the same as in Figure 1.

Figure 11 .
Figure 11.Initial stages of micelle formation at 5 ns of simulation time for 10R5 and L-35 systems at 5 wt % concentration.The panels on the right-hand side are blow-ups showing cross-linked micelles in the 10R5 system and loose PEG segments in the L-35 system.Water molecules were removed for clarity.Color code is the same as in Figure 1.

Table 3 .
Diffusion Coefficient (D cd ), Average Aggregation Number (N a ), Total Number of Micelles (N), Number of Cross-Linked Micelles (N c-l ), and Average Micelle Diameter (⌀) Obtained after 1 μs of Simulation Time for L-35 and 10R5 L-35 10R5

Table 1 .
Molecular Dynamics Simulations Carried out in This Work a

Table 2 .
PEG/PPG Ratio, Average Micelle Diameter (⌀), Aggregation Number (N a ), and Normalized Aggregation Number a Obtained for 1 wt % of Pluronic Concentration at the Cloud Point Temperatures between Brackets a N a * = N a × number of PPG segments.