Free energy profiles for penetration of methane and water molecules into spherical sodium dodecyl sulfate micelles obtained using the thermodynamic integration method combined with molecular dynamics calculations

When a solution dissolves micelles consisting of surfactant molecules, a water-insoluble substance can bind to the micelles and be dispersed in water. This phenomenon is called “solubilization.” One familiar example of solubilization is cleaning, in which a greasy stain is solubilized in micelles. In cosmetics,^1–5 a water-insoluble perfume is dispersed using micelle solubility. In pharmaceuticals production,⁶ a drug and a surfactant are often mixed to smoothly disperse the drug within the human body. Recently, drug-delivery systems using micelles to transport drugs have been actively developed.^7–9 In synthetic chemistry, synthesis is sometimes performed within reactant-solubilizing micelles.^13,14 Solubilization is often used in biochemistry to extract water-insoluble membrane proteins from cells.^10–12 In spite of the utility of micelles for solubilization in a wide range of basic scientific fields and applications, the molecular picture of solubilization is still not clear.

In previous physicochemical studies, the free energy difference between the micelle and bulk water for a water-insoluble molecule was obtained experimentally by measuring the bound fraction of the molecule in the micelle.^17–19 The binding site of the solubilized molecule has also been obtained using ultraviolet absorption spectroscopy^15,20–22 and x-ray diffraction.^23,24 Hydrophobic molecules such as alkanes are found in the center of micelles, whereas amphiphilic molecules such as alcohols are solubilized in the palisade layer of the micelle. The binding structure of solubilized molecules of this kind has been obtained by assuming a model structure for solubilization and analyzing the experimental data. In the present study, molecular dynamics (MD) calculations have been performed to obtain the binding structure of solubilized molecules without using such a model structure.

The phase-separation and mass-action thermodynamic models of solubilization phenomena are well known from previous theoretical studies. In the former, the micelle is regarded as an isolated phase (micellar phase) in bulk water, and solubilization is treated as the distribution of a water-insoluble molecule in the micellar phase.¹⁵ On the other hand, in the mass-action model, solubilization is regarded as a chemical reaction between a water-insoluble substance and the micelle.¹⁵ These models are powerful tools for understanding solubilization phenomenologically. However, few theoretical studies of solubilization phenomena starting from intermolecular interactions and statistical mechanics have been performed. To our knowledge, other than our own previous studies, a theoretical study of solubilization has been performed only by Matubayasi et al.,¹⁶ based on distribution function theory combined with MD calculations.

In our previous work,^23–27 we performed a series of MD calculations and found that the free energy of transfer of a water molecule from bulk water to a sodium dodecyl sulfate (SDS) micelle core has a large positive value, +28 kJ mol⁻¹, and the water molecule seldom solubilizes into the SDS micelle core.²⁵ The free energy of transfer was also calculated²⁶ for a series of alkanes, i.e., methane, ethane, butane, hexane, and octane molecules, from the water phase to the SDS micelle core. The calculations showed that the free energy decreases linearly as the number of carbon atoms in the alkane molecule increases, with a slope of 3.3 kJ mol⁻¹ per carbon atom. The enthalpy, ΔH, and entropy, ΔS, of transfer decrease linearly with increasing numbers of carbon atoms. The obtained free energy of transfer is, however, an average value inside the micelle. No information about the binding site was obtained for the penetrating molecule into the SDS micelle. It is essential to obtain information about both the binding site and the free energy change of the solubilized molecule in the SDS micelle. Furthermore, from the physicochemical and applications point of view, it is very important to know which interaction between the solubilized molecule and the SDS micelle produces the free energy change.

In this work, in order to clarify where the solubilized molecule is bound and how to stabilize it in the SDS micelle, the free energy profiles of penetration of hydrophobic methane and hydrophilic water have been calculated as a function of the distance between the center of mass of the SDS micelle and the penetrating molecule, using the thermodynamic integration method combined with MD calculations. The binding site and free energy barrier have been quantitatively evaluated. Furthermore, the solubilization free energy of the solubilized molecule was divided into four contributions, i.e., from the hydrophobic group, sulfate ions, sodium ions, and solvent water, in order to clarify which part of the solution is important. Finally, cavity distribution in the micelle was analyzed.

Computational details are given in Sec. II. The calculated free energy profile and distributions of the solubilized molecules are discussed in Sec. III. We conclude the paper in Sec. IV.

The method of calculating the free energy profile for penetration of methane and water molecules in SDS micelles and the computational details of the MD calculations are described in this section.

We consider the transfer process of the methane or water molecule from a position in the water phase far from a micelle to the micelle core. The system is composed of an SDS micelle in water and one methane or water molecule to be examined. Here, the methane and water molecules are called “solubilized” or “penetrating” molecules. Let r be the distance between the center of mass of the penetrating molecule and that of the SDS micelle, and let G(r) be the Gibbs free energy of the whole system, where the distance is constrained to r. The Gibbs free energy of transfer ΔG(r) of the penetrating molecule from a reference distance r₀ to r may be calculated as^28–30 

where F(r^′) is the force acting between the center of mass of the SDS micelle and that of the penetrating molecule, separated by a distance r^′, and 〈…〉 represents the isothermal–isobaric ensemble average. F(r^′) can be calculated as

where F_i and F_j are the forces acting on atom i and atom j, respectively, from all surrounding atoms, and m_s and m_m are the total masses of the penetrating molecule and the SDS micelle, respectively; i ∈ M_s and j ∈ M_m represent the ith atom in the penetrating molecule and the jth atom in the SDS micelle, respectively; and u is the unit vector from the center of mass of the penetrating molecule to that of the SDS micelle. The free energy profile can be evaluated from Eq. (2.1) using 〈F(r^′)〉 calculated from the MD.

In the MD calculations, an aggregation number of 60 was adopted for the present SDS micelle since our previous work,³¹ as well as light-scattering experiments^32,33 and time-resolved fluorescence spectroscopy,³⁴ indicated that the aggregation number is thermodynamically most stable at about 60. Thus, 60 SDS molecules, 8360 solvent water molecules, and one penetrating molecule were contained in a cubic simulation box with the periodic boundary condition. The SDS concentration in the present calculation is 50 CMC. TIP4P (Ref. 35) and CHARMM (Refs. 36 and 37) potentials were used for water and the other molecules, respectively. The pressure and temperature were controlled at P = 0.1 MPa and T = 300 K, using the algorithm proposed by Martyna et al.^38–40 The inertial constant of the barostat was set to 6.7 × 10⁻¹⁷ J s², so the relevant time constant is 2.0 ps. A fivefold Hoover chain thermostat was separately connected to particle degrees of freedom and the barostat. The inertial constants of the thermostat for the former were set to 3.6 × 10⁻²¹ J s² for the first chain and 1.2 × 10⁻²¹ J s² for the second to fifth chains, and those for the latter to 4.2 × 10⁻¹⁸ J s² for the first chain and 7.5 × 10⁻²³ J s² for the second to fifth chains; the time constants of the five chains are all 0.5 ps. The time step, Δt, was 2 fs. The particle mesh Ewald (PME) method⁴¹ was adopted to calculate the Coulombic interaction. The cutoff distance for the Lennard-Jones interaction and the real-space part of the PME method was 1 nm. The parameter α and the grid points of the PME method were 0.4 × 10¹⁰ m⁻¹ and 64 × 64 × 64, respectively. SHAKE/ROLL and RATTLE/ROLL methods were used to constrain the bond length between carbon (C) and hydrogen (H) atoms in the methyl and methylene groups of dodecyl sulfate ions, the bond length between oxygen (O) and hydrogen atoms in the water molecule, and the bending angles H–C–H in the methylene groups and H–O–H in the water molecule.

The initial configuration of the spherical SDS micelle was set such that the hydrophobic group of the dodecyl sulfate ion is inside the micelle and the hydrophilic sulfate ion is outside. Water molecules (8360) were set around the SDS micelle in the simulation box. MD calculations were performed for 2 ns to equilibrate the micelle structure and spatial distribution of the sodium ions.

The initial configurations for the mean-force calculation were generated using the above equilibrated configuration of SDS micelles in water. A single penetrating molecule was located at a distance r from the center of mass of the SDS micelle. Here, the mean force 〈F(r)〉 was calculated at seven distances: r = 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 nm. SHAKE and RATTLE methods were used to constrain the distance from the center of mass of the SDS micelle to that of the penetrating molecule. Each MD calculation was performed for 5 ns. The first 1-ns trajectory was excluded from the analysis. At r = 0.1 and 0.5 nm, MD calculations performed for 4 ns for the production run were sufficiently long to obtain the averaged value of 〈F(r)〉 because the trajectory of the penetrating molecule fully covers the whole spherical region of radius r, and satisfactory statistics were obtained from this. At r = 1.0, 1.5, 2.0, 2.5, and 3.0 nm, however, two MD calculations, starting from different initial configurations, were performed at each distance to obtain enough statistics for analysis.

In order to discuss differences between the penetration of the water molecule into the SDS micelle and the penetration of the water molecule into a hydrophobic droplet, as discussed in Sec. III C, the free energy profile of water in a hydrophobic dodecane droplet has also been calculated. The aggregation number of the dodecane droplet was set to 60, as for the SDS micelle. One dodecane droplet, 8397 water molecules, and one solubilized water molecule were contained in a cubic simulation box with the periodic boundary condition. TIP4P and CHARMM potentials were used for the water and dodecane molecules, respectively. The mean force 〈F(r)〉 was calculated at seven distances: r = 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 nm. Each MD calculation was performed for 5 ns. The droplet, which is thermodynamically in the metastable state, was stable during these calculations. The first 1-ns trajectory was excluded from the statistics, and the remaining 4-ns trajectory was used for analysis. Details of the MD calculation were the same as those for the SDS micelles stated above.

Figure 1(a) shows the calculated radial number density profile, ρ(r), of carbon atoms in the hydrophobic group, sulfur and oxygen atoms in the sulfate ion of the SDS micelles, solvent water molecules, and sodium ions. The ρ(r) obtained in the present study is in very good agreement with the results of other studies.^16,42,43 For instance, the peak position of the sulfate ion of the SDS micelle is almost the same as that calculated by Matubayasi et al.¹⁶ Furthermore, as in the previous studies, a low-density region is observed at the center of the SDS micelle in the present work. These results indicate that our calculations have been performed correctly.

In Fig. 2, cumulative averages of F(r) of solubilized methane (Fig. 2(a)) and water (Fig. 2(b)) molecules at r = 1.5 nm, evaluated using Eq. (2.2), are shown as a function of simulation time. The solid and dashed lines in Figs. 2(a) and 2(b) show the results of the MD calculations starting from two different initial configurations. At r = 1.5 nm, convergence is slowest among the various r values because of large fluctuations in F(r) in the boundary region between the hydrophobic groups and sulfate ions of the SDS micelles. It was found that two cumulative averages converged to the same value F(r = 1.5 nm) = −2.0 × 10⁻¹² N for the methane and F(r = 1.5 nm) = 33 × 10⁻¹² N for the water after 4-ns MD calculations. The 4-ns MD calculations are therefore long enough to obtain satisfactory statistics for the mean force 〈F(r)〉. In order to evaluate the error of the mean force 〈F(r)〉, independent intervals were determined using the analysis of Fincham et al.⁴⁴ In this research, the time average of each 200-ps trajectory was regarded as an independent measurement of F(r). Forty samples were obtained from two different 4-ns MD calculations at each distance of r = 1.0, 1.5, 2.0, 2.5, and 3.0 nm, and the statistical errors of 〈F(r)〉 were evaluated. At distances r = 0.5 and 0.1 nm, 20 samples were obtained from 4-ns MD calculations. The 80% confidence interval of 〈F(r)〉 was evaluated, and is shown by the error bar in Fig. 1.

Figure 1(b) shows the mean force 〈F(r)〉 of the solubilized methane and water molecules at distances r = 0.1, 0.5, 1.0, 1.5, 2.0, 2.5, and 3.0 nm. Positive values of 〈F(r)〉 indicate that the mean force acting on the solubilized molecule is a repulsive force, and excludes it from the SDS micelle. Negative values of 〈F(r)〉 indicate that the mean force acting on the solubilized molecule is an attractive force, and draws it into the SDS micelle. The figure shows that the interaction between methane and the SDS micelle is attractive at a range of 1.5 < r < 2.5 nm, and that the interaction between water and the SDS micelle is repulsive at 0.5 < r < 2.0 nm. At 0.1 ≤ r ≤ 0.5 nm, attractive interactions can be found for both methane and water. Details of the behavior are discussed in Sec. III C.

Figure 1(c) shows the free energy profile, ΔG(r), for penetration of methane and water. ΔG(r) can be obtained by integrating 〈F(r)〉 according to Eq. (2.1). In our previous work, the free energy of transfer, ΔG, of methane from bulk water to the SDS micelle core rather than the free energy profile ΔG(r) was calculated using the thermodynamic integration method with respect to the potential parameters. The free energy profile ΔG(r) for methane solubilization was also evaluated approximately by Matubayasi et al.¹⁶ using distribution function theory combined with MD calculations. The results are also plotted in Fig. 1(c). ΔG(r) in the region of the SDS micelle (0.1 ≤ r ≤ 1.5 nm) is consistent with the value of ΔG(r) = −5 kJ mol⁻¹ in our previous work²⁶ and is in reasonable agreement with the ΔG(r) calculated by Matubayasi et al.

A free energy barrier is not observed between the micelle core and the bulk water, so methane can be easily solubilized into SDS micelles in the process of penetrating the micelle. The free energy in the SDS micelle core (r = 0.1 nm) is about 6 kJ mol⁻¹ lower than in the hydrophilic part (r = 1.2–2.5 nm). This indicates that from the viewpoint of excess free energy, the methane is more stable in the SDS micelle core than it is in the vicinity of the hydrophilic sulfate ions. However, the free energy difference (6 kJ mol⁻¹) between the micelle core and the hydrophilic part is not so large, only 2.4 times as large as the thermal energy at temperature T = 300 K. The solubilized methane might therefore be nonlocalized in the micelle core and move around in the whole of the micelle. Furthermore, the free energy difference between the hydrophilic part and the bulk water is also 6 kJ mol⁻¹, so the methane can easily go in and out of the SDS micelle. Distribution of methane in the micelle will be discussed in more detail in Sec. III D.

With respect to the penetrating water molecule, a large positive ΔG(r) value, around 24–35 kJ mol⁻¹, is observed in the region r ≤ 1.0 nm in Fig. 1(c). This value is in good agreement with the free energy of transfer of water from the water phase to the SDS micelle core of 28 ± 4 kJ mol⁻¹ in our previous work.²⁵ It is certain that no water molecule permeates into the micelle core region in the present calculation (Fig. 1(a)), where ΔG(r) is larger than 10 kJ mol⁻¹.

The minimum ΔG(r) of methane and water is found at the center of the SDS micelle in the region r ≤ 0.5 nm. This feature is also observed in the free energy profile of a lipid membrane at the center of the membrane.^28,30 This phenomenon is commonly observed in the low-density region of the hydrophobic groups of surfactant and lipid molecules, i.e., the contact region of the tail end.

We clarify how hydrophobic groups and hydrophilic sulfate ions of SDS micelles, sodium ions, and solvent water affect the penetration of methane and water molecules into the SDS micelles. First, the force F(r) acting on the penetrating methane and water was divided into contributions from the hydrophobic groups and sulfate ions of the SDS, sodium ions, and solvent water. When F_i and F_j in Eq. (2.2) are divided into these four contributions, F(r) can be written as

where the superscripts on F_i, F_j, i.e., pho, phi, sod, and w, stand for the hydrophobic groups and sulfate ions of the SDS, sodium ions, and solvent water, respectively. For example, F_i^pho found in the first term, represents the force acting on atom i of the penetrating methane or water from the hydrophobic groups of the SDS, and F_j^pho represents the force on atom j of the hydrophobic groups of the SDS from the penetrating methane or water. Thus, the first, second, third, and fourth terms in the right-hand side of Eq. (3.1) are the contributions to the mean force 〈F(r)〉 from the hydrophobic groups and sulfate ions of the SDS, sodium ions, and solvent water, respectively.

Formally integrating each contribution over r, we obtain a free energy-related quantity, whose summation gives the correct free energy. Of course, the quantity is not the thermodynamically defined free energy, but it works as a measure of the contribution to the free energy of interest. Here, we describe this quantity using the same symbol, ΔG(r), as that used for the free energy change. It should be noted that an infinite number of paths may exist in an actual penetrating process, taking roundabout routes to the center of the micelle. However, in our analysis, the path dependence is averaged out and, at the same time, the force in only the r direction is investigated. Although contributions to the force acting on a penetrating molecule along an actual path are path dependent, we believe that the decomposition analysis in the present study is very useful in addressing the question of what role each component plays in the permeation.

In Fig. 3(a), the calculated mean force 〈F(r)〉 acting on the penetrating methane is shown with four contributions, i.e., from the hydrophobic groups and sulfate ions of the SDS micelle, sodium ions, and solvent water. The red and brown lines in Fig. 3(a) indicate that contributions from the sulfate ions and sodium ions are small in the whole region of r. This is because methane is a nonpolar molecule and interacts weakly with the sulfate ions and sodium ions. Even in the region 1.2 ≤ r ≤ 2.5 nm, where the sulfate ions and sodium ions are largely distributed, contributions from the sulfate ions and sodium ions are quite small. The weak Lennard-Jones interactions between methane and the sulfate ions and sodium ions indicate that the methane is not in close contact with them. In contrast, contributions from the hydrophobic part and solvent water are dominant at r = 1.5 and 2.0 nm, the interface between the micelles and the water phase. Thus, contributions from the hydrophobic groups and solvent water are important for the solubilization of methane in the SDS micelle.

We now consider the solubilization of methane in the SDS micelle using the schematic diagram shown in Fig. 4. The number densities of the hydrophobic groups of the SDS micelle are represented by light and dark shades of gray; the darker the gray color, the higher the number density. Closed circles represent the methane molecules. The green and blue arrows indicate contributions to the mean forces 〈F(r)〉 on the methane from the interaction with the hydrophobic groups and solvent water, respectively. The inset in Fig. 4 shows the radial number density profile, ρ(r), of hydrophobic carbon atoms of the SDS micelle and solvent water.

At r = 3.0 and 2.5 nm in the water phase, only a very small force from the hydrophobic groups and solvent water was observed, as shown in Fig. 3(a). The small interaction between methane and the hydrophobic groups is simply a result of the distance between them, i.e., the Lennard-Jones force is short range. The small averaged force from the solvent water, in spite of the close contact, must be caused by symmetric hydration around the methane, which is located at a distance from the interface with the micelle. So, no arrow is drawn in Fig. 4. At r = 2.0 nm, however, an attractive force (blue arrow) acts on methane from the solvent water. A repulsive force from the hydrophobic groups is also acting on the methane (green arrow). However, as shown in Fig. 4, the attractive force is stronger than the repulsive one, so the methane is drawn into the SDS micelle. The attractive behavior of the total mean force in this interface region is, of course, caused by the hydrophobic interaction, which avoids the entropy loss of hydrophobic hydration of the methane. Here, we must remember that the hydrophobic interaction is a solvent-induced interaction, where water–water interactions play an essential role. Thus, in the present case, we should understand that, first, the solvent water molecules tend to gather together in order to reduce contact with the methane molecule located at the interface with the micelles. Second, this tendency is observed from the force on the methane from the water molecules (blue arrow), excluding it from the water phase. Third, however, a small amount of work is required to push the methane molecule into the high-density region of the hydrophobic groups (green arrow). At r = 1.5 nm, the situation is similar to that at r = 2.0 nm, but the attractive and repulsive forces are balanced. At r = 1.0 nm, no contribution from solvent water is observed because the methane is not in contact with the solvent water at this distance, as shown in Fig. 1(a). With regard to the contribution from hydrophobic groups, almost no force is acting on the methane at r = 1.0 nm because the density gradient of the hydrophobic groups is small at that point, so the solvation structure may be considered to be symmetric. In the region where r < 1.0 nm, an attractive force associated with the hydrophobic part acts on the methane because of the lower density inside the SDS micelle.

Contributions to the free energy profiles were evaluated by integrating 〈F(r)〉 according to Eq. (2.1). The resulting profiles for four contributions are shown in Fig. 3(b). In this figure, contributions to the free energy profiles from the sulfate ions and sodium ions were found to be small, consistent with the results for the mean force 〈F(r)〉. The methane in the SDS micelle is stabilized by about 25 kJ mol⁻¹ because of hydrophobic interaction induced by the solvent water, but it is destabilized by about 5–15 kJ mol⁻¹ by the hydrophobic groups. Solubilization of methane into the SDS micelle therefore comes from avoidance of hydrophobic hydration of the methane.

Figure 5(a) shows the contributions to the mean force acting on the penetrating water from hydrophobic groups, sulfate ions, sodium ions, and solvent water. These four contributions change in complex ways compared with the case of methane. The contributions from sulfate ions, sodium ions, and solvent water cancel one another. Thus, the total value becomes small, as shown in Fig. 5(a). This is also the case for the free energy profile shown in Fig. 5(b), where the contributions from sulfate ions, sodium ions, and solvent water to the free energy profile also cancel one another, and as a result, the total free energy change is small compared with the component changes.

In order to simplify the discussion, the mean force of water transfer to the dodecane droplet has also been calculated as a reference, using the thermodynamic integration method combined with MD calculations. The radial number density, ρ(r), of the hydrophobic groups of the dodecane and water molecules is presented in Fig. 6(a). Figure 6(b) shows the mean force between the water molecule and the dodecane droplet with two contributions to the mean force of penetrating water from the hydrophobic groups and solvent water. The calculated free energy profile ΔG(r), together with its two contributions relevant to the decomposition of the force in Fig. 6(b), is given in Fig. 6(c). First, we discuss this reference system before we investigate the micelle system of interest. It was found that at r = 1.0 and 1.5 nm, i.e., the interface region, a repulsive force acts on the water both from the hydrophobic groups and solvent water, resulting in the high free energy barrier found in Fig. 6(c). This result is consistent with the well-known view that a water molecule is less stable in a hydrophobic environment than in the water phase because of the energetic disadvantage of breaking the hydrogen bonds formed in the water phase. This view may also be applied to the case of SDS, where a water molecule is unstable in the hydrophobic core of the micelle compared with that in the water phase. The degrees of instability are similar to each other: 24–35 kJ mol⁻¹ for SDS and 25 kJ mol⁻¹ for the dodecane droplet. However, the behavior found at the interface region is very different. An essential difference between SDS micelles and dodecane droplets is the sulfate and sodium ions found in the interface region of the SDS micelles (see Fig. 1(a)). The radial distribution of ions shown in the figure forms a strong dipole moment at the interface, where sodium ions are widely distributed outside the location of sulfate ions. The structure of the solvent water must therefore be strongly influenced by this surface dipole moment or the two ions.

It is reasonable to consider that the force on a penetrating water molecule around the interface is dominated mostly by this dipole moment. The calculated mean forces shown in Fig. 5(a) may then be well understood as follows. At r = 3.0 nm, far from the interface, the mean forces acting on the penetrating water are all almost zero since the water molecule is in a symmetric field, as in bulk water. At r = 2.5 nm, near but outside the surface, forces from the solvent water and sodium ions cancel each other out. At r = 2.0 nm, at the interface, the forces from sulfate ions and sodium ions are great compared with the other components. Although the sign of their forces is different, as a result of the complicated structure of the interface, it is important to find that the two forces, about +70 × 10⁻¹² N and −50 × 10⁻¹² N from sulfate ions and sodium ions, respectively, or force from the dipole moments tend to keep the penetrating water molecules inside the water phase. The small negative force from the solvent water is no more than a secondary one after the ions at the surface have determined the structure, in particular, the orientation of the penetrating water molecule as well as the orientation of the solvent water molecules located across the ions. The orientations must be opposite, resulting in the negative force. At r = 1.5 nm, the dipole moment at the surface strongly draws the penetrating water back to the water phase, where the sign of the force is the same for sulfate ions and sodium ions. Furthermore, high-density hydrophobic tail groups also tend to exclude the water molecule from the micelle core. Together with the counterforce from the solvent water, the total mean force on the water molecule is about 30 × 10⁻¹² N, drawing it back to the water phase. At r = 1.0 nm, the situation is similar to the case at r = 1.5 nm, except for the values. At r = 0.5 and 0.1 nm, in the center of the hydrophobic core of the micelle, the averaged forces are small because of the long distance from the interacting groups or the symmetric force at a distance from the interface.

The radial distribution f(r) of methane from the center of mass of the SDS micelle may be evaluated by

where ΔG(r), R, and T are the free energy profile of methane as a function of r, the gas constant (R = 8.314 J mol⁻¹), and the absolute temperature, respectively. In Fig. 7, f(r) is plotted together with f_t(r), which is the radial distribution of methane evaluated from the 15-ns trajectory of MD calculations in our previous work. At r = 1.5 nm, f_t(r) is reduced to the same value as f(r). As shown in the figure, the shape of f(r) is in good agreement with f_t(r). The radial distribution, f_v(r), of the space that can accommodate a sphere of diameter of 0.3 nm, i.e., a cavity, in the micelle is also plotted in Fig. 7. f_v(r) is also reduced to the same value as f(r) at r = 1.5 nm. A good correlation between f_v(r) and f(r) indicates that both methane and cavities are distributed in the micelle, and the methane is considered to move about among the micelle cavities. As stated in Sec. III B, solubilized methane can come in and out through the surface of the micelle, so a distribution of methane can be seen outside the micelle at r ≥ 2.5 nm in f(r) and f_t(r).

Both f(r) and f_t(r) show a high probability in the vicinity of the center of the micelle. The number of solubilized methane molecules can be calculated by multiplying the Jacobian 4πr² and f(r) together, as shown in the inset of Fig. 7. Although the error included in the value is large, we may qualitatively discuss the actual distribution of methane in the micelle. It is interesting to find that after the volume factor 4πr² is taken into account, the methane is most often observed at the interface.

In the present work, ΔG(r) for penetration of methane and water molecules into SDS micelles have been calculated using the thermodynamic integration method combined with MD calculations.

Methane is about 6 kJ mol⁻¹ more stable in the vicinity of the sulfate ions of the micelle than it is in the water phase. No free energy barrier is observed there, so the methane seems to be easily drawn into the micelle. Furthermore, methane in the center of the micelle is about 6 kJ mol⁻¹ more stable than it is in the vicinity of the sulfate ions. The binding energy is about 2.4 times larger than the thermal fluctuation energy at T = 300 K. This value is relatively small, and therefore the solubilized methane may move about in the micelle core.

The relative contributions of the hydrophobic groups and sulfate ions of the SDS, sodium ions, and solvent water to methane solubilization in SDS micelles have been clarified. The contributions of solvent water and hydrophobic groups are dominant. With respect to the contribution from solvent water, in order to avoid hydrophobic hydration of methane in solvent water, the methane tends to be excluded from the water phase to the micelle. Considering the mean force just from the solvent water, the methane in the micelle is about 25 kJ mol⁻¹ more stable than it is in the solvent water. In contrast, in the crowded hydrophobic groups of the micelle core, the methane is destabilized by 5–15 kJ mol⁻¹ because of the repulsive interactions between methane and the hydrophobic groups. As a result, the methane remains in the micelle because the stability provided by the hydrophobic interaction is greater than the instability caused by the repulsive interactions between the methane and the hydrophobic parts.

The ΔG(r) of the penetrating water is evaluated to be 24–35 kJ mol⁻¹, so penetration of water into the micelle hardly occurs. The mean force 〈F(r)〉 on the water was divided into four contributions. These contributions change in complex ways as a function of r, and cancel one another out. In order to solve this puzzle, we calculated the mean-force profile for transfer of water into the dodecane droplet as a reference, using MD calculations. The essential mechanism of the penetration of a water molecule is the same for the SDS micelle and the dodecane droplet; the water molecule is unstable in the hydrophobic core of the micelle and in the dodecane droplet compared with in the water phase. The complex behavior observed in the SDS micelle is from the strong dipole moment formed at the interface by sulfate ions and sodium ions.

Finally, the radial distribution of solubilized methane has been evaluated from exp{−ΔG(r)/RT}. A good correlation was found between the radial distribution of 0.3-nm cavities in the SDS micelle and that of the methane. The methane is considered not to be localized but to move about in the SDS micelle using cavities in the hydrophobic core.

This work was supported by the Next Generation Super Computing Project, Nanoscience Program, and by TCCI/CMSI in the Strategic Programs for Innovative Research, MEXT, Japan. The work was also supported by a Grant-in-Aid for Scientific Research (No. 18066020) from the Japan Society for the Promotion of Science. The authors thank the Okazaki Research Center for Computational Science, National Institute of Natural Science, the Center for Computational Science, University of Tsukuba, and the Information Technology Center of Nagoya University for the use of supercomputers.