Tristearin bilayers: structure of the aqueous interface and stability in the presence of surfactants

Abstract

We report results of atomistic molecular dynamics simulations of an industrially-relevant, exemplar triacylglycerol (TAG), namely tristearin (TS), under aqueous conditions, at different temperatures and in the presence of an anionic surfactant, sodium dodecylbenzene sulphonate (SDBS). We predict the TS bilayers to be stable and in a gel phase at temperatures of 350 K and below. At 370 K the lipid bilayer was able to melt, but does not feature a stable liquid–crystalline phase bilayer at this elevated temperature. We also predict the structural characteristics of TS bilayers in the presence of SDBS molecules under aqueous conditions, where surfactant molecules are found to spontaneously insert into the TS bilayers. We model TS bilayers containing different amounts of SDBS, with the presence of SDBS imparting only a moderate effect on the structure of the system. Our study represents the first step in applying atomistic molecular dynamics simulations to the investigation of TAG-aqueous interfaces. Our results suggest that the CHARMM36 force-field appears suitable for the simulation of such systems, although the phase behaviour of the system may be shifted to lower temperatures than is the case for the actual system. Our findings provide a foundation for further simulation studies of the TS-aqueous interface.

---

Investigation of the structure and properties of fats has a long and rich history, due to the importance of such molecules both biologically and industrially.^1,2 Fats and oils comprise a complex mixture of triacylglycerols (TAGs), also known as triglycerides. TAGs consist of a tri-ester of glycerol with fatty acid tails, and are neutral lipids insoluble in water.^1,2 TAGs have been proposed for use in the future implementation of solid–lipid nanoparticles (SLNPs) for nanomedicine applications.^3–6 SLNPs have been proposed as delivery agents for drugs, imaging agents and biomolecules. In these systems the drug is encapsulated within the hydrophobic core of the SLNP and are slowly released as the SLNP is broken down. Like polymeric nanoparticle (NP) drug delivery systems,^7,8 SLNPs could offer a biodegradable route for the sustained and targeted release of therapeutics, but with the added advantage of being constructed from non-toxic, physiological molecules.

Naturally-occurring TAGs are generally heterogeneous, possessing fatty acid tails of different lengths and unsaturations.² However, for convenience, both experimental and theoretical past studies have often used the homogenous TAG molecules as model systems. Below the solid–liquid transition temperature TAG molecules form crystalline networks of nanoplatelets.^2,9,10 A schematic illustration of the structure of a fat nanoplatelet is shown in Fig. 1(a). Within the nanoplatelets the TAG molecules are arranged in lamellar structures with the acyl tails aligned approximately perpendicular to the planes formed by the glycerol head groups.^2,11,12 In these lamellae the TAG molecules generally adopt one of two conformations, either the “tuning fork” (where the sn-1 and sn-3 chains—see Fig. 1(b)—have the same orientation relative to the glycerol plane) or the “chair” (where the sn-1 and sn-2 chains have the same orientation relative to the glycerol plane). Another possible conformation is the “trident”, where all three acyl tails pack side by side, which while not predicted to be prevalent in bulk TAG crystals, could play an important role at hydrophilic–hydrophobic interfaces.^2,13–15

Fig. 1 (a) A schematic of the structure of the aqueous interface of a fat nanoplatelet. The chemical structure of (b) tristearin (TS) and (c) sodium dodecylbenzene sulphonate (SDBS). The sn-1, sn-2, and sn-3 acyl tails of TS are marked. Also labelled are the C1, C9 and C18 positions in the acyl tail of TS and the C1 and C12 position of the acyl tail of SDBS.

The aqueous interface of TAG systems is of particular interest to the detergent industry, as the removal of crystallised lipids from soiled fabrics is an energy-intensive process, usually requiring high temperatures even with modern detergents.^16–19 As such, there is considerable interest in developing formulations that are able to break down and dissolve TAG meso-structures at or close to room temperature in aqueous solution. Experimental studies have investigated using nonionic surfactants in combination with a small amount of a lipophilic amphiphile to enhance the removal of lipids at lower temperatures.^20,21 In addition, a recent study reported that the presence of nanodiamonds (NDs) alongside various surfactants could assist in the removal of crystallised lipid for both anionic and nonionic surfactants at low temperatures (15–25 °C), especially in the case of the anionic surfactant sodium dodecylbenzene sulphonate (SDBS).¹⁹ The increase in TAG removal efficiency rendered by the ND–surfactant mixtures was proposed to relate to the ability of NDs to roughen the surface of the lipid particles, enhancing surfactant adsorption.

While there are still a number of open questions regarding the mechanism of removal for the crystallised lipid, we suggest that further investigation of the molecular-level interactions of NPs, including NDs, and other molecules at the aqueous-TAG interface is warranted. However, while experimental techniques can provide valuable information about the effects of molecules on TAGs, it is often challenging to determine the underlying structures of these systems at an atomistic level. Thus, molecular simulation can play an important role in providing complementary insight into the interaction of NPs, including NDs, at the aqueous interface of TAGs. Investigations of the bare aqueous TAG interface (in the absence of additives) forms a first essential step in this realising this goal.

Simulation and theoretical studies have been provided useful information about TAG systems at the mesoscale, elucidating the phase behaviour, structure and aggregation of TAGs nanoplatelets.^9,22,23 At the atomistic level, molecular dynamics (MD) simulations of TAGs have investigated the behaviour of bulk phase TAGs; their lamellar structures and the solid–liquid phase transition.^{11,12,23–25} In general, the results of these past studies have been in agreement with both theory and experimental studies, despite the challenging nature of these simulations, which require long time- and length-scales. Below the solid–liquid transition, the formation of lamellar phases are observed, with the majority of TAG molecules adopting either the tuning-fork and/or chair conformations. In addition, recent work has investigated the solid–liquid/fat–oil interface of TAGs.^26,27

While these studies have shown the ability of molecular simulation to help elucidate the structure of TAGs at the atomistic level, to the authors' knowledge, there have not been any previous MD studies investigating TAG-aqueous interfaces. Therefore, prior to investigating the effect of NPs on TAG surfaces, the characterisation of the TAG-aqueous interface itself is required, in particular the ability of existing force-fields (FFs) to model such an interface must be evaluated. In this work we have investigated the behaviour of bilayers of a model TAG, tristearin (TS), at an aqueous interface, in both the presence and absence of SDBS molecules. The chemical structures of both TS and SDBS are shown in Fig. 1(b) and (c). TS has been used experimentally as a model fat in the investigation of the ability of different processes to remove laundry soil.¹⁹ The actual system would contain fat nanoplatelets, comprising multiple layers of TS molecules in a lamellar structure,^2,9,23 with individual TS molecules in the bulk assuming “chair” or “tuning fork” configurations, and with the TS molecules at the interface present in the “trident” configuration. To simulate such a system in its entirety at the atomistic level for the timescales considered here would demand extremely large computational resources. As such, we believe that a bilayer of TS molecules provides a reasonable first approximation for the structure that such nanoplatelets may present to the aqueous interface. We have characterised the structure of the TS bilayers at both ambient and elevated temperatures. Following this characterisation, we also have investigated the impact on the structure of the bilayers of the presence of SDBS molecules, at both temperatures.

Methods

As discussed previously, to the authors' knowledge, this is the first study to use MD simulations to investigate the behaviour of the TAG-aqueous interface and the interaction of TAGs with surfactant molecules. Because of this, our chosen FF must be able to capture both the TS and SDBS and their mutual interactions, in the presence of water. MD simulations have been widely used to investigate surfactant systems, however, sodium dodecyl sulphate (SDS)^28–34 or hexadecyltrimethylammonium bromide (CTAB)^29,35–38 have been investigated more extensively than sodium dodecylbenzene sulphonate (SDBS).^33,39,40 Part of the reason for this may be ascribed to the very low critical micelle concentration (CMC) of SDBS, ≤3.0 mM under standard conditions,^41,42 compared with ∼8.3 mM for SDS,⁴² meaning that the study of SDBS via simulation typically requires the simulation of systems comprising a large number of water molecules. He et al. developed a set of parameters within the CHARMM FF to model linear alkyl benzene sulphonates (LAS), the family of compounds of which SDBS is a member, in solution.³⁹ This study also investigated how the behaviour of alkyl benzene sulphonates varied with tail length and branching point. To this end, it was decided to model the TS and SDBS molecules via the CHARMM36 FF parameter set,^43,44 with additional parameters taken from the previous study of SDBS.³⁹ The modified version of the TIP3P^45,46 water model was used for the water molecules.

The simulations were performed using the GROMACS code, version 4.6.3.⁴⁷ The LJ-nonbonded interactions switched to zero between 0.9 and 1.0 nm and the electrostatics evaluated using a PME (cutoff 1.1 nm). The simulations were performed in the isothermal–isobaric (NpT) ensemble, with the Nosé–Hoover thermostat^48,49 and Parrinello–Rahman barostat⁵⁰ used to maintain the temperature and pressure, respectively. The bilayer systems were arranged such that the bilayer was parallel to the xy-plane, and the xy and z cell dimensions allowed to vary independently of one another, with a reference pressure, P_ref, of 1 bar. During the simulations the P_xx, P_yy and P_zz components of the pressure tensor were monitored to ensure that P_xx = P_yy = P_zz = P_ref (within statistical variation) and that no artefacts were introduced into the system from the pressure coupling. For all systems a time-step of 1 fs was used. In all simulations summarised herein, we saved configurations at 1 ps intervals.

The initial configurations were generated by taking four different TS molecules in trident conformations and then replicating this set of molecules three times in both the x and y axis to make one leaflet comprising 36 TS molecules. This monolayer was then rotated 180° about the x axis to create the opposite leaflet. After energy minimisation of the bilayer, a brief (10 ps) MD simulation was carried out at 350 K in vacuum, to remove some of the initial order of the system. This procedure was repeated three times to generate three different TS bilayer samples. Table 1 summarises the details of the different systems considered in our study. Sets A and B correspond to the simulations of the aqueous TS bilayers alone (i.e. in the absence of any SDBS molecules). For these runs, after the initial setup, the resulting bilayers (three samples) were each solvated in 4989 water molecules, amounting to ∼70 water molecules per lipid, and subjected to 300 ns of simulation at 350 K. These 300 ns trajectories at 350 K correspond to the Set B simulations. After this stage, each of the three bilayers was then cooled from 350 K to 300 K over a 100 ns timescale, and was followed by a further MD simulation at 300 K for another 300 ns. These 300 ns trajectories at 300 K correspond to the Set A simulations.

Table 1 Parameters of the different TS bilayers simulated: run label, run temperature, the number of SDBS molecules present in the simulation (values in parentheses are indicate the number initially inserted into the bilayer), N_SDBS, the lateral area, A, thickness of the bilayer, D_HH, thickness of the interface, D_Inter and averaged height two peaks, H_max

Run	Temp./K	NSDBS	A/nm^2	DHH/nm	DInter/nm	Hmax/kg m^−3
A1	300	0	21.4 ± 0.1	4.58	0.36/0.36	1305
A2	300	0	21.1 ± 0.2	4.60	0.32/0.34	1357
A3	300	0	21.4 ± 0.1	4.43	0.38/0.38	1300
B1	350	0	23.2 ± 0.3	4.40	0.45/0.45	1196
B2	350	0	23.3 ± 0.3	4.37	0.47/0.44	1183
B3	350	0	23.4 ± 0.3	4.28	0.46/0.46	1165
C1	300	20	21.4 ± 0.1	4.53	0.35/1.07	1314
C2	300	20	21.0 ± 0.1	4.53	1.38/0.37	1304
C3	300	20	21.5 ± 0.1	4.60	0.37/0.42	1326
D1	350	20	23.2 ± 0.2	4.46	0.48/0.48	1147
D2	350	20	23.1 ± 0.3	4.39	0.45/0.43	1207
D3	350	20	23.4 ± 0.3	4.44	0.47/0.47	1160
E1	350	28 (8)	23.1 ± 0.3	4.41	0.48/0.47	1120
E2	350	28 (8)	23.1 ± 0.3	4.41	0.47/0.47	1166
E3	350	28 (8)	23.4 ± 0.3	4.30	0.46/0.46	1146
F1	350	28 (16)	25.0 ± 0.7	4.19	0.46/0.52	1079
F3	350	28 (16)	25.5 ± 0.3	4.24	0.53/0.57	1045
G1	350	28 (24)	25.9 ± 0.3	4.09	0.55/0.50	1013
G3	350	28 (24)	25.7 ± 0.3	4.12	0.55/0.53	1006

---

The simulations in Sets C–G, correspond to those in which SDBS was present. The low CMC of SDBS meant that to mitigate the immediate formation of SDBS micelles in these samples, it was necessary to increase the number of water molecules in our system. Ultimately, it was found that modelling the TS bilayers with 20 or 28 SDBS and ∼13 600 water molecules, plus 20 Na⁺ counter ions, provided a reasonable compromise between a sufficiently low SDBS concentration (∼0.08 mol kg⁻¹) that would allow some SDBS molecules the opportunity to diffuse to the TS bilayer and also provide a TS : SDBS ratio that can facilitate sufficient statistical detail, while economising the computational expense.

For the simulations in Sets C and D the initial bilayer configurations were taken from the final configurations of the simulations of the runs of Sets A and B, i.e. the TS bilayers in the presence of water alone in the absence of any SDBS molecules. To construct the initial configurations for Sets C and D, the SDBS molecules were placed randomly in the space between the TS bilayers in the periodic simulation cell, before each system was solvated. In other words, in Sets C and D, SDBS molecules were present in the samples, but the surfactants were not embedded into the TS bilayer from the outset. Of the six resulting samples, the three samples derived from Set A were then subjected to 300 ns at 300 K, generating Set C. Similarly, the samples derived from Set B were run for 300 ns at 350 K, generating Set D.

To investigate the influence of different concentrations of SDBS on the bilayers in Sets E, F and G, mixed TS–SDBS bilayers were created containing 8, 16 or 24 SDBS molecules, respectively. However, the total number of SDBS molecules in each sample (both in solution and within the bilayer) was kept constant at 28 molecules. Because we had observed very few instances of spontaneous insertion of SDBS molecules into the bilayers at 300 K (see Results and discussion), all of these simulations were performed at 350 K for 200 ns.

In addition to the simulations of the bilayers, control simulations of SDBS in water (in the absence of the TS bilayer) were carried out at both 300 and 350 K. These control systems comprised 20 SDBS molecules solvated in 13249 water molecules and 20 Na⁺ counter ions. The SDBS molecules were initially randomly distributed in a cubic cell, and subjected to 100 ns of MD simulation in the NpT ensemble with isotropic pressure coupling. All the remaining simulation parameters were the same as those used for the TS bilayer and TS bilayer–SDBS simulations.

For the bilayer systems the final 60 ns of the trajectory was used for analysis, while analysis of the SDBS control simulations was performed over the final 20 ns.

Results and discussion

Tristearin bilayers

The lateral area, A, membrane thickness D_HH (measured as the distance between the two peaks of the TS mass density profile), interface width D_Inter (measured as the distance over which the water density decays from 90 to 10% of the bulk value), and peak height H_max (measured as the average of the density maximum of the two leaflets in the TS density profile) of the different TS bilayers are given in Table 1. On the 300 ns timescale considered here the bilayers are stable at both temperatures and the systems appear to have equilibrated. As an example, Fig. S1 in the ESI† shows the evolution of A as a function of time for one run at 300 and 350 K. Fig. 2 shows the area per lipid, A_L for each of the TS bilayers in Sets A and B, as well as the experimentally measured A_L of a TS monolayer at the air–water interface at 293 K (0.62 ± 0.01 nm²).¹⁵ There is no statistically significant difference in A_L amongst the three 300 K runs; the same is also true for the 350 K runs. As one would expect, the bilayers at 350 K have a greater A and lower D_HH than the bilayers at 300 K, i.e. the bilayers have expanded in the plane parallel to the bilayer and have contracted in the direction perpendicular to the bilayer. Compared to the bilayers at 300 K, those at 350 K showed an increase in A of ∼9% and a decrease in D_HH of ∼4%. The average A_L value obtained from simulations of bilayers at 300 K, 0.59 nm², is slightly lower than the experimental value but not unreasonably so.

Fig. 2 Area per lipid, A_L, of the TS bilayers for all three bilayers at 300 and 350 K, compared against the experimental value for a TS monolayer.

Exemplar density profiles of bilayers at 300 and 350 K are shown in Fig. 3. The classical lipid bilayer profile is observed, with peaks in the lipid profile at the water–lipid interface corresponding to where the head groups were located, and a minimum at the bilayer centre, indicating that there was little interdigitation between the tails in the two leaflets. As discussed above, at 350 K the TS bilayers were thinner than those at 300 K and the height of the peaks in the TS profile was decreased. Despite these differences, the overall shape of the profiles is very similar at the two temperatures, indicating a lack of significant change in the structure of the bilayers at these two temperatures.

Fig. 3 Exemplar density profiles of the TS bilayers at (a) 300 and (b) 350 K.

As can be seen from Fig. 4(a) the overwhelming majority of TS molecules were present in a trident type of conformation with the glycerol head group at the TS-aqueous interface, in agreement with the experimental evidence.^13–15 At both temperatures the TS bilayers are in a gel phase with the acyl tails of the lipids packing in a hexagonal pattern, as can be seen in Fig. 4(b). The same packing behaviour of acyl tails is observed in simulations of ceramide bilayers.^51–56 While nearly all of the TS molecules are in a trident conformation, a number of different variations in this conformation were noted. The majority of TS molecules are arranged into a tripod type structure, such that the angle ∠T₁T₂T₃ (where T_x refers to the position in the x/y-plane of the sn-x tail) was approximately 60°. However, some TS molecules assumed a linear arrangement, ∠T₁T₂T₃ ≈ 180° or a kinked arrangement, ∠T₁T₂T₃ ≈ 120°, as highlighted in Fig. 4(b). In addition, there were a few molecules that did not fall into one of these three possible categories, effectively becoming “defects” within the hexagonal packing arrangement. The relative percentage of the three categories was similar at 300 and 350 K, being; 46 (42)%, 20 (18)% and 28 (26)% at 300 K (350 K) for the tripod, linear and kink categories, respectively. However, the systems at 300 K had fewer molecules in “defect” configurations than those at 350 K. This may be a result of either the annealing period the 300 K system underwent (see Methods) and/or the greater density of the bilayer forcing defects from the system.

Fig. 4 Snapshots taken from simulations of TS bilayer at 300 K. (a) The tristearin bilayer, tristearin oxygen atoms, carbon atoms and water molecules are coloured red, grey and cyan, respectively, the tristearin hydrogen atoms are not shown for clarity. (b) A cross-section of the bilayer showing the carbon atoms in the acyl tails of one monolayer. The acyl tails pack in a hexagonal arrangement. Examples of the tripod, linear and kink configurations are coloured orange, blue and yellow, respectively.

As evidence of the lateral ordering in the TS bilayer leaflets, in Fig. S2 in the ESI† we show the g(r) in the x/y-plane (i.e. the lateral plane) for the first and ninth carbons of the acyl chains within one leaflet. From the position of the peaks in this plot it is possible to determine the short-ranged “lattice constant” for the hexagonally packed tails, which is ∼0.46 nm and ∼0.48 nm at 300 and 350 K, respectively. In contrast to the well-packed tails, the head groups of the TS molecules showed little ordering. The lateral g(r) of the TS head groups (defined as centre of mass of the three carbons in the glycerol head group) is shown in Fig. S2(b).† The g(r) plots show that head groups are typically located ∼0.8 nm from each other, but with a few head groups at a closer separation of ∼0.6 nm.

The degree of ordering in the acyl tails of the phospholipids can be determined from the order parameter, S_z. S_z for atom C_n is calculated from

where θ is the angle between the z-axis and the vector between C_n−1 and C_n+1. The order parameters for the carbon atoms in the lipid tails for the sn-1 chain are shown in Fig. 5, the profiles for the sn-2 and sn-3 are similar (see Fig. S3 in the ESI†). At 300 K the order parameters for the C6–C13 carbons show little variation with position. The order at start of the chain (C2–C4), where the presence of the aqueous interface may perturb the ordering of the alkyl groups, and the end of the chain (C16–C17), where the alkyl groups have greater conformational freedom due to the lower density, is lower than seen for the central alkyl groups. At 350 K there was a general decrease in the value of the order parameters, but with the shape of the order parameter profile showing only slight changes (the profile became more curved) compared to those of the systems at 300 K.

Fig. 5 The order parameter of the sn-1 and acyl tails at 300 and 350 K.

It is of interest to compare the behaviour of the TS bilayers to ceramide 2 bilayers, as they are similar to each other in other ways. Both types of lipids have long fatty acid tails with relatively small polar headgroups, compared with phospholipids or even sphingomyelin. This gives rise to a narrow lipid–water interface width, 0.42–0.54 nm and ∼0.35 nm at 300 K for ceramide 2 (ref. 51) and TS, respectively. In addition, the hexagonal packing of the lipid tails in both systems have the acyl tails tilted at a slight angle relative to the bilayer normal. In the case of the ceramide bilayers this tilt angle is ∼17–24° (depending on the lipid modelled and FF used), while for the TS bilayers we predict a tilt angle of ∼12°.

Simulations of ceramides have indicated that the bilayer can undergo a phase change with the acyl tails no longer packing in a hexagonal arrangement but being essentially liquid.^51,56 In the case of TS, our simulations predicted the bilayer to be stable at 350 K. However, when the temperature was increased to 370 K, a phase change was observed. However, unlike in the case of ceramides, which retain a stable bilayer structure (at least on the timescales simulated),^51,56 once the lipids were in the liquid phase, the TS bilayer became unstable. To investigate this phase transition in more detail we performed annealing runs on each of the three bilayers simulated at 350 K. Similar simulations have been reported for phosphocholine bilayers to investigate the gel–liquid crystalline phase transition.^57,58 To accomplish this, we took the final configuration from each of our Set B runs, and heated these samples from 350 to 370 K at a rate of 0.1 K ns⁻¹ (i.e. corresponding to a further simulation duration of 200 ns). The lateral area of the bilayers as a function of time/temperature is shown in Fig. 6 (snapshots of the bilayer during the course of this simulation are shown in the ESI, Fig. S4†). In general, these profiles show that in the gel phase, A increases linearly with temperature, with a sharp increase in A upon transition to the liquid crystalline phase, as the acyl tails of the TS molecules melt (ESI, Fig. S4(a)†). Here, we define the transition temperature as the temperature average taken over the last 10 ps of trajectory in the gel phase and the first 10 ps of trajectory in the liquid crystalline phase. We identified the transition temperature to lie in the range of ∼365.5–367 K (Fig. 6). However, after only a short time above the transition temperature, the bilayer structure becomes unstable as liquid TS molecules move away from the bilayer–water interface into the centre of the bilayer. This leads to a sharp contraction in the lateral area of the “bilayer” (ESI, Fig. S4(b)†), which, due to the periodic boundary conditions imposed on it, drives the system towards a configuration where x/y → 0 and z → ∞. As TS molecules diffuse away from the aqueous interface, a few water molecules are able to enter the bilayer. The head groups of the TS molecules pack around these water molecules, such that there are regions within the TS system that are dense in head groups (see Fig. S4(b) in the ESI†). While such behaviour should not be used to infer information about the real TS system at elevated temperature, these data provide useful information regarding the force-field used to describe these aqueous TS systems.

Fig. 6 The A of the TS bilayers as a function of temperature/time. Dashed lines indicate the transition temperature for each sample (see text for definitions).

The results of the simulations of TS bilayers are in general agreement with available experimental data, as well as previous simulations of similar lipid systems. However, the paucity of experimental data at the atomistic level of resolution makes it difficult to determine how closely the modelled bilayers represent the actual system. The fact that A for the bilayers at 300 K is lower than the experimental value, though not unreasonably so, and that the bilayers are stable at 350 K and do not break down until 370 K (experimentally, the TS nanoplatelets will start to break apart at around 350 K (ref. 19)), suggests that the CHARMM36 FF may over-stabilise the TS bilayers slightly with respect to temperature. This in effect artificially “shifts” the phase behaviour of the system to higher temperatures for this FF. However, further studies, both experimental and theoretical, are needed before this behaviour can be definitively resolved.

Tristearin bilayers in the presence of SDBS

Table 1 gives the A, D_HH, D_Inter, H_max of the different simulations of TS bilayers in the presence of SDBS molecules (Sets C, D, E, F and G). In the control simulations comprising SDBS solvated in water, the SDBS molecules aggregated and formed a single agglomerate at both 300 K and 350 K. The solvent accessible surface area of the agglomerates was 55 ± 1 nm and 68 ± 3 nm at 300 and 350 K, respectively. In the simulations for Sets C and D, where SDBS molecules where randomly distributed in solution, the SDBS molecules showed a propensity to aggregate, either in solution or at the surface of the lipid bilayer (Fig. S5 in the ESI† shows an exemplar snapshot from one simulation, illustrating the formation of aggregates of SDBS molecules both in solution and at the interface). In addition, however, some SDBS molecules also spontaneously embedded within the TS bilayer. Fig. 7 shows an example of a TS bilayer structure where a number of SDBS molecules had spontaneously entered the bilayer structure. The number of SDBS molecules embedded within the bilayer as a function of time for the Set C and D simulations is shown in Fig. S6, ESI.† The number of embedded SDBS molecules is greater at 350 K than 300 K, probably due to the fact that the bilayer was not packed as tightly at 350 K as it was at 300 K, thus facilitating easier ingress of the SDBS molecules into the bilayer structure.

Fig. 7 Snapshot of a TS bilayer into which SDBS molecules have spontaneously inserted into the bilayer (taken from run D1). The carbon and oxygen atoms of TS are coloured grey and red, respectively, the SDBS molecules are green and the water is coloured cyan.

While the simulations corresponding to Sets C and D showed that SDBS molecules were able to spontaneously insert into the bilayer at the TS-aqueous interface, the simulation timescale over which the equilibrium distribution of bilayer-inserted SDBS molecules could be determined is not practicable for standard atomistic MD simulations. Therefore, to better understand the effect of SDBS on the TS bilayer structure, mixed SDBS–TS bilayers were constructed, initially containing 8, 16 or 24 SDBS molecules. The simulations corresponding to these mixed bilayer systems comprises Sets E, F and G, respectively. The total number of SDBS molecules in the simulation cell was kept constant at 28 molecules. There were no restrictions placed on the remaining solution-phase SDBS molecules from entering the mixed bilayer, nor were there restrictions on any pre-inserted SDBS molecules from exiting the mixed bilayers. Fig. 8 shows the lateral area of the different systems averaged over all the bilayers in that set of simulations. The bilayers showed no appreciable increase in their lateral area when eight SDBS molecules were initially incorporated. However, a lateral expansion of the bilayers containing sixteen or more SDBS molecules was noted. This expansion was also accompanied by a transverse thinning of the bilayer and an increase in D_Inter (Table 1).

Fig. 8 The average lateral area of the bilayers at 350 K containing 0, 8, 16 and 24 SDBS molecules, corresponding to simulation Sets B, E, F and G, respectively.

The SDBS molecules embedded within the bilayers had their hydrocarbon tails located within the body of the TS bilayer. It appeared less favourable for the bulky phenyl rings to enter the tightly packed structure of the bilayer, and these were generally left somewhat exposed to the solvent. This means that the peak in the SDBS density profile, shown in Fig. 9 for run G1 (with 24 SDBS initially placed in the mixed SDBS–TS bilayer) was positioned further from the centre of the bilayer compared with that of TS. The density profiles of all simulations in Sets C, D, E, F and G are provided in the ESI, Fig. S7–S11.† While the incorporation of SDBS molecules into the TS bilayers does cause changes in the density profiles of TS, with the height of the TS peaks being reduced and the two leaflets sometimes developing an asymmetry, the overall nature of the profile remained largely similar. Even at a TS : SDBS ratio of 3 : 1 the bilayers showed no evidence of structural instability at 350 K, over the 200 ns timescale considered.

Fig. 9 The density profile of a bilayer containing 24 SDBS molecules (run G1).

Despite the addition of the SDBS molecules into the bilayer, the hexagonal packing arrangement of the fatty acid tails remained intact. Fig. S12 of the ESI† shows a cross-section through a single leaflet of two bilayers, containing 16 and 24 SDBS molecules. Overall, the hexagonal lattice did not appear particularly perturbed, however, there did appear to be some localised packing defects around the SDBS molecules, particularly in the case where a number of SDBS molecules have clustered together.

The average order parameter profiles of the acyl chains of the TS molecules in the mixed SDBS–TS bilayers are shown in Fig. S13 in the ESI.† There is a slight decrease in the ordering of the tails as the amount of SDBS in the bilayer is increased, but the systems remained highly ordered. The order parameters of the carbon atoms in the SDBS tail, for those SDBS molecules embedded in the bilayer are presented in Fig. S13(d), ESI.† The overall ordering of the SDBS acyl tails was found to be on the same level as that of the TS molecules. However, the top and bottom of the SDBS tails appeared more ordered than the TS molecules. This can be explained by the fact that the SDBS molecules are significantly shorter than the TS molecules (a hydrocarbon chain of 12 carbon atoms as opposed to 18 carbon atoms), meaning that the ends of the SDBS molecules are more constrained, resulting in a higher degree of ordering.

Table 2 gives the radius of, and coordination number of water molecules within, the first solvation shell of the sulphur and phenyl ring of SDBS molecules at 350 K. These data have been calculated for SDBS molecules in different environments: those that were located in the agglomerate in the control (SDBS–water) simulations, those located in the aggregates that were formed in the mixed SDBS–TS bilayer simulations (data taken from Set E), and those SDBS molecules that were initially embedded in the mixed SDBS–TS bilayers (simulations in Sets E, F and G). The radial distribution functions for the control and embedded cases are shown in Fig. S14 of the ESI.† The values in Table 2 for the control agglomerates and the Set E aggregates (i.e. n_Con and n_Agg-E, etc.) were similar. However, upon insertion into the bilayer, both the sulphonate head group and the phenyl ring were found to be desolvated with respect to the control and Agg-E cases. As the number of SDBS molecules incorporated in the bilayer was increased from 8 to 16/24 a further small desolvation effect was noted. However, this desolvation did not affect the number of hydrogen bonds formed between the sulphonate head group and water molecules. The number of SDBS–water hydrogen bonds formed per SDBS molecule was calculated to be 5.5 ± 0.1, 5.4 ± 0.1 and 5.1 ± 0.1 for the SDBS in the control micelles, Set E aggregates, and when embedded in the bilayer (Sets E, F and G), respectively. The number of SDBS–water hydrogen bonds formed by the embedded SDBS molecules calculated for Sets E, F and G was found to be similar.

Table 2 Radius of the first solvation shell, r, for the sulphur atom and phenyl ring of SDBS molecules, and the coordination number, of water molecules for the SDBS molecules located in aggretates in the control run, n_Con, those located in aggregates in mixed SDBS–TS bilayer simulations (Set E), n_Agg-E, and those embedded within the mixed SDBS–TS bilayers, n_Emb-X, where X = E, F and G. All data were generated at 350 K

	S-water	Ring-water
r/nm	0.48	0.38
nCon	11.19 ± 0.05	1.95 ± 0.03
nAgg-E	11.15 ± 0.04	1.93 ± 0.02
nEmb-E	9.59 ± 0.08	1.09 ± 0.04
nEmb-F	8.96 ± 0.14	1.04 ± 0.04
nEmb-G	9.19 ± 0.13	1.08 ± 0.01

---

The calculation of diffusion coefficients for lipid bilayers is a much discussed topic, even when the bilayer is in the liquid crystalline phase.^59–61 As TS bilayers are only stable in the gel phase, the diffusion of the lipids is extremely slow⁶² which means it is not possible, on the 200–300 ns timescales considered in our study, to determine if the SDBS molecules have any meaningful effect on the diffusion of the lipid molecules.

To investigate if the incorporation of SDBS molecules into the bilayer reduced the phase transition temperature, annealing runs of the mixed bilayers containing 24 SDBS molecules (Set G) were conducted, where the final configurations from the Set G runs were heated from 350 to 370 K using the same procedure as described above for pure TS bilayers. Fig. 10 shows the change in the area of the bilayers as a function of temperature, for two independent runs (derived from runs G1 and G3). The presence of the SDBS molecules has reduced the phase transition temperature of both bilayers below 365 K, with the melting temperature of run G1 being markedly lowered. Thus, while the incorporation of SDBS molecules did not affect the structural stability of TS bilayers at 350 K, it did appear to lower the phase transition temperature of these bilayers, relative to the pure TS bilayer case. The lateral size of the bilayer systems simulated, ∼23.3 nm², is such that we believe that the general trends regarding the phase behaviour of the system will not be significantly affected by system size effects. It is possible that system size effects could cause the melting temperature of the TS bilayer system to shift by a modest amount. However, the absolute melting temperatures are of less importance in this study than the fact that the presence of SDBS molecules causes a depression of the melting temperature. As the system size of the bilayers only varied modestly for Set B and Set G we believe that this trend will hold true regardless of system size.

Fig. 10 The area of the mixed SDBS–TS bilayers as a function of temperature/time. Dashed lines indicated the transition temperature for each run. The solid line indicates the average transition temperature calculated for the pure TS bilayer samples.

The interaction of SDBS molecules with TS bilayers noted from our simulations is consistent with the insights obtained from experiments, in that the SDBS molecules were able to spontaneously insert themselves into the TS structure.¹⁹ The incorporation of SDBS molecules into the TS bilayers does cause small changes in the bilayer structure, causing the membrane to expand laterally and thin in the transverse direction, as well making the acyl tails of the TS molecules more disordered. However, the general structural features of the mixed SDBS–TS bilayer were unchanged at 350 K. The most significant impact of SDBS insertion was the slight reduction of the gel–liquid crystalline phase transition temperature.

Conclusions

We have reported results of MD simulations of a triacylglycerol (TAG), namely tristearin (TS), in bilayer form in water. We investigated the structure and stability of the TS-aqueous interface at different temperatures and in the presence of an anionic surfactant, SDBS. At temperatures of 350 K and below, the TS bilayers were predicted to be structurally stable and in a gel phase, with the lipid tails laterally ordered in a hexagonal arrangement. In addition to the highly-ordered fatty acid tails, the bilayer showed a narrow interface width, indicating the relatively weak hydrophilic character of the TS headgroup. Between 365 and 370 K the TS bilayer melted; however, unlike other lipid bilayers, the TS bilayer could not sustain a stable liquid–crystalline phase bilayer. In simulations of TS bilayers in the presence of SDBS molecules, the SDBS molecules were found to spontaneously insert into the TS bilayers. The effects on the bilayer structure, imparted by different ratios of embedded SDBS molecules to TS chains in mixed bilayers, were investigated. A TS : SDBS ratio of 9 : 1 had no significant effect on the structure of the bilayer. Doubling the amount of SDBS incorporated in the bilayer did cause the bilayer to expand laterally, by ∼8%, while at the same time leading to a slight thinning in the transverse direction. Despite this, even bilayers with a TS : SDBS ratio of 3 : 1 showed only small changes in their overall structure at 350 K, compared with the pure TS case. The most significant impact of the incorporation of SDBS within the TS bilayer was the depression of gel–liquid crystalline phase transition temperature to below 365 K. Our study represents the first step in applying atomistic molecular dynamics simulations to the investigation of TAG-aqueous interfaces. Our results suggest that the CHARMM36 FF is suitable for the simulation of such systems, although the phase behaviour of the system may be shifted to lower temperatures than is the case for the actual system. Our findings provide a foundation for the further study of the effects of the presence of nanoparticles on/at the TS-aqueous interface, as well as a wider investigation of the effects of other surfactant molecules, consideration of different TAGs, and application to the study of inhomogeneous TAG bilayers.

Acknowledgements

This research was undertaken with the assistance of resources from the National Computational Infrastructure (NCI), which is supported by the Australian Government. ZEH and TRW thank veski for funding and an Innovation Fellowship for TRW.

References

1. K. Sato, Chem. Eng. Sci., 2001, 56, 2255–2265.
2. A. G. Marangoni, N. Acevedo, F. Maleky, E. Co, F. Peyronel, G. Mazzanti, B. Quinn and D. Pink, Soft Matter, 2012, 8, 1275–1300.
3. A. Almeida and E. Souto, Adv. Drug Delivery Rev., 2007, 59, 478–490.
4. S. A. Wissing, O. Kayser and R. H. Müller, Adv. Drug Delivery Rev., 2004, 56, 1257–1272.
5. W. Li and F. C. Szoka, Pharm. Res., 2007, 24, 438–449.
6. T. Lobovkina, G. B. Jacobson, E. Gonzalez-Gonzalez, R. P. Hickerson, D. Leake, R. L. Kaspar, C. H. Contag and R. N. Zare, ACS Nano, 2011, 5, 9977–9983.
7. M. D. Krebs, O. Jeon and E. Alsberg, J. Am. Chem. Soc., 2009, 131, 9204–9206.
8. P. Guo, O. Coban, N. M. Snead, J. Trebley, S. Hoeprich, S. Guo and Y. Shu, Adv. Drug Delivery Rev., 2010, 62, 650–666.
9. D. A. Pink, B. Quinn, F. Peyronel and A. G. Marangoni, J. Appl. Phys., 2013, 114, 234901.
10. F. Peyronel, J. Ilavsky, G. Mazzanti, A. G. Marangoni and D. A. Pink, J. Appl. Phys., 2013, 114, 234902.
11. A. Hall, J. Repakova and I. Vattulainen, J. Phys. Chem. B, 2008, 112, 13772–13782.
12. A. Brasiello, S. Crescitelli and G. Milano, Faraday Discuss., 2012, 158, 479–492.
13. A. N. Zdravkova and J. P. J. M. van der Eerden, J. Cryst. Growth, 2006, 293, 528–540.
14. A. N. Zdravkova and J. P. J. M. van der Eerden, J. Cryst. Growth, 2007, 307, 192–202.
15. A. N. Zdravkova and J. P. J. M. van der Eerden, Cryst. Growth Des., 2007, 7, 2778–2787.
16. C. A. Miller and K. H. Raney, Colloids Surf., A, 1993, 74, 169–215.
17. A. W. Sonesson, T. H. Callisen, U. M. Elofsson and H. Brismar, J. Surfactants Deterg., 2007, 10, 211–218.
18. K. Thirunavukarasu, N. G. Edwinoliver, S. Durai Anbarasan, M. K. Gowthaman, H. Iefuji and N. R. Kamini, Process Biochem., 2008, 43, 701–706.
19. X. Cui, X. Liu, A. S. Tatton, S. P. Brown, H. Ye and A. Marsh, ACS Appl. Mater. Interfaces, 2012, 4, 3225–3232.
20. C. Tongcumpou, E. J. Acosta, L. B. Quencer, A. F. Joseph, J. F. Scamehorn, D. A. Sabatini, S. Chavadej and N. Yanumet, J. Surfactants Deterg., 2003, 6, 205–214.
21. P. Tanthakit, A. Nakrachata-Amorn, J. F. Scamehorn, D. A. Sabatini, C. Tongcumpou and S. Chavadej, J. Surfactants Deterg., 2009, 12, 173–183.
22. R. W. Corkery, D. Rousseau, P. Smith, D. A. Pink and C. B. Hanna, Langmuir, 2007, 23, 7241–7246.
23. D. A. Pink, C. B. Hanna, C. Sandt, A. J. MacDonald, R. MacEachern, R. Corkery and D. Rousseau, J. Chem. Phys., 2010, 132, 054502.
24. A. K. Sum, M. J. Biddy, J. J. de Pablo and M. J. Tupy, J. Phys. Chem. B, 2003, 107, 14443–14451.
25. W.-D. Hsu and A. Violi, J. Phys. Chem. B, 2009, 113, 887–893.
26. M. Greiner, A. M. Reilly and H. Briesen, J. Agric. Food Chem., 2012, 60, 5243–5249.
27. M. Greiner, B. Sonnleitner, M. Mailänder and H. Briesen, Food Funct., 2014, 5, 235–242.
28. H. Dominguez and M. Berkowitz, J. Phys. Chem. B, 2000, 104, 5302–5308.
29. E. Clavero, J. Rodriguez and D. Laria, J. Chem. Phys., 2007, 127, 124704.
30. M. Sammalkorpi, M. Karttunen and M. Haataja, J. Phys. Chem. B, 2007, 111, 11722–11733.
31. M. Sammalkorpi, M. Karttunen and M. Haataja, J. Am. Chem. Soc., 2008, 130, 17977–17980.
32. M. Sammalkorpi, M. Karttunen and M. Haataja, J. Phys. Chem. B, 2009, 113, 5863–5870.
33. F. Palazzesi, M. Calvaresi and F. Zerbetto, Soft Matter, 2011, 7, 9148–9156.
34. S. Storm, S. Jakobtorweihen, I. Smirnova and A. Z. Panagiotopoulos, Langmuir, 2013, 29, 11582–11592.
35. M. Tarek, D. J. Tobias and M. L. Klein, J. Phys. Chem., 1995, 99, 1393–1402.
36. S. Yuan, L. Ma, X. Zhang and L. Zheng, Colloids Surf., A, 2006, 289, 1–9.
37. G. F. Catá, H. C. Rojas, A. P. Gramatges, C. M. Zicovich-Wilson, L. J. Álvarez and C. Searle, Soft Matter, 2011, 7, 8508–8515.
38. L. Ferrer-Tasies, E. Moreno-Calvo, M. Cano-Sarabia, M. Aguilella-Arzo, A. Angelova, S. Lesieur, S. Ricart, J. Faraudo, N. Ventosa and J. Veciana, Langmuir, 2013, 29, 6519–6528.
39. X. He, O. Guvench, A. D. MacKerell Jr and M. L. Klein, J. Phys. Chem. B, 2010, 114, 9787–9794.
40. M. Suttipong, N. R. Tummala, B. Kitiyanan and A. Striolo, J. Phys. Chem. C, 2011, 115, 17286–17296.
41. J. W. Gershman, J. Phys. Chem., 1957, 61, 581–584.
42. H. M. Vale and T. F. McKenna, Colloids Surf., A, 2005, 268, 68–72.
43. J. B. Klauda, R. M. Venable, J. A. Freites, J. W. O'Connor, D. J. Tobias, C. Mondragon-Ramirez, I. Vorobyov, A. D. MacKerell Jr and R. W. Pastor, J. Phys. Chem. B, 2010, 114, 7830–7843.
44. R. W. Pastor and A. D. MacKerell Jr, J. Phys. Chem. Lett., 2011, 2, 1526–1532.
45. W. L. Jorgensen, J. Chandrasekhar, J. D. Madura, R. W. Impey and M. L. Klein, J. Chem. Phys., 1983, 79, 926–935.
46. E. Neria, S. Fischer and M. Karplus, J. Phys. Chem., 1996, 105, 1902–1921.
47. B. Hess, C. Kutzner, D. Van Der Spoel and E. Lindahl, J. Chem. Theory Comput., 2008, 4, 435–447.
48. S. Nosé, J. Chem. Phys., 1984, 81, 511–519.
49. W. Hoover, Phys. Rev. A, 1985, 31, 1695–1697.
50. M. Parrinello and A. Rahman, J. Appl. Phys., 1981, 52, 7182–7190.
51. R. Notman, W. K. den Otter, M. G. Noro, W. J. Briels and J. Anwar, Biophys. J., 2007, 93, 2056–2068.
52. C. Das, P. D. Olmsted and M. G. Noro, Soft Matter, 2009, 5, 4549–4555.
53. C. Das, M. G. Noro and P. D. Olmsted, Biophys. J., 2009, 97, 1941–1951.
54. M. I. Hoopes, M. G. Noro, M. L. Longo and R. Faller, J. Phys. Chem. B, 2011, 115, 3164–3171.
55. R. Notman and J. Anwar, Adv. Drug Delivery Rev., 2013, 65, 237–250.
56. S. Guo, T. C. Moore, C. R. Iacovella, L. A. Strickland and C. McCabe, J. Chem. Theory Comput., 2013, 9, 5116–5126.
57. S. Leekumjorn and A. K. Sum, Biochim. Biophys. Acta, Biomembr., 2007, 1768, 354–365.
58. J. M. Rodgers, J. Sørensen, F. J. M. de Meyer, B. Schiøtt and B. Smit, J. Phys. Chem. B, 2012, 116, 1551–1569.
59. J. P. Jämbeck and A. P. Lyubartsev, J. Chem. Theory Comput., 2012, 8, 2938–2948.
60. T. J. Piggot, Á. Piñeiro and S. Khalid, J. Chem. Theory Comput., 2012, 8, 4593–4609.
61. D. Poger and A. E. Mark, J. Chem. Theory Comput., 2012, 8, 4807–4817.
62. S. J. Marrink, J. Risselada and A. E. Mark, Chem. Phys. Lipids, 2005, 135, 223–244.

---

Footnote

† Electronic supplementary information (ESI) available: The area of example bilayers from the Set A and B simulations; the 2D radial distribution function of the TS acyl tails; the order parameter of the acyl tails of the TS molecules in simulations of Set A and B; snapshots showing the gel–liquid crystalline phase transition; snapshots showing the presence of SDBS aggregates; the number of SDBS molecules embedded in the bilayer during the Set C, D, E, F and G simulation runs; density profiles of the Set C, D, E, F and G simulation runs; snapshots showing the packing of the SDBS molecules in the bilayers; order parameters of the acyl tails of TS and SDBS molecules; the sulphur/phenyl ring–water radial distribution functions of SDBS. See DOI: 10.1039/c5ra09192f

---