New Coarse-Grained Models to Describe the Self-Assembly of Aqueous Aerosol-OT

TL;DR

This paper introduces new simplified models to study how AOT surfactant forms structures in water, enabling faster and more efficient simulations.

Contribution

The development of three coarse-grained models for AOT using the Martini 3 force field to simulate self-assembly behaviors.

Findings

The models accurately reproduce AOT self-assembly at different concentrations.

Vesicle formation from bicelles above the critical vesicle concentration is demonstrated.

The models enable studying large-scale structures and long time scales not feasible with all-atom simulations.

Abstract

Aerosol-OT (AOT)
is a very versatile surfactant that
exhibits a
plethora of self-assembly behaviors. In particular, due to its double-tail
structure, it is capable of forming vesicles in water. However, the
size of these structures, and the time scales over which they form,
make them difficult to study using traditional all-atomistic molecular
dynamics simulations. Here, three coarse-grained models are developed
for AOT with different levels of detail. The models take advantage
of the Martini 3 force field, which enables 2:1 mappings to be employed
for the tail groups. It is shown that these models are able to reproduce
the self-assembly behavior of AOT in water at three concentrations:
below the critical vesicle concentration (CVC), above the CVC, and
in the lamellar phase. The results also demonstrate the formation
of vesicles from bicelles above the critical vesicle concentration,
which is an important milestone for the continued study of vesicle
behavior.

Figures

Skeletal structure of Aerosol-OT.

Coarse-grained mappings of AOT at different resolutions: (a) Finest, (b) Mixed and (c) Coarsest. Sets of highlighted atoms and bonds are mapped to a single bead. Symmetrically equivalent beads are highlighted in the same colors; for a given mapping, bonds between pairs of beads of the same colors have the same parameters. For example, the bonded parameters between the red and purple tail groups of the Coarsest model are identical.

Illustration of the coordinate-pair eccentricity (CPE) values for different types of ellipsoids. Adapted from Gale et al. Copyright 2022 American Chemical Society.

Normalized aggregation numbers (fraction of total number of surfactants in a cluster) for each aggregate (aggregation number ≥5) in the dilute isotropic systems. Only aggregates with an aggregation number of at least 5 are shown.

Final configurations of the Mixed dilute isotropic systems at (a) 0.27 wt % and (b) 1 wt % AOT. Yellow beads represent surfactant headgroups and gray beads represent tail groups. The box in 5­(a) is isotropic, of side ≈ 37 nm.

Mean surface area to volume ratio, surface area per surfactant, volume per surfactant, and aggregation number of aggregates (aggregation number ≥5) in 100 ns bins for each of the 1 wt % systems. Shaded area represents the 95% confidence interval calculated using bootstrapping.

Stages of bicelle-to-vesicle transition for the Mixed model. The initial bicelle (a/e) begins to curve in on itself (b,c/f,g), eventually closing to become a vesicle (d/h). The top row shows the side view and the bottom row shows the top view at the same time step.

Mean surface area to volume ratio, surface area per surfactant, volume per surfactant, and aggregation number of aggregates (aggregation number ≥5) in 100 ns bins for each of the 0.27 wt % systems. Shaded area represents the 95% confidence interval calculated using bootstrapping.

Coordinate-pair eccentricities of every aggregate (aggregation number ≥5) across the entire trajectory for each of the coarse-grained simulations. See also Figure for a schematic interpretation of this mapping.

Combined kernel density estimate and box-and-whisker plots showing the distributions of (a) thickness and (b) area per headgroup of the AOT bilayer for each CG model. The distributions are aggregated across t > 30 ns. Dashed lines indicate experimental values. Experimental values were derived by Fontell using X-ray diffraction.

Tables

1. Number of AOT Molecules and Water Beads for the Coarse-Grained AOT Simulations

Model	Concentration (wt % AOT)	Number of AOT	Number of Water Beads	Initial Box Size (nm)
Finest	0.27	200	452043	38
	1.0	200	128903	25.2
	18.0	300	8453	
Mixed	0.27	200	453128	38
	1.0	180	111937	24
	18.0	340	9667	
Coarsest	0.27	200	453348	38
	1.0	180	112134	24
	18.0	350	9965	

2. Time-Averaged Mean Properties of the Systems (Averaged across Aggregates with Aggregation Number ≥5) in the 0.27 Wt % Simulations Compared to Experimental Values

Model	Final Convergence Time (μs)	Mean Surface Area to Volume Ratio (Å–1)	Mean Surface Area per Surfactant (Å2)
Finest	12.6	0.288 ± 0.004	180 ± 4
Mixed	10.3	0.303 ± 0.002	188 ± 3
Coarsest	3.08	0.265 ± 0.015	238 ± 16
Experiment	–	0.15	194.3

3. Mean Bilayer Thickness and Area Per Headgroup for Each CG Model

Model	Thickness (Å)	Area per headgroup (Å2)
Finest	18.96	55.71
Mixed	19.32	60.25
Coarsest	22.20	70.00
Experiment	20.0	64.0

Full text

Introduction

Aerosol-OT (sodium bis­(2-ethylhexyl) sulfosuccinate) has attracted much interest because of its ability to self-assemble into robust reverse micelles in a plethora of solvents. It has been used to study confined water behavior? to template nanoparticle synthesis? and to stabilize microemulsions.? (We refer to Gale et al.? for a review on Aerosol-OT and its applications.)

The structure of Aerosol-OT (AOT) is shown in Figure. It has been suggested that the “cone-like” molecular structure of AOT is responsible for its phase behavior, but Nave et al.? argue that this is a myth; AOT simply has a convenient phase boundary around room temperature, and other sulfosuccinates are just as good at stabilizing microemulsions, but at different temperatures or using different oils. ?,?−? ?

Skeletal structure of Aerosol-OT.

One thing that each of the sulfosuccinate surfactants have in common is their doubled tail groups. According to the packing parameter theory of Israelachvili,? such surfactants lend themselves to forming vesicles. AOT is no exception; indeed, it has a characteristic critical vesicle concentration (CVC) of 7.5 mM in water, above which its micelles transition into vesicles.?

Vesicles have been studied extensively in a biological context, as they play an essential function in several physiological processes? and they show promise in enhancing drug delivery systems.? More recently, another important application of vesicles has been discovered: the control of wetting and bouncing behavior of droplets. Cryo-TEM images of aqueous systems containing AOT above and below the CVC demonstrate that vesicular AOT facilitates a wetting transition on hydrophobic leaf surfaces, thereby reducing droplet bounce.?

To elucidate the mechanisms responsible for these observations, and clarify whether vesicles are indeed essential for these types of phenomena, it is desirable to couple experiments such as those in the literature with molecular simulations. In fact, because of its wide range of applications, AOT has attracted the interest of many computational studies, which have been, for the most part, conducted at the atomistic level. ?−? ? ? ? ? ? ?

Vesicles are relatively large self-assembled structures, whose size is highly variable. For example, Fan et al.? reported an average radius of AOT vesicles in water at 10 mM of approximately 150 nm, as obtained using dynamic light scattering (DLS) and transmission electron microscopy (TEM) measurements, while Velázquez et al.? reported a radius of 40 nm under the same conditions using cryo-TEM. Irrespective of the exact radius, the size range observed experimentally imposes the use of large simulation boxes and a very large number of self-assembling structures, should one desire to observe the formation of vesicles during the course of unbiased simulations. Because of the related computational constraints, coarse-grained (CG) models are required, rather than atomistic ones.

Coarse-grained models allow one to reduce the number of particles in a molecular dynamics simulation by treating groups of atoms as a single unit (a bead). The bead represents the collective motion of the group of atoms, as well as the net interaction with other beads. The scaling of the number of particles with the simulation dimensions is significantly reduced by coarse-graining, allowing larger systems to be studied. Furthermore, because the degrees of freedom are reduced, CG simulations can study longer time scales compared to atomistic ones within the same computational time. However, care must be taken to ensure the models are accurate.

Here, we develop a CG model for AOT and show that it is capable of reproducing experimental self-assembly behavior. While the level of detail achieved with our method is less than atomistic models, ?−? ? ? the coarse-grained model presented here is capable of reproducing several experimental observables such as bilayer thickness and the transition of micelles to vesicles as the surfactant concentration increases.

Method

To enhance the atomistic models available in the literature, our goal is to develop a realistic coarse-grained model for AOT. For this purpose, we used the recently released Martini 3 force field to establish bead interaction parameters, as it boasts more accurate predictions of molecular packing than its predecessor,? as well as the option to include varying ratios of atom-to-bead mapping, which enables chemically similar molecules to be distinguished when a 4-to-1 mapping would be imprecise.? Furthermore, the Martini 3 force field provides more detailed models of ions through the use of beads of different sizes.

Others have attempted to derive coarse-grained models for AOT surfactants. For example, Negro et al.? developed one such model using Martini 2.? However, the authors admit that the model is “somewhat generic” because the two tail group beads could represent any alkyl chain between 6 to 10 carbons long. Because the tail groups have a pronounced effect on behavior ?,? it is desirable to develop a model that is coarse grained enough to provide economy of computational power, but detailed enough to allow us to study the effect of variations of the molecular architecture on its properties.

In our approach, CGBuilder? was used to map atoms to coarse-grained beads. Three mappings were generated using different scales of mapping for the carbon tails: Coarsest uses exclusively 4:1 mappings (R beads), Finest uses 2:1 mappings (T beads) and Mixed uses a combination. The resulting mappings are shown in Figure. Bead types were selected based on the suggestions by Souza et al.?

Coarse-grained mappings of AOT at different resolutions: (a) Finest, (b) Mixed and (c) Coarsest. Sets of highlighted atoms and bonds are mapped to a single bead. Symmetrically equivalent beads are highlighted in the same colors; for a given mapping, bonds between pairs of beads of the same colors have the same parameters. For example, the bonded parameters between the red and purple tail groups of the Coarsest model are identical.

Simulations of a single AOT molecule in water were performed for each of the three mapping schemes. The distribution of bond lengths and angles in the resulting trajectories were compared to a 50 ns all-atomistic (AA) simulation of the same system using the force field developed by Abel et al.? The AA trajectory was mapped to an effective coarse-grained (CG) trajectory by determining the center-of-geometry (mean coordinates) of the constituent atoms for each bead.

The bonded parameters must be fit to the bond length and angle distributions from the AA trajectory. This optimization was achieved using SwarmCG.? SwarmCG implements swarm optimization coupled with Boltzmann inversion to fit the bonded parameters to the AA target distribution. Boltzmann inversion has been successfully employed in parametrizing CG models, such as in the case of nanoparticles? and swarm optimization acts to efficiently fine-tune the parameters. Because of the symmetry of the AOT molecule, symmetrically equivalent bonds in the tails were assigned the same parameters and optimized jointly (see Figure).

In the Martini family of force fields, it is standard to exclude nonbonded interactions between directly bonded particles.? However, in AOT both carboxylate beads are close to the sulfonate group, but only one of them is bonded to it. The nonbonded, intramolecular interaction between these beads caused simulation instability and led to a poor match in the bond angle distribution between the three groups. The nonbonded, intramolecular interactions between the sulfonate and both carboxylate beads were therefore excluded. The resulting bond angle and bond length distributions are shown in the Supporting Information. The resulting topology files, in GROMACS ITP format,? are available in the Supporting Information.

Validation

Because we are ultimately interested in studying vesicles, which are self-assembled structures, it is pertinent to validate the models based on their self-assembly behavior. To that end, three sets of simulations were performed to study the properties of different types of structures: bilayer, micelles, and vesicles.

Dilute

Isotropic Phase

To study the properties of self-assembled micelles and vesicles, systems of 1 and 0.27 wt % AOT in water were generated for each CG model. These systems were created by randomly dispersing AOT molecules in a cubic box of water. The number of molecules in each system and the initial box dimensions are given in Table.

1: Number of AOT Molecules and Water Beads for the Coarse-Grained AOT Simulations

After setting up the initial configuration, the following procedure was used for the coarse-grained simulations:

• 1.An energy minimization step using steepest descent integration, for 100000 steps or until the maximum force in the system was no greater than 100 kJ mol^–1^ nm ^–1^.
• 2.For the 0.27 wt % systems only, an NVT simulation was conducted to increase the temperature from 0 K to 298 K across 5 ns (Δt = 10 fs).
• 3.An NPT equilibration for 10 fs (Δt = 10 fs) using a stochastic cell rescaling barostat? (τ = 4 ps).
• 4.Finally, a production run using an NPT ensemble (Δt = 20 fs) using a Parrinello–Rahman barostat? (τ = 12 ps). The 1 wt % systems were run for 2.5 μs, while the 0.27 wt % systems were run for 23 μs.

For each of steps 2 to 4, a V-rescale thermostat? (τ = 1 ps) was used to control the temperature. For the 0.27 wt % systems, the solvent has many more degrees of freedom than the solute and therefore it is more sensitive to stochastic velocity rescaling. This makes these systems susceptible to the “hot-solvent, cold-solute” effect when coupling the entire system using one thermostat, whereby the average temperature of the solute is much lower than that of the solvent. ?−? ? To overcome this artifact, it was necessary to thermostat the AOT and the solvent (water and Na^+^ ions) separately for the 0.27 wt % systems.

In order to characterize the aggregates formed during the simulations, a clustering algorithm was applied to the final production run. If the distance between two surfactant tail group beads was smaller than a threshold value, r cut, the parent molecules were considered to be part of the same aggregate. This threshold value was set to $[eqn]$ , where σ_LJ_ is the distance at which the Lennard-Jones potential between the smallest tail group beads is zero (r cut = 5.875 Å for the Coarsest model and r cut = 4.250 Å for the Mixed and Finest models). This choice of threshold respects the difference in the equilibrium distance between the tail groups of the different models, and was empirically observed to yield reasonable results.

Several properties of the aggregates were calculated. The radius of gyration, R g, was determined using

$$[eqn]$$

where N is the number of surfactants in the aggregate, r⃗_ i _ is the displacement of the ith surfactant from the center of mass of the aggregate and *m_i_

• is the mass of the ith surfactant.

The Willard–Chandler surface? of each aggregate was calculated using the Pytim? and PyVista? Python packages. The kernel density estimate of each aggregate’s surfactants’ beads was computed using a Gaussian kernel with a bandwidth of 4.0 Å. This can be thought of as a smoothing parameter. If it is too high, the geometry of the aggregate’s surface may be lost. If it is too low, voids may appear in the surface, corresponding to regions of lower density within the aggregate.

The code was amended to accelerate the computation of the density estimate using JAX.? Rather than using the minimum image convention, a supercell of 3 × 3 × 3 was created and the kernel was evaluated using all positions within this supercell. This allowed the use of JAX’s GPU acceleration when calculating the density estimate, as it does not support periodic boundary conditions and instead requires the input coordinates to exist in linear space.

The density estimate was then evaluated on a grid with spacing 2.0 Å. The Willard–Chandler isosurface was then generated using a density cutoff of $[eqn]$ , where ρ_max_ was the maximum computed density for the aggregate. The enclosed volume of the surface is negligibly larger (approximately 5% on average) using the JAX-accelerated method, but the computation time is improved significantly. The choice of cutoff and bandwidth were found empirically to yield surfaces that were very close to the locations of the headgroups, while avoiding artifacts (like gaps in the surface) that could be caused by fluctuations of the density within the core of the micelles. The volume and surface area of each Willard–Chandler surface were calculated to determine the volume and surface area per surfactant of each aggregate.

Following the procedure of Gale et al.,? the coordinate-pair eccentricities (CPEs) of the aggregates were also determined. An ellipsoid has three principal axes, a ≥ b ≥ c. The moments of inertia of the ellipsoid are related to the principal axes by

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$

The principal axes can therefore be solved from the moments of inertia algebraically.

In order to characterize this shape without losing information, two eccentricity values can be calculated:

$$[eqn]$$ $$[eqn]$$

Depending on the lengths of the three principal axes, and therefore the eccentricity values, four types of ellipsoids can be characterized:

Spherical ellipsoids have a ≈ b ≈ c and e _ ab _ ≈ e _ ac _ ≈ 0.

Prolate ellipsoids have a ≫ b ≈ c and e _ ab _ ≈ e _ ac _ ≫ 0.

Oblate ellipsoids have a ≈ b ≫ c and e _ ab _ ≪ e _ ac _ ≫ 0.

Triaxial ellipsoids have a b ≫ c and e _ ab _ < e ac ≫ 0.

These different types of ellipsoids can be distinguished graphically by plotting the CPE values; see Figure for a schematic example.

Illustration of the coordinate-pair eccentricity (CPE) values for different types of ellipsoids. Adapted from Gale et al. Copyright 2022 American Chemical Society.

Each of these properties was calculated in 2 ns intervals across the simulation trajectory for aggregates with aggregation number ≥5. The simulations were considered equilibrated once these values had reached a steady state. As discussed in the Dilute isotropic phasein Results section, for some systems it was not possible to achieve a steady state despite sampling across several microseconds.

In order to assess the convergence of the simulations, pyMSER ? was used to analyze the timeseries of the mean system properties, averaged across simultaneously coexisting aggregates. That is, the properties were calculated for each aggregate at each time step, then reduced to one mean system property per time step, for each of the following: volume per surfactant in aggregate, surface area per surfactant in aggregate, radius of gyration, and aggregation number. The evolution of these mean properties was compared across simulation time. pyMSER implements an algorithm to determine the optimal truncation point, k, for timeseries data in order to minimize the marginal standard error (MSE) of the remaining data. The MSE is defined as

$$[eqn]$$

where n is the total number of timesteps, Y _ i _ is the value of the property at time step i, and $[eqn]$ _ n,k _ is the mean of the property from time step k to n. Intuitively, the prefactor (n – k)^−2^ biases the metric toward maximizing the remaining sample size, therefore improving confidence, while the sum term penalizes the deviation of the remaining data from the mean. The data from k to n is then used to calculate the equilibrium mean and standard deviation.

In addition, pyMSER employs the augmented Dickey–Fuller test to verify that the remaining data is stationary.? Briefly, this method tests the null hypothesis that a unit root is present in the data, which would indicate that the data is nonstationary. The test statistic is compared to a critical value, which is dependent on the number of observations and the significance level. If the test statistic is less than the critical value, the null hypothesis is rejected and the data is considered stationary. For each of the dilute isotropic phase systems, convergence was tested by confirming that the equilibrium data for the aforementioned properties was stationary with at least 99% confidence.

We note that the geometrical properties measured here do not necessarily reflect internal molecular rearrangements within the aggregates, which have been shown to be significant particularly in the case of multicomponent systems with large amphiphiles, such as those found in biological systems. ?−? ? ? ? However, there is insufficient experimental data regarding the internal conformation of AOT aggregates to validate the behavior of the models in this respect. Furthermore, because AOT is a comparatively small surfactant, and the tail groups are chemically homogeneous, it is expected that the conformational space within an aggregate is limited.

Bilayers

The properties of spontaneously assembled bilayers were studied to compare to experimental data. A system of 18 wt % AOT in water was created for each CG model, which is well above the lamellar phase transition of aqueous AOT at 298 K.?

Because the mapping ratio is different for each model, so too is their molecular volume. During the solvation step, each bead was given the same van der Waals radius, so that the coarser models occupied a smaller volume and water was packed more densely around them in the initial configuration. This difference in density was corrected during the pressure coupling stage.

The number of molecules in the starting configurations are given in Table. It is important to note that the mass of the beads in the Martini parametrization is not necessarily equal to the total mass of the atoms that they are supposed to represent; the mass depends only on the size of the bead, i.e., how many non-hydrogen atoms it represents, regardless of their atomic masses.? The concentrations reported here were determined by converting the mole fractions using the real molar mass of AOT (M r = 444.6).

To simplify the analysis, the system was extended in the z-axis, so that it was more likely the bilayer would form in the xy-plane. Specifically, a cubic box (10 nm sides) of randomly dispersed AOT was solvated in a rectangular box of solvent (10 nm × 10 nm × 15 nm).

The simulation procedure was as follows:

• 1.An energy minimization stage, for 100000 steps or until the maximum force in the system was no greater than 100 kJ mol ^–1^ nm^–1^.
• 2.An NPT ensemble, using isotropic pressure coupling, for 50 ns (Δt = 20 fs). This was in order to allow a bilayer to form in the xy-plane. A stochastic cell rescaling barostat (τ = 4 ps, κ = 3 × 10^–4^ bar^–1^)? was used to couple the pressure.
• 3.A further 150 ns (Δt = 20 fs) NPT simulation using a semi-isotropic stochastic cell rescaling barostat (τ = 4 ps, κ_ xy _ = κ_ z _ = 3 × 10^–4^ bar^–1^).

A stochastic velocity rescaling thermostat was used to control the temperature? (τ = 1 ps, T = 298 K). We found that using only one thermostat was sufficient to control the temperature. This is likely because the concentration of surfactants is relatively large and therefore the difference in the number of degrees of freedom between the solute and solvent is less significant, so that the temperature of the solvent and solute remain the same.

During the semi-isotropic pressure coupling stage, the xy-area of the simulation box changes in dimension in order to minimize the in-plane tension of the bilayer. After this area had converged, two properties were measured: bilayer thickness and area-per-headgroup.

The Willard–Chandler surface of the bilayer was computed using the method described in the Dilute Isotropic Phase section. In this case there are two layers: one upper and one lower. The marching cubes algorithm used to identify the Willard–Chandler surface does not automatically partition the vertices into these layers.? For these calculations, the Pyvista software library? was used to identify which vertices belong in which layer based on their connectivity. This yields two surfaces (one for each layer).

The area per headgroup was calculated by dividing the combined area of the two surfaces by the number of surfactants, while the thickness was determined by calculating the nearest-neighbor distance between the two surfaces at each vertex.

Results and Discussion

Dilute

Isotropic Phase

Figure shows the normalized aggregation numbers for each dilute isotropic system as a function of simulation time. Each of the 1 wt % systems forms a monolithic aggregate after approximately 1 μs, while the 0.27 wt % systems exhibit a much more gradual growth and present a distribution of aggregation numbers at equilibrium.

Normalized aggregation numbers (fraction of total number of surfactants in a cluster) for each aggregate (aggregation number ≥5) in the dilute isotropic systems. Only aggregates with an aggregation number of at least 5 are shown.

It is important to note that the types of aggregates obtained differ depending on the surfactant concentration. For each model, the 1 wt % systems form a vesicle shortly after the surfactants have assembled into a monolithic aggregate, while the 0.27 wt % systems form micelles; see Figure for a snapshot of each of the Mixed model systems in their final configurations.

Final configurations of the Mixed dilute isotropic systems at (a) 0.27 wt % and (b) 1 wt % AOT. Yellow beads represent surfactant headgroups and gray beads represent tail groups. The box in 5­(a) is isotropic, of side ≈ 37 nm.

Several studies report that the critical vesicle concentration (CVC) of AOT occurs at 7.5 mM, which is approximately 0.33 wt %. ?,?,? By definition, above this concentration the surfactants spontaneously self-assemble into vesicles. Our results are consistent with these data; we observe micelles below the CVC and vesicles above it.

For completeness, it should be noted that some studies have used small-angle neutron scattering (SANS) to measure micelle properties at concentrations well above the CVC. ?−? ? Furthermore, prior all-atomistic molecular dynamics simulations found not vesicles, but micelles at 1 wt % AOT in water.? It is possible that micelles coexist alongside vesicles at these concentrations, though the results of Velázquez et al.? suggest that the population of micelles must be very small.

We also note that, as discussed in the Introduction, the dynamics of coarse-grained simulations are generally much faster than for all-atomistic simulations. For our Coarsest 1 wt % system, vesicle formation took approximately 600 ns. A typical scaling factor to convert the effective time from a coarse-grained to all-atomistic system is around 4,? which implies that the transition may take more than 2 μs to occur when all-atomistic simulations are used to simulate the self-assembly process. This is longer than the runtime used by Bhat et al.,? which may explain why they did not observe vesicle formation despite having considered concentrations above the CVC.

Figure shows the evolution of the properties of the aggregates in the 1 wt % systems. Visual analysis of these results suggests that, for our models, there is an increase in the time taken to form a vesicle as the degree of coarse-graining decreases. This is explained by the fact that the potential energy surface becomes rougher as the degrees of freedom increase. Therefore, the dynamics of the Finest model are much slower than those of the Coarsest model.

Mean surface area to volume ratio, surface area per surfactant, volume per surfactant, and aggregation number of aggregates (aggregation number ≥5) in 100 ns bins for each of the 1 wt % systems. Shaded area represents the 95% confidence interval calculated using bootstrapping.

The size and shape of the vesicles is very consistent between our models. The main difference is that there is a monotonic increase in the surface area and volume per surfactant as the degree of coarse-graining increases. This implies that the more coarsely grained models are unable to pack as densely, which we attribute to the larger average van der Waals radius of the beads. This is consistent with the results obtained for the bilayer simulations, as discussed later.

We note that the size of the vesicles is likely constrained by the box size and the number of surfactants in the simulation. In order to establish the equilibrium size of a vesicle, larger simulations must be performed. This will be the subject of a future paper. Some researchers find that the average hydrodynamic radius of AOT 10 mM vesicles is approximately 40 nm 70 nm,? while others find it may be as large as 150 nm.? Considering that vesicles are likely dynamic systems not necessarily at equilibrium, this difference in reported size of vesicles is not surprising. Nevertheless, it should be pointed out that these values exceed the size of the simulation boxes that can be considered using our current computational resources.

Notably, there is a brief jump in aggregate volume in Figure once the surfactants have formed a monolithic aggregate. This jump corresponds to a transitional period, before the aggregate becomes a vesicle. Analysis of the simulation snapshots at this point reveals the presence of disk-like structures, which are referred to as “bicelles” in the literature.? They are essentially like small bilayers, in that they are planar, but they do not percolate across the periodic boundaries of the simulation box.

During this transition, the bicelle begins to curve into a concave shape, until finally the edges meet to form a vesicle. This process is well described in the literature and has been observed experimentally for other surfactants ?−? ? ? but to our knowledge it is the first time this process has been demonstrated for simulations of AOT. This suggests that our CG models, coupled with the long simulation times implemented in this work, are appropriate for identifying key steps in the self-assembly of AOT. Snapshots from this transition are shown for the Mixed system in Figure and a video showing the full transition is available in the Supporting Information.

Stages of bicelle-to-vesicle transition for the Mixed model. The initial bicelle (a/e) begins to curve in on itself (b,c/f,g), eventually closing to become a vesicle (d/h). The top row shows the side view and the bottom row shows the top view at the same time step.

Figure shows the evolution of aggregate properties for the 0.27 wt % systems. These are much larger than the 1 wt % systems, and the growth process observed in our simulations is not monotonic; surfactants will often leave the aggregates. This highlights their metastable nature. Because of these competing processes, it takes far longer to reach equilibrium. The stochastic nature of the growth process suggests that it may be necessary to repeat the simulations several times in order to more reliably characterize the micelles. Furthermore, despite the fact that these simulations are quite large, they are still subject to finite size effects. Fajalia et al.? report that the mean intermicellar distance for this concentration of AOT in water is 223.2 Å, which is much larger than the maximum possible distance in our simulation boxes.

Mean surface area to volume ratio, surface area per surfactant, volume per surfactant, and aggregation number of aggregates (aggregation number ≥5) in 100 ns bins for each of the 0.27 wt % systems. Shaded area represents the 95% confidence interval calculated using bootstrapping.

At 0.27 wt %, although the Mixed and Finest systems’ properties appear to reach a plateau, with only one or two coexisting micelles, the Coarsest system exhibits a much greater degree of polydispersity and the properties fluctuate a lot more, even once converged. This is likely due to the faster dynamics that are characteristic of the Coarsest model, which facilitate surfactant exchange across shorter time scales. The large variance in aggregate properties across the equilibrium regime is also explained by the metastable nature of the aggregates. In any case, the simulation results show that, qualitatively, the three models yield very consistent results.

The resulting time-averaged mean property values and their standard deviation are reported in Table. The experimental value for the radius of gyration, R g, was calculated using the experimentally determined ratio of the major to minor semiaxes, ε = 1.22:?

$$[eqn]$$

2: Time-Averaged Mean Properties of the Systems (Averaged across Aggregates with Aggregation Number ≥5) in the 0.27 Wt % Simulations Compared to Experimental Values

For consistency, eq employs the same assumptions made by Fajalia et al:? that the micelles are oblate ellipsoids (a ≈ b ≫ c), with a semiaxis c = 12.57 Å, corresponding to the length of a fully stretched alkyl tail. Fajalia et al.? used these assumptions to infer the volume of their micelles. Our results suggest that these assumptions may not be accurate in all conditions, as our values for CPE indicate that the micelles are predominantly prolate or triaxial ellipsoids. Care must be taken, therefore, when comparing the computational values of R g and the surface area to volume ratio to the experimental ones. For completeness, we note that in both of these cases, the Coarsest model is in closer agreement with the experimental values. This correlates with the fact that the Coarsest model exhibited the largest average aggregation number.

The surface area per surfactant for the Mixed and Finest models seems to be in very good agreement with the experimental value of 194.3 Å^–2^. We consider the CG models derived here reasonable approximations of reality, although one should be careful in extrapolating results obtained for the size of AOT micelles.

Finally, Figure shows the coordinate-pair eccentricities (CPEs) of every aggregate across the trajectories of the simulations. The vesicles, which correspond to a normalized aggregation number of 1, are all mostly spherical. Conversely, micelles are almost always prolate ellipsoids. Generally, larger micelles have a greater e _ ac _ than smaller micelles. The transitional bicelles can also be identified in the CPE plots for the 1 wt % systems: they are the points with a normalized aggregation number of 1 that exhibit a very high e _ ac _ value and an intermediate e _ ab _ value, indicating that they are triaxial ellipsoids.

Coordinate-pair eccentricities of every aggregate (aggregation number ≥5) across the entire trajectory for each of the coarse-grained simulations. See also Figure for a schematic interpretation of this mapping.

Bilayer Properties

To assess the convergence of the bilayer simulations, the xy-area of the simulation box was monitored. After 30 ns, all of the systems’ areas had converged. Notably, there was a monotonic increase in the equilibrium area and the time taken to reach it as the degree of coarse-graining increased. This is due to the fact that the coarser systems were initialized with more AOT molecules. The timeseries is shown in the Supporting Information.

The distributions of the bilayer thickness measurements and the area per headgroup were calculated for time t > 30 ns for each system. The results are shown in Figure and the mean values are reported in Table. The thicknesses obtained from the simulations using the Finest and Mixed models are quite similar; both are ∼5% thinner than the experimental value. The Coarsest model bilayer, however, is over 10% thicker on average than the experimental data.

Combined kernel density estimate and box-and-whisker plots showing the distributions of (a) thickness and (b) area per headgroup of the AOT bilayer for each CG model. The distributions are aggregated across t > 30 ns. Dashed lines indicate experimental values. Experimental values were derived by Fontell using X-ray diffraction.

3: Mean Bilayer Thickness and Area Per Headgroup for Each CG Model

There is a more significant difference in the area per headgroup. In this case, the Mixed and Coarsest models are both within ∼10% of the experimental value. The Finest model’s mean area per headgroup is ∼13% lower than the experimental value. The area per headgroup exhibits a monotonic increase with respect to the degree of coarse-graining. This mirrors the trend observed for the dilute isotropic systems, whereby the more coarsely grained models are unable to pack as densely.

Based on these results, the Mixed model seems to be the most consistent with the experimental data, though all of the models derived here are reasonably accurate. This is in spite of the fact that the Finest model is “closer” to the atomistic model in terms of the number of degrees of freedom. This may be due to the way that Martini was parametrized: the finer degrees of mapping are only recommended when necessary, at the end of chains or in the case of branched chains, where a 4-to-1 mapping represents too large of a group.?

Conclusions

In this work, three coarse-grained models for AOT were developed with different degrees of mapping from the atomistic scale: Finest, Mixed and Coarsest. The models were validated by simulating systems of different concentrations; two systems (0.27 and 1.0 wt %) below the solubility limit of AOT in water, which represent the dilute isotropic phase, and one (18 wt %) which represents an anisotropic phase, in order to study the properties of self-assembled bilayers.

It was found that all of the models were able to broadly reproduce the expected self-assembly behavior of AOT in water, forming micelles at 0.27 wt % (below the critical vesicle concentration) and vesicles at 1.0 wt % (above the CVC). The 0.27 wt % system was characterized by very slow dynamics, with the Finest model requiring up to 13 μs to converge. The systems at higher concentrations converged faster. Of note, the simulations at 1.0 wt % showed, for the first time, in the case of AOT, the transition from the large, disc-like micelles known as “bicelles” to vesicles. The vesicles were found to be spherical, while the micelles were prolate ellipsoidal.

There are two principal differences in the behavior of the models. First, the more coarsely grained models are unable to pack as densely, which is reflected in the surface area per surfactant and the volume of the vesicles, as well as the area-per-headgroup of the bilayers. Second, the dynamics of the models are faster for the more coarsely grained systems, as evidenced by the rate of micelle growth.

The Mixed model appears to be the most reliable choice based on the agreement between our simulations and experimental bilayer properties. This model promises a new way to interrogate the behavior of AOT vesicles on the molecular scale and, simultaneously, to study their effect on macroscopic behavior like droplet bouncing and spreading. This will be the subject of future work.

Supplementary Material