Composition effect on the aggregate/solution interface of a nematic lyotropic liquid crystal

Abstract

Lyotropic nematic liquid crystals (lnlc) have attracted attention due to their resemblance to natural membranes and several technological applications. The effect of tetradecyltrimethyl ammonium chloride (TTAC) in the viscosity of a new lnlc containing also decanol, NaCl and natural lipids, was studied. In this article we study how the concentration of phospholipids, decanol and NaCl affects the bilayer and the liquid crystal/solution interface structures. ²H-NMR quadrupole splittings and longitudinal relaxation times of several deuterated species were measured. Polarized light microscopy, transmission electron microscopy and kinematic viscosity measurements were also obtained. The structure of the bilayer is dominated by the surface charge density of TTAC. Water rotation, translation and ion solvation depend on the charge and distance from the interface. With increasing TTAC, MD predicts hydrocarbon chain inter-digitation and lipid rafts formation. This was rationalized using energy arguments.

---

1 Introduction

Different lyotropic liquid crystalline phases result from the self-assembly of amphiphilic molecules, such as lipids and surfactants in water when their concentration achieves a certain value because of hydrophobic and electrostatic interactions.^1–4 In this context, water plays a crucial role in regulating the forces that drive the formation of these molecular aggregates. Different properties such as ²H-NMR quadrupole splitting (Δv_q), the longitudinal relaxation time (T₁), water re-orientation, the translational diffusion coefficient and ion solvation, among others, can be modified by changes in the composition of the aggregate/solvent interface. For instance, different types of water molecules, in terms of their physico-chemical properties, are expected in the solution, depending on their distance from the interface, and it is very likely that the structure of the bilayer itself is also affected by its composition.

A new nematic lyotropic liquid crystal formed by tetradecyltrimethyl ammonium chloride (TTAC), decanol (DeOH), sodium dodecyl sulphate (SDS), sodium chloride (NaCl) and a natural mixture of lipids extracted from soybean (PL), constituted mostly by phosphatidylcholine (DOPC) and phosphatidylethanolamine (POPE), all dissolved in water, has been previously prepared. This nematic liquid crystal shows that variations of about 30% w/w in TTAC content continuously increase the magnitude of the viscosity by more than 24 times, making them suitable for certain biomedical applications as lubricants for artificial implants in hips, shoulders and knees, and certain therapies against osteo-arthritis or rheumatoid-arthritis diseases.⁵ Therefore, since the effects of TTAC content on the structure and dynamics of the aggregate itself were already studied,⁵ in this work we explore the consequences of introducing variations in DeOH, PL and NaCl content on the structures of the bilayer and the interface. Δv_q and T₁ from DHO and Δv_q from DeOH-α-d₂, were measured to explore the dynamics near the interface, mainly to locate the phase transitions. Δv_q and T₁ from DHO and SDS-d₂₅ were measured in selected samples, to probe the dynamics at the interface and towards the interior of the hydrophobic core. Transmission electron microscopy (TEM) and polarized light microscopy (PLM) images of representative compositions were also obtained. Because the structure of the molecular aggregates that form the liquid crystal appears to be very dependent on the surface charge density provided by TTA⁺, here we study how it affects solvent properties going from the interface to the bulk aqueous solution. Effectively, water dynamic properties such as the diffusion coefficient, rotational correlation time and ion solvation, were obtained from previously simulated trajectories. To validate the molecular dynamics simulations, T₁ values of DHO were measured and the results compared with values obtained from simulation. Finally, the effect of variations in TTAC composition on the structures of the bilayer and interface were examined as well.

2 Materials and methods

2.1 Mesophase preparation

TTAC, SDS, SDS-d₂₅, D₂O, DeOH and the PL mixture were purchased from Aldrich. Water of HPLC grade and NaCl at 99.8% purity were purchased from Merck. All reagents were used as received.

DeOH-α-d₂ was synthesized by reduction of ethyl decanoate ester with LiAlD₄ and purified by vacuum distillation. 5% w/w D₂O was added to the solvent to provide ²H-NMR signal from DHO. For a detailed characterization of the dynamics towards the interior of the hydrophobic core, SDS was replaced by SDS-d₂₅ in selected samples.

In order to explore concentration effects on the interface properties, from components not studied before, we have prepared three series of samples, one for DeOH, one for PL and a third series extending the range of the NaCl concentration. Consequently, a series of 4 samples containing 0.2253 mg TTAC, 0.0912 mg NaCl, 0.0481 mg DeOH, 0.012 mg SDS and 0.0466, 0.0816, 0.1165 and 0.1400 mg PL, were dissolved in 1 ml of water each. Similarly, a second series of 4 samples containing 0.1165 mg PL, 0.2253 mg TTAC, 0.0912 mg NaCl, 0.012 mg SDS and 0.0351, 0.0442, 0.0506 and 0.0547 mg of DeOH was prepared. Finally, a third series of 4 samples containing 0.1165 mg PL, 0.2253 mg TTAC, 0.0481 mg DeOH, 0.012 mg SDS and 0.0623, 0.0771, 0.0912 and 0.1070 mg NaCl each, was also prepared. All mentioned samples contained DeOH-α-d₂ and D₂O. Selected samples, representative of the different phases found, were reproduced by using DeOH and SDS-d₂₅ instead of DeOH-α-d₂ and SDS.

Before any measurement was performed, all mesophases were allowed to equilibrate at least 48 hours at either 300 K or 298 K, for NMR or kinematic viscosity (η), respectively.

2.2 ²H-NMR spectra

All NMR spectra were registered using a Bruker Avance 400 NMR spectrometer, located at the Universidad de Santiago de Chile, employing a broadband probe tuned to 61.425 MHz. The ²H 90° pulse was 19 μs long and more than 1000 transients from a spectral window of 40 kHz were accumulated in 32 kB files. T₁ were obtained using the standard T1IR experiment.

2.3 Polarized light microscopy (PLM) and kinematic viscosity (η)

A Motic series B microscope equipped with crossed polarizers was used to obtain the textures. To orient the aggregates, the samples were mounted in a concave glass microscope and exposed to a 1.4 T magnet for 15 minutes; more time in the magnet did not modify the textures. Photographs were obtained at 300 K, along the direction of the field, with a 20-fold magnification, from the center of the concavity, to assure that all textures arise from the same sample thickness.

Kinematic viscosities of the solutions, η, were measured using a standard Ubbelohde viscometer at 298 K and 1 atmosphere.

2.4 Transmission electron microscopy

TEM images were collected with a JEOL JEM-1400PLUS instrument operating at 120 kV, using carbon coated 400 square mesh copper grids. Samples were drop casted onto the TEM grids. All TEM images were obtained in conditions identical to those previously reported.⁵

2.5 Molecular dynamics calculations

The GROMACS package 4.5.3 was used to calculate the molecular dynamics trajectories, and most of the analyses of this work were carried out using our own code. The integration time step was 4 fs. van der Waals forces were simulated using the Lennard-Jones potential (LJ), and electrostatic interactions were calculated with the particle mesh Ewald (PME) method.^6,7 In both cases, a cut-off of 1.2 nm was employed. The bond lengths were restrained using LINCS.⁸ All simulations were carried out in an NPT ensemble at 300 K and 1 atm, coupled to weak temperature and pressure bath algorithms,⁹ with time constants of 0.1 ps and 1 ps for temperature and pressure, respectively. The Ryckaert–Bellemans potential¹⁰ was used in all torsions along the aliphatic chains, for a better reproduction of cis–trans transitions. Simple point charge (SPC)¹¹ was the water model considered in all simulations. Atom numeration and charge distribution, can be found elsewhere.⁵

Finally, since most properties of the aggregate depend on the surface charge from TTA⁺, a series of calculations employing previously reported trajectories of different systems were performed. The calculated trajectories correspond to experimental samples containing 0.2253 g ml⁻¹ TTAC (system 1), 0.295 g ml⁻¹ TTAC (system 2) and 0.2253 g ml⁻¹ TTAC including 1 M NaCl (system 3).

3 Results and discussion

3.1 Interface composition

In order to modify the composition of the interface, PL, DeOH and NaCl contents were modified and the splittings from DeOH-α-d₂ and DHO measured. The results of these experiments are displayed in Fig. 1 of the ESI.† These experiments were performed mainly to locate the compositions at which the characteristic spectrum of the nematic anisotropic liquid crystal is lost, an unequivocal sign that the mesophase was destroyed. Once these limits had been identified, samples representative of the different nematic mesophases found were prepared incorporating SDS-d₂₅, to probe the order and dynamics not only at the interface but also towards the interior of the hydrophobic core.¹² Fig. 1 shows the ²H-NMR spectra from SDS-d₂₅ and DHO, Fig. 2 displays the PLM textures and Fig. 3 corresponds to TEM images arising from selected samples. Fig. 1A, 2A and 3A were reproduced from a previous work.⁵ Measured values of Δv_q and T₁ from SDS-d₂₅ and DHO are presented in Table 1, and values of η for these samples are listed in Table 2.

Fig. 1 ²H-NMR spectra from SDS-d₂₅ and DHO of selected samples. *Previously reported results.⁵

Fig. 2 PLM textures of selected samples. *Previously published results.⁵

Fig. 3 TEM images of selected samples. *Previously published results.⁵

Table 1 Experimental values of ²H-NMR quadrupole splittings (Δv_q, Hz) and longitudinal relaxation time (T₁, ms) of SDS-d₂₅ and DHO from selected samples. *Previously reported results⁵

Carbon	*0.225 g ml^−1 TTAC		*0.293 g ml^−1 TTAC		0.047 g ml^−1 PL		0.140 g ml^−1 PL
	Δvq/Hz	T1/ms	Δvq/Hz	T1/ms	Δvq/Hz	T1/ms	Δvq/Hz	T1/ms
1	25 371	34.2	20 153	45.9	22 966	40.4	25 792	35.4
2	25 371	34.2	20 153	45.9	22 966	40.4	25 792	35.4
3	25 371	34.2	20 153	45.9	22 966	40.4	25 792	35.4
4	25 371	34.2	20 153	45.9	22 966	40.4	25 792	35.4
5	25 371	34.2	20 153	45.9	22 966	40.4	25 792	35.4
6	25 371	34.2	20 153	45.9	22 966	40.4	25 792	35.4
7	22 734	36.7	20 153	45.9	22 966	40.4	22 812	38.7
8	20 750	50.2	17 386	45.9	19 432	66.9	20 988	43.3
9	18 016	61.5	15 054	69.5	16 907	76.3	18 101	56.2
10	14 815	81.9	12 371	186.8	13 890	94.7	14 919	71.7
11	10 852	112.4	9046	120.4	10 106	125.6	10 893	103.4
12	2889	308.3	2383	309.9	2635	352.1	28 885	289.4
DHO	18.8	338.9	28.7	317.7	5.8	365.5	26.2	331.4

---

Carbon	0.040 g ml^−1 DeOH		0.050 g ml^−1 DeOH		0.073 g ml^−1 NaCl		*0.104 g ml^−1 NaCl
	Δvq/Hz	T1/ms	Δvq/Hz	T1/ms	Δvq/Hz	T1/ms	Δvq/Hz	T1/ms
1	23 142	37.9	25 269	33.4	24 983	34.6	25 232	31.1
2	23 142	37.9	25 269	33.4	24 983	34.6	25 232	31.1
3	23 142	37.9	25 269	33.4	25 232	34.6	25 232	31.1
4	23 142	37.9	25 269	33.4	24 983	34.6	25 232	31.1
5	23 142	37.9	25 269	33.4	24 983	34.6	25 232	31.1
6	23 142	37.9	25 269	33.4	24 983	34.6	24 489	31.1
7	21 182	41.2	22 482	33.9	22 468	33.5	22 564	36.7
8	20 988	54.6	20 714	41.9	20 562	50.3	20 663	42.2
9	18 101	59.0	17 980	56.1	17 831	58.5	17 951	52.1
10	14 919	75.3	14 779	69.2	14 689	80.3	14 710	70.0
11	10 893	111.9	10 809	109.1	10 723	110.3	10 772	99.1
12	2888	289.5	2883	313.2	2865	307.4	2854	294.4
DHO	23.5	337.0	24.5	339.5	28.1	329.8	18.9	334.8

---

Table 2 Kinematic viscosity of selected samples. *Previously reported results⁵

Sample	Viscosity/mPa s
*0.225 g ml^−1 TTAC	12.9 ± 0.1
*0.293 g ml^−1 TTAC	310 ± 4
0.047 g ml^−1 PL	106.7 ± 0.8
0.140 g ml^−1 PL	8.0 ± 0.2
0.040 g ml^−1 DeOH	59.8 ± 0.4
0.050 g ml^−1 DeOH	7.7 ± 0.1
0.073 g ml^−1 NaCl	29.5 ± 0.2
*0.104 g ml^−1 NaCl	10.3 ± 0.1

---

An inspection of Table 1 reveals that, as expected, Δv_q is greater for the first CD₂ of the chain and decreases towards the interior of the bilayer. It is also observed that increasing either the content of PL or DeOH has the same effect on the Δv_q and T₁ of SDS-d₂₅; the first increases and the second decreases for all deuterium atoms towards the interior of the aggregate. This indicates that the degree of alignment of the aliphatic chains with the magnetic field increases and, consequently, their re-orientational dynamics becomes slower. This trend is opposite to that observed when increasing TTAC,⁵ and is very likely to arise from charge screening and spacing effects experienced by the TTA⁺ head-groups present at the interface when increasing the concentration of non-charged amphiphiles.

On the other hand, from Table 2 it is observed that augmenting the content of either PL or DeOH decreases η by about an order of magnitude in both cases. This is presumably because the aggregates become smaller. In effect, an increase of either DeOH or PL appears to be equivalent to a decrease in TTAC content, which in accordance with previous results,⁵ decreases the size of the aggregate. This makes sense if we observe that Δv_q from DHO increases in both cases. Decreasing the aggregate’s size generates a greater interfacial volume, able to accommodate more oriented solvent molecules, increasing the observed DHO splitting. An inspection of Table 1 reveals that the increase in DHO splitting is 20.4 Hz for PL and only 1 Hz for DeOH. Calculated radial distribution functions for the solvation of all head-groups (see Section 3.3) show that solvation capability of PL is greater than for DeOH. Effectively, the ratio between the first solvation layer of (DOPC + POPE)/DeOH is 2.1, supporting this interpretation.

Increasing the ionic strength does not affect significantly the order and dynamics of the amphiphiles.⁵ Despite the above, η decreased to almost one third of the original value with the increase in NaCl, presumably because of a size decrease. A plausible explanation arises from the forces involved in the aggregation process. The higher screening of charges at the interface, at higher ionic strength, allows smaller aggregates to form, since less aliphatic chain association is necessary to balance the electrostatic repulsions.

As previously reported,⁵ a TTAC increment of 30% w/w leads to significant differences in the SDS-d₂₅ ²H-NMR spectrum and solution η. These observations were explained because higher surface charge density from TTA⁺ increased repulsions among the head-groups, incrementing the mobility, decreasing the alignment and modifying the values of Δv_q and T₁. However, to explain the unusual η increment, another phenomenon, related to modifications in the size and structure of the aggregate itself, going from circular mono-axial to larger elliptic bi-axial, was invoked.

PLM images from selected samples (Fig. 2) shows the appearance of the same images observed when modifying TTAC content,⁵ oily streaks and schlieren textures, an indication that the previously observed phase transitions are also occurring with variations in the content of the other components of the aggregate. Previously we concluded that this corresponded to structural changes preserving the bilayer arrangement and nematic character, as evidenced by the ²H-NMR spectra. The oily streak texture is characteristic of homeotropic order, more likely to be found in smaller and better aligned aggregates, but is more usually observed in lamellar phases. Schlieren textures are usually found in nematic mesophases.¹³ Despite the similarities observed, there is an obvious difference appreciated in Fig. 2: schlieren textures appear at low concentrations of either PL or DeOH, however at higher concentrations oily streaks appear. This order is reversed in the case of TTAC. This observation reveals that it is effectively the surface charge density, i.e. the content of TTA⁺, which controls the morphology of the aggregates, since increasing the concentration of either PL or DeOH appears to be equivalent to decreasing the concentration of TTA⁺, inducing the same phase transitions observed before. To provide more experimental evidence of these changes, TEM images of the selected compositions were obtained and they are displayed in Fig. 3. It is possible to observe similar images to those previously obtained when modifying TTAC content, but in reverse order, corroborating the previous hypothesis.

Since many properties of the aggregate are mostly dependent on the surface charge of TTAC, a detailed theoretical study using previously reported MD trajectories⁵ is presented below. As mentioned before, in these simulations a variation in TTAC concentration and ionic strength were introduced. Therefore, the following study is oriented to explore surface charge effects on the bilayer itself, and to reveal how the charge and distance from the surface affect some dynamic properties of the solvent, such as the translational diffusion coefficient, rotational correlation time and ion hydration.

3.2 Rotational dynamics of water

To obtain information on the rotational dynamics of the water electric dipole moment near the liquid crystal interface, rotational correlation times at different distances from the center of the bilayer were estimated, defining parallel planes of equal thickness at different distances from the interface. Once the thickness of this planes was defined, rotational correlation times were calculated on the basis that each solvent molecule can be assigned to different slabs. To do this, the 200 ns trajectory length was split into sub-trajectories of 20 ps each and the correlation time calculated in each slice. Therefore, from a single trajectory the values going from the interface to the aqueous solution were calculated for systems 1, 2 and 3.

In general, water dipole moment re-orientational dynamics can be described by the rotational correlation function which can be formulated in terms of the first Legendre polynomial:

where correspond to the water dipole moment orientation separated by a time interval t.

〈P₁(t)〉 can be fitted to a multi-exponential function, as follows:

here a_i and τ_i correspond to different amplitudes and relaxation times associated with water dipole reorientation.

Since a_i and τ_i are not single valued, a τ_eff, independent of the fit, must be defined. Thus, τ_eff is defined as follows:¹⁴

Fig. 4 shows the best fit of 〈P₁(t)〉 as a function of time for different distances from the center of the bilayer in systems 1 and 2. Thus, from the multi-exponential fit just described, a τ_eff at different distances from the interface can be calculated. Fig. 5 displays τ_eff of the water dipole moment for systems 1, 2 and 3 as a function of the distance from the center of the bilayer, and shows how it abruptly increases when water approaches the interface. This is due to the fact that water at the interface coordinates the head-group atoms of PL and TTAC, introducing mobility restrictions. From this Figure, we observe how τ_eff decreases by roughly one order of magnitude when going from the interface to the bulk aqueous solution. However, for the case corresponding to system 2, this increase is less noticeable than in system 1. This decrease in τ_eff with the increase in TTAC content is associated with the relative decrease in PL content, since water coordinates preferentially to PL. Hence, a relative decrease in PL content is reflected by in a decrease in τ_eff as well. System 3 shows a similar behavior to system 1.

Fig. 4 〈P₁(t)〉, as a function of time for (a) system 1 and (b) system 2 at different distances from the center of the bilayer aggregate. Points represent the values obtained from simulation, and solid lines the regression obtained from the multi-exponential fit.

Fig. 5 Relaxation time of water dipole moment for the three studied cases, as a function of the distance from the center of the bilayer, Z = 0.

To assure that the simulations provide a good representation of the solvent dynamics, T₁ measurements of DHO from samples 1, 2 and 3 were carried out and compared with values obtained from simulations. The Redfield theory relates T₁ to the rotational correlation time and the electric field gradient experienced by the nucleus as follows:¹⁵

where J(ω₀) is the spectral density. Eqn (4) can be evaluated using τ_eff, calculated above. The electric field gradient is axially symmetric along the O–D bond, (η = 0.00), the term in square brackets is the quadrupole coupling constant in angular frequency units, 308 kHz for the O–D bond, and ω₀ is the Larmor frequency. Furthermore, assuming that the rotational correlation function is a single exponential decay, the value of T₁ can be easily calculated.¹⁵ T₁ from DHO were measured and values of 339 ms, 318 ms and 334 ms were obtained for samples 1, 2 and 3 respectively. These values are in reasonable agreement with the values 290 ms, 265 ms and 270 ms estimated for bulk water in systems 1, 2 and 3, respectively. Therefore, the simulations appear to represent the solvent dynamics reasonably well.

3.3 Translational diffusion coefficient of water, D_t

Translational diffusion coefficient of water inside a slice in the xy-plane can be calculated from the mean square displacement as follows:where n is the number of degrees of freedom, 2 in this case, r² is the square displacement parallel to the xy-plane and t is the time. To provide insights into the variations of water translational dynamics when going from the interface plane to the aqueous bulk solution, D_t at parallel planes of equal thickness from the interface were calculated using a methodology similar to that employed in Section 3.2. These results are displayed in Fig. 6. System 3 did not show differences with system 1 and is not included in the Figure.

Fig. 6 Water translational diffusion coefficient, D_t (solid line), water density (dotted line) and TTAC nitrogen distribution (slash line), from bulk water to the interior of the bilayer. (S1) System 1, (S2) system 2.

Fig. 6 shows the calculated values of D_t, ranging from about 1.0 × 10⁻⁵ cm² s⁻¹ near the interfaces to 5.8 × 10⁻⁵ cm² s⁻¹, at the bulk of the aqueous phase. These values are in agreement with both experimental and simulation results. In bulk water, D_t = 7.5 × 10⁻⁵ cm² s⁻¹ has been reported from the SPC model,¹⁶ and values ranging from 2.3 × 10⁻⁵ cm² s⁻¹ to 2.74 × 10⁻⁵ cm² s⁻¹ have been measured using ¹H-NMR¹⁷ and ²H-NMR.¹⁸

At deeper zones of the bilayer interface, the water self-diffusion coefficient becomes one fifth the value at the bulk aqueous solution. This decrease is mainly associated with the fact that near the interface, water molecules more resemble solvation water than free water. Furthermore, Fig. 6 also shows that there is an important decrease in the width of the hydrophobic core when increasing TTAC, which makes the profiles displayed in Fig. 6 significantly different. As previously observed,⁵ this reduction has its origin in the inter-digitation of TTAC aliphatic chains when increasing TTAC content and is discussed in Section 3.5.

Besides, a study of water coordination around the atoms of PL, DeOH, and TTA⁺ head-groups was performed. The radial distribution function, g(r) is defined by:

here N(r) is the number of atoms in a spherical shell at distance r and thickness ∂r, and ρ is the macroscopic number density, taken as the ratio between the mass of all atoms and the volume of the simulation box. In effect, Fig. 7 shows some of these plots and is observed that water is mainly coordinating PL head-groups, and as a consequence, a decrease in the self diffusion coefficient of water in the vicinity of the interface is expected. The head groups of PL and DeOH in system 1 can hold 5.7 and 1.5 water molecules in the first solvation layer, respectively. These values, weighted by the mole fraction of each species, allowed us to estimate the relative degree of solvation. Therefore we can estimate the ratio between the solvation water of PL vs. DeOH to be 2.1. This explains the significant difference observed in the splitting of DHO when increasing either DeOH or PL.

Fig. 7 Radial distribution function of water around the phosphate oxygen of DOPC corresponding to systems 1, 2 and 3.

3.4 Sodium and chloride hydration shell

Another aspect related to the stability of the liquid crystal is the structure of ions and the dehydration process when approaching the bilayer/solution interface. In this regard, a study on the solvation of sodium and chloride ions at different distances from the center of the bilayer, similar to that of Sections 3.2 and 3.3, has been performed. From integration of the radial distribution function of water around them, the hydration of these ions can be estimated at the previously defined slabs.

Fig. 8 shows the results of this study for sodium. Systems 1 and 2 give the same hydration number in the bulk aqueous solution, 5.00 ± 0.002 and 4.90 ± 0.01 respectively. These values agree with the hydration numbers of sodium ions measured by X-ray and from simulations, in which sodium hydration numbers in a range from 4 to 6 water molecules have been reported.¹⁹ Fig. 8 also shows how in both cases, a dehydration of almost 50% takes place for the sodium ions, going from the bulk aqueous solution to the interface.

Fig. 8 Sodium hydration number from bulk water to deeper zones of the hydrophobic core for two bilayers, corresponding to systems 1 and 2.

Fig. 9 shows the variation of the chloride hydration number from the bulk aqueous solution to the center of the bilayer, for systems 1 and 2. In bulk solution, hydration numbers of 7.98 ± 0.04 and 8.06 ± 0.03 were obtained for systems 1 and 2, respectively. These values agree with the experimental data measured by X-ray spectroscopy²⁰ and neutron diffraction,^21–23 as well as that calculated from simulations,^19,24 for which coordination numbers from 5.2 to 7.4 were measured. In all cases studied in this work, the counterions completely released their hydration shell when moving from bulk water towards the center of the bilayer. Effectively, a variation from 7 water molecules in the hydration state in bulk water to almost naked ions was estimated for chloride when embedded in the middle of the hydrophobic core. Furthermore, Fig. 9 also reproduces the reported experimental behavior in which the coordination number decreases with the increase in ionic strength.²⁵ Thus, from simulation, we obtained a 9% decrease in the chloride hydration number when the salt concentration increased to 1 M, following the same behavior as that reported by Cummings et al.,²⁵ in which a decrease of 11% was measured for the chloride hydration number with the increase in ionic strength. This significant decrease in the hydration shell of chloride ions is compensated by chloride coordination to the positive charges of TTA⁺, as evidenced by the calculated radial distribution function of Cl⁻ around TTA⁺ (see Fig. 2 ESI†).

Fig. 9 Chloride hydration number from bulk water to the deepest zone of the hydrophobic core for systems 1 and 2.

3.5 Bilayer structure

In this section we explore in more detail effects of TTAC content on the structure of the bilayer itself. In this regard, Fig. 10 shows the density profiles of TTA⁺ across the Z axis, perpendicular to the bilayer surface. It is observed that system 1 shows slightly overlapping hydrocarbon chains, as deduced from the valley located at the center of the plot. This profile contrast with that of system 2, in which a maximum in TTA⁺ density is observed at the center of the bilayer. Furthermore, system 2 displays a narrower profile. These observations clearly indicate that inter-digitation of the aliphatic chains of TTA⁺ from both leaflets is occurring. This can be understood on the basis of electrostatic and Lennard-Jones interactions. Increasing TTAC content by 30% w/w produces an increment in the ammonium electrostatic repulsions, from 113 kJ mol⁻¹ to 312 kJ mol⁻¹. This is compensated with a better compaction of the hydrocarbon chains, increasing TTA⁺–TTA⁺ Lennard-Jones attractions from −4157 kJ mol⁻¹ to −6460 kJ mol⁻¹.

Fig. 10 TTAC density across the molecular aggregate for systems 1 and 2. The origin of the Z-axis was taken at the center of the bilayer.

Furthermore, a visual inspection of the simulated bilayer interfaces reveals the formation of lipid rafts by phospholipid segregation. Fig. 11 shows a snapshot of system 1 after 160 ns simulation and is possible to distinguish the existence of clusters formed by PL only. To provide a thermodynamic explanation for the formation of lipid rafts, we have calculated the potential of mean force (PMF) corresponding to the insertion of DOPC into two ideal systems. Thus, in the first, the insertion of DOPC into a bilayer with a composition identical to system 1, except that all TTAC were removed, was simulated. In the second, DOPC was inserted into a bilayer of TTAC only. The PMF at different depths of the bilayer can be used to estimated the spontaneity of the insertion process in both ideal systems as follows:

here ΔG(z) is the difference between the free energy from the bulk aqueous solution to a certain position z, along an axis perpendicular to the bilayer interface. R is the gas constant, T the absolute temperature, c_i(z) the concentration of DOPC at position z, and c⁰_i the concentration of DOPC in the bulk aqueous solution. Due to sampling incompleteness, c_i(z) profiles cannot be obtained straightforwardly from normal MD simulations, however, using Umbrella²⁶ and WHAM²⁷ methodologies, such as they have been described in previous works,^14,28–30 the PMF associated with the insertion of DOPC into the two systems proposed above was estimated. To improve the sampling, two DOPC molecules were introduced into the system, one placed at the bulk water, at a distance of −3.9 nm from the center of the bilayer, and the other placed at the center of the hydrophobic core. Thus, both DOPC molecules were simultaneously shifted 0.1 nm each maintaining a distance of 3.9 nm between them. The motion of the center of mass of each molecule was subjected to a restrain along the Z axis using a harmonic potential of 3000 N nm⁻² mol⁻¹, within windows of 0.1 nm thick, but maintaining the freedom of movement in the xy plane. For each of the 40 subsystems, 30 ns of trajectory length was calculated.

Fig. 11 Snapshot of the xy-plane of system 1 after 160 ns of trajectory. Darker areas correspond to clusters of lipids.

Fig. 12 shows the PMF associated with the DOPC insertion from bulk water to the interior of the two models described above. This Figure evidences that for any distance along the Z axis, DOPC prefers to be surrounded by other DOPC molecules rather than by TTAC. Therefore, according with this result, the formation of lipid rafts in the aggregate is expected, since DOPC shows greater affinity than other phospholipid molecules. Finally, the formation of lipid rafts is favored by about 5 kJ mol⁻¹ with respect to insertion in a TTA⁺ bilayer.

Fig. 12 Potential of mean force (PMF) associated with the insertion of DOPC from the bulk aqueous solution into the two molecular bilayers mentioned above.

4 Conclusions

The use of MD simulations in the interpretations of the experimental results arising from complex systems, such as lyotropic nematic liquid crystals, has proven to be very useful. In this particular case, the dynamics of the solvent appears to be adequately represented.

All evidence strongly suggests that the surface charge density, provided by TTA⁺, is responsible not only for the size and shape of the aggregates, but also has a strong influence in the structure of the bilayer itself, inducing inter-digitation and lipids rafts formation.

Dynamical properties of the solvent, such as translation diffusion coefficient and rotational correlation time, are strongly dependent on the surface charge density and the distance from the interface. This effects disappear at a distance of 3.5 nm away from the interface.

Finally, lipid rafts formation arises from a greater affinity between PL molecules rather than PL–TTA⁺.

Acknowledgements

The authors are pleased to acknowledge financial support from Fondecyt-Chile (Grant No. 1150138). A. R. R. F. acknowledges a fellowship from Conicyt. The authors also acknowledge the Staff of the Computing Center at the Polytechnic University of Cartagena for their support and facilities. J. J. L. C. and F. G. D. B. acknowledge the financial support from Fundación Séneca (Grants No. 19353/PI/14 and 19499/PI/14, respectively). J. J. G. C. acknowledges the Ministry of Economy and Competitiveness for a Ramon y Cajal contract (#RyC-2014-14956).

References

1. C. Tanford, The Hydrophobic Effect: Formation of Micelles and Biological Membranes, J. Wiley, 2nd edn, 1980.
2. J. Moreno-Razo, E. Sambriski, N. Abbott, J. Hernandez-Ortiz and J. de Pablo, Nature, 2012, 485, 86–89.
3. F. Biedermann, W. Nau and H. Schneider, Angew. Chem., Int. Ed., 2014, 53, 11158–11171.
4. X. Wang, D. Miller, E. Bukusoglu, J. de Pablo and N. Abbott, Nat. Mater., 2015, 15, 106–112.
5. A. Ruiz-Fernández, J. López-Cascales, J. Giner-Casares, R. Araya-Maturana, F. Díaz-Baños, D. Muñoz-Gacitúa and B. Weiss-López, RSC Adv., 2016, 6, 7455–7464.
6. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
7. U. Essmann, L. Perera, M. Berkowitz, T. Darden, H. Lee and L. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
8. B. Hess, H. Bekker, H. Berendsen and J. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
9. H. Berendsen, J. Postma, W. van Gunsteren, A. DiNola and J. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
10. J. Ryckaert and A. Bellemans, Faraday Discuss. Chem. Soc., 1978, 66, 95–106.
11. H. Berendsen, J. Postma, W. van Gunsteren and J. Hermans, Intermolecular Forces, D. Reidel Publishing Company, 1981.
12. M. F. Brown, J. Magn. Reson., 1979, 35, 203–215.
13. H. Stegemeyer and H. Behret, Liquid Crystals, Steinkopff, 1994.
14. V. E. Bahamonde-Padilla, J. J. López-Cascales, R. Araya-Maturana, M. Martínez-Cifuentes and B. E. Weiss López, ChemPhysChem, 2014, 15, 1422–1431.
15. A. Abragam, Principles of Nuclear Magnetism, Oxford University Press, 2011, p. 314.
16. J. López Cascales, J. García de la Torre, S. Marrink and H. Berendsen, J. Chem. Phys., 1996, 104, 2713–2720.
17. D. O’Reilly, J. Chem. Phys., 1974, 60, 1607–1618.
18. K. Tanaka, J. Chem. Soc., Faraday Trans. 1, 1975, 71, 1127–1131.
19. J. Chandrasekhar, D. Spellmeyer and W. Jorgensen, J. Am. Chem. Soc., 1984, 106, 903–910.
20. N. Skipper, S. Cumming, G. Neilson and J. Enderby, Nature, 1986, 321, 52–53.
21. M. Yamagami, H. Wakita and T. Yamaguchi, J. Chem. Phys., 1995, 103, 8174–8178.
22. D. Powell, A. Barnes, J. Enderby, G. Neilson and P. Salmon, Faraday Discuss. Chem. Soc., 1988, 85, 137–146.
23. D. Powell, G. Neilson and J. Enderby, J. Phys.: Condens. Matter, 1993, 5, 5723–5730.
24. J. López Cascales and T. Otero, J. Chem. Phys., 2004, 120, 1951–1957.
25. S. Cummings, J. Enderby, G. Neilson, J. Newsome, R. Howe, W. Howells and A. Soper, Nature, 1980, 287, 714–716.
26. G. Torrie and J. Valleau, J. Comput. Phys., 1977, 23, 187–199.
27. S. Kumar, D. Bouzida, R. Swensen, P. Kollman and J. Rosemberg, J. Comput. Chem., 1992, 13, 1011–1021.
28. J. L. MacCallum and D. P. Tieleman, J. Am. Chem. Soc., 2006, 128, 125–130.
29. J. Maccallum, W. Drew Bennett and D. Tieleman, J. Gen. Physiol., 2007, 129, 371–377.
30. J. Lopez Cascales, S. D. Oliveira Costa and R. Porasso, J. Chem. Phys., 2011, 135, 135103–135110.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: 10.1039/c6ra13597h

---