Production run

After the equilibration run, we want to remove the fictitious pulling on the terminal carbon and let the micelle be on its own to observe its properties. We instead apply another harmonic potential for the entire center of mass (COM) of the micelle to restraint the whole micelle to float near the center of the box. This is to prevent the micelle from float to the “wall” of the box, making the snapshots of the simulation confusing and harder to follow. We want the COM to be at about the box center with 1Å fluctuation, hence k_spring is set at 250 kJ/mol.nm². The pull-code is:

pull = yes
pull-ngroups = 1
pull-ncoords = 1
pull-group1-name = SDS
pull-coord1-type = umbrella
pull-coord1-k = 250
pull-coord1-geometry = distance
pull-coord1-groups = 0 1
pull-coord1-dim = Y Y Y
pull-coord1-origin = 3.25 3.25 3.25

Production run lasted 10.0ns with 1.0fs time step and output 200 frames for analysis.

Results analysis

The first fundamental question we can ask is how close the micelle configuration is compared to a “theorist’s drawing” that was shown in the beginning. Does the SD monomer indeed forms a “perfect” spherical or how close can it get? Does the monomer chain stay straight and align and point themselves toward the micelle center? MD analysis techniques used to address these questions are shown below.

Shape

To quantify the sphericality, we can first compute the micelle’s eccentricity \(e\) [1]:

\begin{equation} e=1- {I_{min}}/{I_{avg}} \end{equation}

Where \(I_{min}\) is the smallest moment of inertia (m.o.i) along either the x, y or z-axis and \(I_{avg}\) is the average of all 3 m.o.i. An eccentricity \(e\) approaching 0 means that the micelle is perfectly spherical in shape and a small deviation from 0 indicates some degree of ellipticity.

Next, the mean micellar radius \(R_s\) can be inferred from [1]:

\begin{equation} R_s = \sqrt{5/3} R_g \end{equation}

Where \(R_g\) is the gyration radius of the DS micelle. Both the three m.o.i and the gyration radius can be obtained from MD frames using gmx gyrate. Figure 5. below shows the time series of the three m.o.i and the gyration radius from the 10ns production run.

We found that \(e=0.203\) with \(R_g=17.1Å\) and \(R_s=22.1Å\). These values indicate the micelle is not quite spherical but has some ellipticity with an average radius of about 22Å.

Another way to infer the shape of our micelle is to look at the radial distribution function \(g(r)\) for certain atom types with respect to the center of mass (COM) of the micelle. The atom types we can look at are carbons (C1 to C12), sulfur in head groups, counter-ion Na^+, and oxygen of water. Intuitively, function gmx rdf -f md.trr -s md.tpr -n index.ndx -o rdf.xvg -selrpos whole_res_com -seltype atom should be able to produce such \(g(r)\), but our attempts to do so have failed so far, the reader can possibly explore more by themselves on how to use this function properly.

Instead, we used gmx traj to get the coordinates of all relevant atoms as well as of the COM of the micelle and import to python to compute \(g(r)\) manually using this formula:

\begin{equation} g(r) = \frac{\rho_{local}}{\rho_{avg}}=\frac{N_{local}/4\pi\Delta r^2}{N_{total}/4\pi r_{max}^2} \end{equation}

Where \(N_{local}\) is the count of the atom in question in a circle shell area with radius difference \(\Delta r\), which we chose to be 0.15Å. \(N_{total}\) is the total atom count in a circle with radius \(r_{max}\), which we chose to be 3.5nm since that is observed to be where the distribution trail off to 0 (fig.6).

It can be seen in Figure 6 (left) that the sulfur peak is very symmetric at about 20Å, this peak can also be used to infer the micelle radius. The Na⁺ peak at about 23Å indicates the counter-ion binding near the sulfur headgroup as expected. The water peak indicates that water does not really penetrate into the hydrophobic core of the micelle (< 15Å) and water mostly gathers outside the sulfur as a hydration shell, which stabilizes the micelle.

Figure 6 (right) can be used to visualize how the carbons in DS distribute in the micelle. It is interesting to note that [1] reports that carbons near the center (C12 is the nearest) have more mobility than the ones near the headgroup, hence C12 can be observed as having more noise and a broader peak in their g(r) profile.

Monomer chain direction/chain stiffness

Another analysis we can make is how well the 72 DS monomers are packing in the micelle. As seen in Figures 3 and 4, some DS chains are decently bent and do not point quite toward the center of micelle. This can be confirmed by looking at the C1-COM-C12 angle distribution and the C1-C12 distance distribution. Both of these can be found from the coordinates data we already obtained above.

As seen in figure 7 left, the angle C1-COM-C12 does not quite center around 0^o but rather around 40^o, which means that the DS monomers are not pointing perfectly straight at the COM.

Furthermore, the word “pointing straight” implicitly assumes that the C12-C1 bond through the other carbons is straight itself, which is not always the case by looking at the right figure. 13.9Å corresponds to the C12-C1 distance of a straightly built DS monomer, which is a C12-C1 trans dihedral. The distribution shows that the C12-C1 distance spends a lot more time in the <13.5Å region which corresponds to the gauche C12-C1 dihedral, or a bent carbon chain.

[Page 4: Appendix and reference]