Micellization of Pluronic P123 in Water/Ethanol/Turpentine Oil Mixed Solvents: Hybrid Particle-Field Molecular Dynamic Simulation.

The hybrid particle-field molecular dynamics simulation method (MD-SCF) was applied to study the self-assembly of Pluronic PEO20-PPO70-PEO20 (P123) in water/ethanol/turpentine oil- mixed solvents. In particular, the micellization process of P123 at low concentration (less than 20%) in water/ethanol/turpentine oil-mixed solvents was investigated. The aggregation number, radius of gyration, and radial density profiles were calculated and compared with experimental data to characterize the structures of the micelles self-assembled from P123 in the mixed solvent. This study confirms that the larger-sized micelles are formed in the presence of ethanol, in addition to the turpentine oil-swollen micelles. Furthermore, the spherical micelles and vesicles were both observed in the self-assembly of P123 in the water/ethanol/turpentine oil-mixed solvent. The results of this work aid the understanding of the influence of ethanol and oil on P123 micellization, which will help with the design of effective copolymer-based formulations.

1. Introduction

The self-assembly of amphiphilic block copolymers has attracted considerable attention for many decades because they can form many novel structures, including spheres, cylinders, bicontinuous phases, lamellae, vesicles, and many other complex or hierarchical assemblies [1,2,3,4]. Amphiphilic block copolymers have been widely used in many industrial applications, such as detergents, dispersion stabilizers, foaming agents, emulsifiers, and lubricants [5]. Triblock copolymers composed of hydrophilic (polyethylene oxide) −amphiphilic (polypropylene oxide)− hydrophilic (polyethylene oxide) blocks are commercially known as Pluronics or Poloxamer [6]. By changing the relative composition of the polyethylene oxide (PEO) and polypropylene oxide (PPO) blocks (general formula: (PEO)m−(PPO)n−(PEO)m), it is possible to tune the hydrophilic–lipophilic balance (HLB), which makes the main structural differences between Pluronics. Since the mid-1960s, Pluronics have been widely and successfully used in biomedical applications, due to their biocompatibility and non-toxicity, for example, micellar nanocarriers for drug delivery, biological reaction modifiers in chemistry, non-viral gene therapy agents, and wound dressings for the treatment of thermal burns [7,8,9,10,11,12,13]. Pluronics P123 (PEO20-PPO70-PEO20) is one of the most used block-copolymers for several applications. For these reasons, the micellar assembly of P123 in water has been widely studied by a variety of experimental techniques, such as static light scattering, fluorescence spectroscopy, and surface tension measurements [13,14,15,16]. Furthermore, the effect of mixed solvents on the aggregation of P123 [17] and the phase behavior of P123 in mixed ethanolic solvents has been the subject of several studies [18,19]. Very recently, Chat et al. successfully investigated the aggregation features of P123 in the presence of lavender oil and ethanol using dynamic light scattering, small angle neutron scattering, and rheology [20]. They showed that lavender oil assisted aggregation characteristics of Pluronic P123 in the mixed solvents of ethanol and water [20]. Furthermore, Ganguly et al. studied tea tree essential oil (TTO) on the properties of aqueous solutions of Pluronic P123 and the effect on micellar solubilization. They showed that the solubilization of TTO brings about spherical-to-worm-like micelles-to-vesicular structural transitions in aqueous solutions of Pluronic P123 [21].

Although experimental data can help with the understanding of the micellization mechanism of P123 in mixed solutions, it is hard to obtain good information of the phenomenon at molecular level. Meanwhile, using computer simulation methods to study the self-assembled structure and properties of block copolymers in aqueous solution has gradually attracted scientists’ attention [21,22,23,24]. Computer simulations can not only explain the underlying self-assembly mechanism, but also offer a route for designing and predicting future experiments. In particular, the standard all-atom molecular dynamics (AAMD) method cannot be directly performed to investigate the self-assembly phenomenon because of the time and length scale limits. Therefore, mesoscopic models, such as dissipative particle dynamics (DPD) [25,26,27,28,29,30], MesoDyn [31], and the Monte Carlo method [32,33], which have access to larger time and length scales, have successfully been used to investigate the micellization phenomenon. For example, Huang et al. used DPD to simulate the entire phase diagram of Pluronic L64 and water. Furthermore, they also proved that non-Newtonian behaviors can be predicted by non-equilibrium DPD simulation [25]. Recently, a different computational approach, the hybrid particle–field MD (MD-SCF) method, was developed, in which the MD and self-consistent field (SCF) method are combined [34,35]. This hybrid particle–field MD method not only has similar accessibility to large time and length scales of pure SCF methods, but also has a chemical specificity that is similar to the MD method. In particular, in our previous works, we presented coarse-grained (CG) hybrid particle–field models to study the phase behavior of the mixture of phospholipid-water [36], Triton TX-100 surfactant [37,38], Pluronic-water [39], and nanocomposite systems [40,41]. However, to the best of our knowledge, studying the effect of oil on the aggregation behavior of P123 in the mixed solvents using a computer simulation method is still lacking.

In this paper, we present a study of the aggregation features of Pluronic P123 in the presence of water, ethanol, and turpentine oil, using a MD-SCF simulation. This paper aims to understand the self-assembly mechanism of micellization of the P123 in the mixed solvents at the molecular level. Furthermore, the influence of turpentine oil on the aggregation behavior of P123 is also discussed.

2. Materials and Methods

2.1. Simulation Method

In 2009, a computer simulation method of hybrid particle–field molecular dynamics method was developed [34,36,39]. The main feature of this method is to combine the molecular dynamics method based on particle description with the average field theory of self-consistent field. That is, the non-bonded interaction between particles of different molecules can be replaced by the evaluation for each particle of those with an external potential that depends on the local density at position r, and can solve the multi-body problem [35] using the molecular dynamics simulation method. In the framework of the SCF theory, one chain is regarded to interact with the surrounding chains not directly but through a mean field. Its calculation method is as follows: the Hamilton operator consists with two parts composed of M molecules, which can be expressed as: is the Hamiltonian of a reference ideal system composed of M molecules but with all the bonded and non-bonded interaction in the standard MD simulations. Γ specifies a point in the phase space, and the hat symbol indicates that the relevant physical quantity is a function of the microscopic states described by the phase space Γ. On the other hand, the deviation from the reference system, which is induced by non-bonded interactions between molecules, is accounted for by the term .

Then, the total partition function of the system is written as:where .

From the microscopical level, the density distribution of particles is defined as the sum of functions centered at the center of mass of each particle as: where is defined as the particles’ coordinate of the Pth molecule, is defined as the number of particles contained in the Pth molecule, and M is the total number of molecules in the system. We assume that the particle density of the phase space determines the external field:

Here, we assume that depends on only through the particle density . Then, we can write the single-molecule partition function as:In this expression, Z is partition function of the single molecule, is a conjugate field of , and is the external potential.

Next, we assume that the density dependent interaction potential , where each species is specified by the index , takes the following form in Equation (6): where is the coarse-grained density of the species at position r, and represents the mean field parameters for the interaction of a particle of type with the density fields due to particles of type . In particular, is the compressibility that is assumed to be sufficiently small to obtain the relaxed incompressibility condition. Then, the external potential of type can be given by:

Finally, the density derivatives of Equation (7) are required to calculate the forces, further details can be found in [34]. Moreover, the values of the coarse-grained density at lattice points are not updated at every time step but only at every prefixed density update time. Between two updates, the values of the densities on the lattice are used to interpolate density derivatives and, hence, forces. Details concerning the OCCAM code used to perform the MD-SCF simulations are reported in [35].

In the hybrid particle–field molecular dynamics simulations, the box length of the system was 64.25 nm. The total number of particles in the system was 192,000, and the diameter of the particle was set at 1.046 nm. The number density was 0.725. The time step used for the integration of the equations of motion was set to 0.03 ps. The grid size was set to 2.1 nm, and the density field was updated every 300 time steps (Update frequency has been tested in several of our previous papers and guarantees reasonable results with respect to finer methods [35,36,37,38,39,41].) The NVT ensemble was chosen, the temperature was kept constant at 303K, the same as in the experiments of Reference [20], using the Andersen thermostat with a collision frequency of 20 ps−1. The total simulation step was 8 × 107 steps, which corresponds to the real time 2.4 µs. The detailed information of the simulation systems is reported in Table 1 and Table 2.

Table 1

The composition of systems including block copolymer P123 and water.

VolumeFraction(vol %)	P123No. of Atoms	WaterNo. of Atoms	ParticleNo. of Atoms	Box Length(nm)
2	3834	188,166	192,000	64.25
5	9594	182,406	192,000	64.25
10	19,206	172,794	192,000	64.25
15	28,800	163,200	192,000	64.25
20	38,394	153,606	192,000	64.25

Table 2

The composition of systems including block copolymer P123, water, ethanol, and turpentine oil.

System	P123No. of Atoms	WaterNo. of Atoms	EthanolNo. of Atoms	Turpentine OilNo. of Atoms
I	9576 (5%)	153,624	28,800 (15%)	0
II	9576 (5%)	134,424	48,000 (25%)	0
III	28,800 (15%)	115,200	48,000 (25%)	0
IV	9576 (5%)	150,744	28,800 (15%)	2880 (1.5%)
V	9576 (5%)	131,544	48,000 (25%)	2880 (1.5%)
VI	28,800 (15%)	112,320	48,000 (25%)	2880 (1.5%)

2.2. Coarse-Grained (CG) Model

The coarse-grained model was used in our simulations, in which the mapping scheme was similar to our previous work [29,30]. According to the molecular weight and experimental bulk densities, the monomer volumes were obtained. In particular, we took 20 water molecules and grouped them together into one bead with a volume of 600 Å3. The same bead volume corresponds approximately to 9 PEO monomers, 5 PPO monomers, 6 ethanol monomers or 2 turpentine oil monomers. With the proposed mapping scheme, a single chain of P123 could be represented by EO2PO14EO2. The mapping scheme of the CG model used in this paper is reported in Figure 1.

Figure 1

Mapping scheme of the coarse-grained models of P123 (EO2PO14EO2), water (corresponds to 20 H2O monomers), ethanol (corresponds to 6 ethanol monomers), turpentine oil (corresponds to 2 monomers).

According to our previous work [29], the bead-bead interaction parameters in the MD-SCF simulation method can be obtained from the Flory–Huggins parameters, which can be calculated by Equation (8): where = 600 Å3 in our MD-SCF simulations. The Hildebrand solubility parameter is the square root of the cohesive energy density. The cohesive energy density is the amount of energy needed to completely remove a unit volume of molecules from their neighbors to infinite separation (an ideal gas), which is equal to the heat of vaporization divided by the molar volume.

As PEO is miscible with water at almost any concentrations, was adapted [42,43,44]. In [44,45,46], corresponding to a bead volume of 150 Å3, the authors used the group contribution method to estimate , and . Thus, , which can be deduced based on the bead volume of 600 Å3. According to [47], , and can be deduced from their DPD simulations. Furthermore, according to [48], the solubility parameters of PEO, PPO, and ethanol are 19.9 , 18.5 , and 26.5 , respectively. Moreover, the solubility parameter of turpentine oil is 16.6 , which can be obtained in the polymer handbook [49]. Thus, , , and can be calculated. Then, the interaction parameters of each component are multiplied by RT in the MD-SCF simulations and shown in Table 3.

Table 3

The interaction parameters (kJ·mol−1) between different particles in the MD-SCF simulation method.

	PEO	PPO	Water	Ethanol	Turpentine Oil
PEO	0.0	30.3	0.76	16.0	21.0
PPO	30.3	0.0	21.2	25.4	2.3
Water	0.76	21.2	0.0	0.76	42.3
Ethanol	16.0	25.4	0.76	0.0	36.0
TurpentineOil	21.0	2.3	42.3	36.0	0.0

For the P123 model, the bonds between two successive beads are described by a harmonic potential: where the equilibrium distance is equal to 1.12 nm for the block copolymer. is the bond force constant. The CG model of BCP contains no angular constraints between two consecutive bonds. All the simulations were performed by OCCAM code, which was developed by us [35]. The various analyses were computed from a single ensemble snapshot after the equilibration for each system of 2.4 µs. The calculated diffusion coefficient from the mean square displacement (MSD) of the MD-SCF simulation is for a polymer chain of 20 beads. This value can be compared with experimental (or calculated by atomistic simulations) diffusion data from different polymers that have a similar chain size. Therefore, we can consider an order of magnitude difference in the time scaling factor~1000 from [41]. According to this, the simulations reported in our paper can be considered on the timescale of 2.4 milliseconds.

3. Results and Discussion

3.1. Effect of Concentration on the Micelles Structure of P123 in Water

Firstly, we consider a binary system composed of P123 (with different volume concentrations) in water. For each concentration, a certain number of P123 chains were randomly placed in the box, previously filled with water CG beads (see Table 1). Then, the dilute aqueous solutions were simulated at 303K with a parallel version of the OCCAM code [35].

For all systems, starting from approximately 200 ns, the formation of aggregates of P123 was observed. In particular, the mean aggregation number of Pluronic micelles increased with simulation time up to an equilibrium value (after 2 ), as reported in Figure 2. In order to characterize the composition of the micelles self-assembled from P123 dilute solution, the mean aggregation number was calculated. In particular, was defined as: , where was the fraction of the aggregates whose aggregation number equaled and [50]. Figure 2 shows the time evolution of for micelles formed by the P123 dilute aqueous solution with different concentrations. As can be seen from Figure 2, the average composition of P123 micelles is the function of the concentration of the P123 in water. In general, the P123 micelles contain from 5 to 22 P123 chains.

Figure 2

The time evolution of of micelles self-assembled by P123 dilute aqueous solution with different volume concentrations of 2%, 5%, 10%, 15%, and 20%, respectively.

To describe the morphology of the micelles, we calculated the eigenvalues of the radius of the gyration tensor [50,51,52], as defined by And where a and b indicate x, y, or z components of the bead’s coordinates, respectively. Here n is the total number of beads in the aggregate. acm and bcm are the centers of mass of the aggregate. The radius of gyration R can be simply calculated by Equation (11): where the three eigenvalues can be obtained via diagonalization of the matrix in Equation (11). Finally, the asphericity parameter Q is calculated through the formula as follows: In Figure 3a,b, the mean radius of gyration <> and the asphericity parameter <Q> of the overall micelles of P123 (at the volume concentration of 2%, 5%, 10%, 15%, and 20%) are reported.

Figure 3

The average , , , and of the overall micelle (a), the average asphericity parameter Q of overall micelle (b), R of the largest micelle and R of the largest micellar core (c), and asphericity parameter Q of the largest micelle (d) corresponds to the P123 dilute aqueous solution with the volume concentration of 2%, 5%, 10%, 15%, and 20%, respectively.

Moreover, we concentrated on the largest micelle of each system at the equilibrium state. Thus, the radius of gyration R, the radius of gyration of the micellar core R, the asphericity parameter Q, and the morphology of the largest micelle are shown in Figure 3c,d, respectively.

In general, we observed an increase in the size of the micelles when the P123 concentration increased (see the red curve in Figure 3a), which is in agreement with the experimental behavior [13,15]. In detail, by visual inspection of simulation trajectories, we found that for a low P123 concentration, the micelles were almost spherical, while, for a higher concentration, the aggregates were much more elongated. To quantify the qualitative observation, we computed the asphericity parameter <Q> (in Figure 3b). We observed that for a low P123 concentration (up to 15%), <Q> was close to 0.5, confirming a spherical shape of the P123 aggregates, which is consistent with the experimental observations [15,53] and previous simulation studies [29,39]. This indicates that our coarse-graining procedure and our hybrid MD-SCF simulation method is appropriate for this P123 dilute solution system. From Figure 3d, we can see that Q equals 0.0069 and 0.035 at 5% and 10% P123 dilute aqueous solution, respectively, which means that the largest micelles at these two concentrations are closer to the regular sphere. At higher concentration of 15% and 20%, anisometric micelles can be found with Q values of 0.13 and 0.41, respectively. This can also be found in previous experiments [15], and is due to the larger aggregation number of P123 and larger micellar cores of the micelles.

3.2. Effect of Ethanol on the Micelles Structure of P123 in the Mixed Ethanol-Water Solvent

On the basis of the coarse-graining and mesoscale simulation strategy, we then consider the self-assembly of P123 in the mixed water and ethanol solutions. The main goal of the study was to find the effect of ethanol on the micellization of Pluronic P123. To compare the simulation results with the experimental results of [20], the volume concentration of P123 was chosen by 5 % and the temperature of systems was kept at 303K. Furthermore, ethanol, at the two different concentrations of 15% and 25% was chosen, and the two different systems are shown in Table 2 (I), and (II), respectively.

The radial density profiles () of PPO, PEO, ethanol, and water in these two systems were calculated in order to better describe the morphologies of the micelles self-assembled by P123 in the mixed water–ethanol dilute solution. First, we calculated the center of mass (COM) of each micelle. Then, the number density as a function of distance r to the COM was calculated. After that, the number density of each component was divided by the shell volume. Finally, the mean radial density profiles () were obtained by averaging the number density over all observed micelles, as shown in Figure 4. From Figure 4, we can see that the peak of block PPO (red line) corresponds to the smallest r, which means that block PPO chains form the core of the micelle. Moreover, the sequence of position r corresponding to the main density peaks of different components is block PPO, block PEO, ethanol, and water, respectively. This indicates that block PEO (green line) wraps block PPO as the shell of the micelle, and most ethanol is outside the micelle. Most water is mixed with the ethanol due to the similar peak location of water and ethanol. However, through comparing Figure 4b,d,f, we could find the small density distribution points formed by water (blue line) and ethanol (yellow line) were smaller than that of the block PEO, and larger than block PPO, which indicates that water and ethanol molecules are in the corona of the micelles. Furthermore, the small density peak of ethanol (yellow line) are shifted left in Figure 3f compared with Figure 3d, which means that increasing the concentration of ethanol causes ethanol molecules to move inside the PEO block.

Figure 4

The radial density profiles () of PPO, PEO, ethanol, and water correspond to the systems of 5% P123 in water in Table 1, (I) and (II) in Table 2. (a) the radial density profile of block PPO in system of 5% P123 in water; (b) the radial density profile of block PEO and water in system of 5% P123 in water; (c) the radial density profile of block PPO in system of 5% P123 and 15% ethanol; (d) the radial density profiles of block PEO, ethanol, and water in the system of 5% P123 and 15% ethanol; (e) the radial density profile of block PPO in the system of 5% P123 and 25% ethanol; (f) the radial density profiles of block PEO, ethanol, and water in the system of 5% P123 and 25% ethanol. Block PEO: green; block PPO: red; ethanol: yellow; and turpentine oil: purple.

3.3. Effect of Turpentine Oil on the Micelles’ Structure of P123 in the Mixed Solvent

To evaluate the turpentine oil effect on the micelles’ structure of P123 in the presence of ethanol, we studied three different complex systems including P123 block copolymer, ethanol, oil, and water by hybrid particle–field molecular dynamics simulation method. The system composition are reported in Table 2 (IV), (V), and (VI). In order to understand the effect of the turpentine oil on the self-assembly of micelles, the size distribution of the micelles (5% P123 in pure water, 5% P123 in 15% ethanol and water, 5% P123 in 15% ethanol, 1.5% oil, and water mixed solvent) were calculated and compared, as shown in Figure 5. Micellar Parameters of mean aggregation number