water system

Chemical Physics 483–484 (2017) 103–111

Contents lists available at ScienceDirect

Chemical Physics journal homepage: www.elsevier.com/locate/chemphys

Influence of water concentrations on the phase transformation of a model surfactant/co-surfactant/water system Raju Lunkad, Arpita Srivastava, Ananya Debnath ⇑ Department of Chemistry, Indian Institute of Technology Jodhpur, Jodhpur 342011, India

a r t i c l e

i n f o

Article history: Received 8 July 2016 In final form 23 November 2016 Available online 11 December 2016 Keywords: Water concentration Phase transition Micelle Lamellar phase

a b s t r a c t The influence of water concentrations on phase transformations of a surfactant/co-surfactant/water system is investigated by using all atom molecular dynamics simulations. At higher water concentrations, where surfactant (behenyl trimethyl ammonium chloride, BTMAC) to co-surfactant (stearyl alcohol, SA) ratio is fixed, BTMAC and SA self-assemble into spherical micelles, which transform into strongly interdigitated one dimensional rippled lamellar phases upon decreasing water concentrations. Fragmentation or fusions of spherical micelles of different sizes are evident from the radial distribution functions at different temperatures. However, at lower water concentrations rippled lamellar phase transforms into an LbI phase upon heating. Our simulations reveal that the concentrations of water can influence available space around the head groups which couple with critical thickness to accommodate the packing fraction required for respective phases. This directs towards obtaining a controlling factor to design desired phases important for industrial and medical applications in the future. Ó 2016 Elsevier B.V. All rights reserved.

1. Introduction Surfactants, co-surfactants, water systems self assemble into a wide variety of topologically distinct mesoscopic structures, among which micellar and lamellar structures have been studied extensively since long [1–9]. Micelles are aggregates in aqueous solution where hydrophilic heads of surfactants are in contact with water and hydrophobic tails form the core in the micelle center. Surfactants self-assemble into micelles at concentrations higher than their critical micellization concentration. The driving force of the formation of these macro-phases is an intricate balance between hydrophobic association at the surfactant-water interface and electrostatic repulsion of head groups each as a function of micelle size [10–12]. In a lamellar phase, surfactants adopt various interesting topologies to form the ripple (P b0 ), interdigitated gel phase (LbI ) due to a complex interplay between the head group interactions of the surfactant molecules and chain packing [13]. Topological features of micelles and lamellae comprising surfactants, polymers, peptides and proteins often mimic biological systems and are used as novel biomedical materials in pharmaceutical sciences, drug delivery and molecular imaging applications [14–17]. A better understanding of factors controlling the phase transformations of these macro-structures can help in developing ⇑ Corresponding author. E-mail address: [email protected] (A. Debnath). http://dx.doi.org/10.1016/j.chemphys.2016.11.014 0301-0104/Ó 2016 Elsevier B.V. All rights reserved.

the technological perspective related to cosmetic formulations, waste water treatment, detergents, secondary and tertiary oil recovery processes and more advanced nano-structure formation concepts [18–22]. Depending on the packing parameter defined by the ratio of the hydrophobic core volume to the area and the critical chain length of surfactants, surfactant-water systems can adopt the most favored structures ranging from micelles, rods, bilayers and vesicles [10,23–25]. Different experimental techniques such as X-ray diffraction, X-ray scattering, electron spin resonance, small-angle X-ray-scattering, electron microscopy, electron paramagnetic resonance, laser light scattering, nuclear magnetic resonance, differential scanning calorimetry and so on are employed to characterize different morphologies and respective physico-chemical properties of self-assembled micelles and bilayers [8,26–29]. Complementary to these experimental studies, computer simulations provide insights into the molecular level interactions responsible for transformations of one macro-structure into another. In recent years computer simulations have opened up a new direction on capturing the self assembly of surfactants into micelles or bilayers using multi-scale simulations. Several coarse grained (CG) models have been developed using all atom (AA) simulations or experimentally known thermodynamics behaviors to derive a complete phase diagram [30–36]. Although these models can rightly access the relevant time and length scale of micellization or lamellae formation, all atom (AA) molecular

104

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111

dynamics simulation is particularly useful to probe the atomic level picture. The kinetics of micelle formation via self-assembly is investigated by AA simulations although formation of finite size large aggregates is associated with it [37]. The self assembly of Triton-X micelles are studied at different concentrations and temperatures [38]. The effects of molecular parameters like the size, the chain length and the form of the head-group on cationic and anionic micelles are investigated [39]. Similarly, an asymmetric sawtooth 1D ripple phase consisting alternate non-interdigitated gellike domain, and a fully interdigitated (LbI ) domain, is reported [40,41] for the dipalmitoylphosphatidylcholine (DPPC) lipid bilayer using a united atom model (UA), where the ripple phase is found to be dependent on hydration and temperature. However, none of these studies shed light on how one micellar phase can be transformed to a lamellar phase by controlling a parameter such as hydration. In this paper, we report different phases obtained for aggregates made up of surfactant and co-surfactant molecules using all atom MD simulations. In contrast to the well studied lipid bilayer systems, reports on mixed surfactant micelles have not received much attention. Here we present simulations for aggregates made up of the cationic surfactants, behenyl trimethyl ammonium chloride (BTMAC) with co-surfactant, stearyl alcohol (SA) and illustrate that varying the water concentrations offers a flexible route to control the macro-structure. The MD simulations of the mixed surfactant aggregates reveal the presence of few micellar aggregates at a given water content. With decreasing water concentration at a fixed BTMAC to SA ratio, the micellar phase can be transformed to a one dimensional interdigitated ripple phase which disappears at higher temperature. The atomistic perspective emerging from the MD study indicates that decreasing the water concentration induces specific head-group area which is coupled with the critical chain length or thickness to obtain a packing parameter associated with a given phase. The study suggests a compositional route which can be exploited to control and design phases with desired structural functionalities. The organization of the paper is as follows, Section 2 describes the details of the simulation method, Section 3 describes the results and discussions for phase transformations with different water concentrations. In Section 4, we summarize our main results and conclude.

2. Simulation details Simulations are carried out for systems made up of the single tailed cationic surfactant BTMAC (chemical formula is CH3–(CH2)21–N+Cl-(CH3)3) with the co-surfactant SA (chemical formula is CH3-(CH2)17-OH). The force-field parameters for BTMAC and SA are used from our earlier study [42] with the SPC model for water [43]. Simulations are carried out by keeping BTMAC to SA ratio 2:1 at 283 K using gromacs-4.6.5 [44–46]. Although several simulations are carried out by varying water content from 84.56 wt% to 62.35 wt% (Fig. 1 in Supplementary Materials (SM) A), only two extreme concentrations respective to two different phases are reported since rests resulted either similar phases or co-existence of two phases. However, the critical water concentration below which a bilayer is formed is 67.28 wt%. The initial configurations of two systems with high (84.56 wt%) and low water (62.35 wt%) content are referred to as S1 and S2 (shown in Fig. 1(a) and (b)) respectively. These are created by randomly placing 100 BTMAC and 50 SA molecules in a cubic box and solvated with 16373 and 4950 water molecules for S1 and S2 respectively. 100 Cl ions are added to make the system neutral in both cases. The random initial configurations are energy minimized with the steepest descent algorithm [48] followed by a 100 ps NVT run. Next 20 ns NPT runs are carried out using Berend-

Fig. 1. Snapshots of BTMAC, SA, water and Cl ions for initial random configurations of (a) S1 at high water content, (b) S2 at low water content, final configurations of (a) A1 (replicated S1), (d) A2 (replicated S2). Color codes: Ice blue, BTMAC; blue, nitrogen (N+) head group for BTMAC; red, SA; yellow, oxygen (O) head group of SA; green, water; purple, (Cl) ions. All water molecules of A1 and A2 are not shown for clarity. The snapshots are created using vmd [47].

sen thermostat [49] with a coupling constant of 0.5 ps and 1 ps for S1 and S2 respectively. Berendsen semi-isotropic pressure coupling method [49] is used to maintain the pressure at 1 bar and trajectories are calculated in every 100 ps for both cases. A time step of 2 fs is used for all simulations. Van der Waals and Columbic interactions are cut off at 1.3 nm using Ewald method [50]. To overcome the system size effect, the final configuration of the 20 ns run for S1 is replicated 4 times along z direction and is referred to as A1 in Table 1. Similarly, the final configuration of S2 after 20 ns NPT runs is replicated twice each along x and y directions and is referred to as A2 in Table 1. Note, replicating in other directions does not change the results reported here. Next 300 ns NPT runs are carried out for A1 and A2 with the same set of parameters used for S1 and S2. The final configurations of A1 and A2 are shown in Fig. 1(c) and (d) respectively. Last 120 ns run of 300 ns NPT runs are analyzed for all results reported here. We have performed simulations with different box dimensions for both A1 and A2 to check the system size effect and reported in table S3 within SM. Radial distributions and density distributions of the systems show that the results reported in the paper do not differ with larger box lengths (see SM). We have investigated the effect of temperature on the micellar and lamellar phase by carrying out simulations in the temperature ranging from 283  345 K. We have used the final configuration of the small system (S1 and S2) obtained at 283 K as an initial configuration for the system at 300 K with the annealing rate of 29.4 ps/K. The system is equilibrated by a 100 ps NVT run followed by a 25–45 ns NPT run. Similarly, the simulations at remaining temperatures (323 K, 330 K, 338 K and 345 K) are carried out where the final configuration of the previous temperature is used as an initial configuration for the succeeding temperature with the same annealing rate for both S1 and S2. The equilibrated S1 and S2 at all temperatures are replicated in a similar way as at 283 K to obtain A1 and A2 at respective temperatures. 100 ns NPT runs are carried out for each of the systems with the same set of simulation parameters mentioned for the 283 K run.

105

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111

Table 1 Number of behenyl trimethyl ammonium chloride (BTMAC), stearyl alcohol (SA) and water molecules used in different compositions investigated. The ratio of BTMAC to SA are 2 : 1 in A1 and A2. The percentage of water decreases from A1 to A2. Note, bilayer normal is along Lx. Simulations with larger boxlengths are reported in SM. Bilayer A1 A2

Number of molecules BTMAC SA

Water

400 400

65492 19800

200 200

Water % 84.56 62.35

Boxlength (nm) Lx

Ly

Lz

NPT (ns)

11.45 17.26

11.45 17.26

18.08 3.26

The co-ordinates of surfactants are analyzed from simulations at different temperatures using following procedure to obtain the surfaces of head groups and Voronoi diagrams.

350

 Density distributions of head groups of BTMAC and SA are calculated from which a geometric criterion is used to select surfactants belonging to upper or lower leaflets (Lu or Ll ). Using the ðx; yÞ coordinates of the head-groups of surfactants, Delaunay triangulation is performed on sets of Lu and Ll separately. The surfaces zu ðx; yÞ and zl ðx; yÞ corresponding to the sets Lu and Ll respectively are obtained by interpolating the triangulated surface on the regular nx  ny grid. The mean positions of the upper and lower interface (zu0 and zl0 respectively) are defined as,

250

250

200

200

150

150

100

100

50

50

zu0 ¼

1 X ðiÞ z ; N i2L h

350

(a)

0

0

2

(b)

M1 M2 M3

300

ð1Þ

300 300

4

6

8

B1

300

0

10

0

2

4

8

6

10

u

and

zl0 ¼

1 X ðiÞ z : N i2L h

ð2Þ

Fig. 2. Number of aggregations present in (a) A1 and (b) A2. For A1, the number of molecules are multiplied by a factor of 3 for the sake of clarity. The number decrease from M1 to M2 to M3. For A2, the number of molecules are same for both aggregates (B1).

l

The height fluctuations of the upper and lower leaflets are defined as,

hu ðx; yÞ ¼ hzu ðx; yÞ  zu0 i;

ð3Þ

and

hl ðx; yÞ ¼ hzl ðx; yÞ  zl0 i:

ð4Þ

 To obtain the specific area per head group of surfactant molecules, a Voronoi diagram is created using MATLAB for each leaflet including the head-group locations for both BTMAC and SA. 3. Results and discussions

decreases from A1 to A2, the micelles transform into two bilayers (B1) each with same number of surfactants (Fig. 2(b)). To find the shape and size of three classes of aggregates, M1, M2, and M3 in A1, radial distribution function (rdf, gðrÞ) of the head beads of BTMAC and SA is calculated with respect to the center of mass (COM) of last tail beads for every individual aggregates with a bin-width of 0:02 nm. Fig. 3(a) shows the rdf for M1, M2 and M3 after averaging over similar class of aggregates where r is the distance between the COM and the head beads. Since all three rdfs are symmetric with respect to the location of the peaks, the shape of M1, M2 and M3 are depicted to be spherical micelles which is evident from the snapshots as well (Fig. 1(c)). The locations of the peaks in Fig. 3(a) can be chosen as radii of the respective spheres which clearly show that the size of the spheres for M3 is

3.1. Shape and size of aggregates 120

(a)

1400

M1 M3 M2

(b)

1200 -3

80

B1

B1

)

100

1000 800

60

g(r)

Fig. 1 shows that the random initial configurations of two systems with high and low water content self-assemble to micellar and lamellar structures respectively. The equilibration of the self-assembled micellar and lamellar structures in A1 and A2 is monitored by calculating the time evolution of solvent accessible surface area (SASA, shown in Fig. S1 of Supplementary Materials A). It has been found that, the values of SASA converge after 180 ns of the 300 ns production run and hence the last 120 ns runs are analyzed for all results reported here. A1, with higher water content, self-assembles into few numbers of micelles, whereas A2, with lower water content self-assembles into a lamellar structure with two bilayers. The number of molecules in each aggregate is calculated for A1 and A2 and shown in Fig. 2(a) and (b) respectively. The plot shows that A1 self-assembles into 10 micelles, which can be classified into three categories named as M1, M2 and M3 depending upon the number of surfactants present in each micelle. The number of BTMAC-SA molecules in each micelle decreases from M1 to M2 to M3 for A1. As the water content

600 40 400 20 0

200 0 0.5 1 1.5 2 2.5 3 3.5 4

0

0

5

10

15

20

Fig. 3. (a) Radial distribution function (gðrÞ) of the head beads with respect to the center of mass of the last tail beads of the surfactants for M1, M2 and M3. (b) Density distributions of A2 along the bilayer normal. Head group densities are multiplied by a factor of 7.5 for clarity.

106

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111

less than the size of M2 which is again less than that of M1. This is consistent with Fig. 2(a) where it is shown that M1 has the highest number of BTMAC and SA molecules, whereas M3 with the smallest size has least number of surfactants present. To understand the structure of A2, a density profile of BTMAC and SA is calculated with respect to the normal of the lamellar phase with a bin-width of 0.01 nm. The density profile shows that there are two aggregates of similar density referring to an equal number of surfactants in each aggregate which is consistent with Fig. 2(b). Density profiles of head beads in Fig. 3(b) clearly demonstrate the presence of two asymmetric and periodic bilayers in A2 where the lower layers have higher density than the upper layers with a strong interdigitation between two layers. To understand the molecular structure of the surfactants present in A1, gðrÞ of different beads of BTMAC and SA are calculated with respect to the COM of last tail beads. Fig. 4(a)–(c) show the gðrÞ for M1, M2 and M3 respectively where the head beads are located far from the center of the spherical micelles, the tail beads reside closest to the center and the middle beads are located in between. Therefore, the surfactant tails are likely to be stretched and not folded back. To confirm the configuration of the surfactants in A1, we have calculated angle distributions of two vectors consisting tail to middle bead and middle to head bead. Location of the most probable distribution shows that both BTMAC and SA chains are stretched (see Fig. 5 in Supplementary Materials A). To monitor the compactness of the aggregates, radius of gyration (Rg ) is calculated for A1 using the formula below,

ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi sP m r2 P i i; Rg ¼ mi

ð5Þ

where mi is the mass of atom i at a distance r i from the COM of the aggregate core. Thus, Rg measures the distribution of the mass of the surfactants constituting the aggregate with respect to its COM and estimates the size of the aggregated structure. Fig. 5(a)–(c) show Rg of M1, M2 and M3 after averaging over surfactants present in respective class of aggregates for the last 120 ns production runs. Rg is calculated using both Eq. (5) and Rh from respective gðrÞ. The sizes of the error bars obtained from Eq. (5) are comparable to the

symbol sizes for individual M1, M2 and M3. Thus, their error bars are not shown. Since the location of the most probable g(r) of M1, M2 and M3 in Fig. 3(a) are similar to their hydrodynamic radius (Rh ), as the micelles are spherical, one can get an estimate of Rg pffiffi using Rg ¼ ð3=5Þ  Rh . Thus the average values of Rg from gðrÞ calculation turn to be 1.16 nm, 1.55 nm and 1.94 nm for M1, M2 and M3 respectively which are consistent with the Rg values calculated using Eq. 5 and shown in Fig. 5. The error bars shown in Fig. 5(a)–(c) are from the Rh calculations. However, Rg averaged over M1, M2 and M3 (Fig. 5(d)) from both calculations have larger error bars due to the size differences in M1, M2 and M3. The average Rg of A1 is 1:62  0:29 nm (shown in Table 2) which is close to the average Rh of the micellar core ( 2 nm, obtained from the average location of gðrÞ for all aggregates). Since Rg for A1 is constant with respect to time, the micellar structures do not change their shape in course of simulation time depicting the stability of the structures. To estimate the deviation of A1 and A2 from a perfect spherical shape, eccentricity (e) is calculated using the formula,

e¼1

3.2. Molecular shape in aggregates To understand the shape of molecules constituting the aggregates, deuterium order parameter (SCD ) is calculated for BTMAC and SA in A1 and A2 using the expression below,

100 80

60

60

40

40

g(r)

g(r)

120

(a)

80

20 0

ð6Þ

where Imin is the minimum of moment of inertia along x; y and z and Iav g is the average of the moment of inertia along these three axes. The value of e is zero for a perfect sphere and thus any deviation from a sphere can be measured from the value of e. Fig. 5(e) shows e of A1 and A2 with respect to simulation time after averaging over all aggregates. The value of e for A1 is close to zero (0:11  0:04) due to its spherical nature as seen in Fig. 1(c). Since the shape of A2 deviates from a sphere (evident from the snapshot of A2 shown in Fig. 1(d)), the value of e for A2 (0:60  0:01) is higher than that of A1. Time evolution of e for both A1 and A2 is constant due to the stability of the respective shapes.

120 100

Imin ; I av g

0 0.5 1 1.5 2 2.5 3 3.5

(b)

Head Middle Tail

20 0

0 0.5 1 1.5 2 2.5 3 3.5

120 100

(c)

80

g(r)

60 40 20 0

0 0.5 1 1.5 2 2.5 3 3.5

Fig. 4. Radial distribution function (gðrÞ) of the head bead, the middle bead and the last tail bead with respect to the center of mass of the tail beads of the surfactants for (a) M1, (b) M2 and (c) M3.

107

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111

Fig. 5. Time evolution of radius of gyration (Rg ) of (a) M1, (b) M2, (c) M3, (d) A1, (e) eccentricity (e) for A1 and A2 after averaging over all aggregates. e for A1 is less than A2 due to the more compact and spherical nature of the micelles in A1. Error bars of Rg calculated from Eq. 5 are not shown for individual M1, M2 and M3 since their sizes are of the sizes of the symbols. Error bars of Rg calculated from Rh are shown.

Table 2 The table shows different parameters of A1 and A2 at different concentrations of water. As the concentration of water decreases from A1 to A2, eccentricity (e), order parameter (SCD ) and critical chain length (lc ) increase and the area per head group (a0 ) decreases. System

Wt % of water

e

SCD

a0 nm2

lc nm

A1 A2

84.56 62.35

0.11 ± 0.04 0.60 ± 0.01

0.05 ± 0.031 4.19 ± 0.004

1.16 ± 0.31 0.38 ± 0.06

2.62 ± 0.12 3.14 ± 0.09

SCD ¼

1 ð3Cos2 h  1Þ; 2

ð7Þ

where h is the angle between the principal axis of each aggregate and the molecular axis of BTMAC or SA. SCD is calculated for every segment of CH2 along the chain of BTMAC or SA after averaging over 0.5

0.5

(b)

(a) 0.4

0.3

0.3

3.3. Influence of temperature on size and topology A1 A2

S CD

S CD

0.4

0.2

0.2

0.1

0.1

0

0

0 20 5 10 15 Atom position along the chains

all BTMAC or SA in A1 and A2 and shown in Fig. 6(a) and (b). For both BTMAC and SA, SCD of A1 is close to zero (0:05  0:03) due to the isotropic orientational order of the spherical micelles. SCD of A2 (4:19  0:004) is higher than that of A1 due to a fully stretched order of BTMAC and SA along the interface normal.

0 20 5 10 15 Atom position along the chains

Fig. 6. Deuterium order parameter (SCD ) of (a) BTMAC (b) SA in A1 and A2 after averaging over all aggregates. SCD of A1 is close to zero for both BTMAC and SA due to the spherical nature of the micelles in A1.

Influence of temperature on the size of the micelles and the topology of the bilayers is investigated by carrying out simulations of A1 and A2 at different temperatures (283 K, 300 K, 323 K, 330 K, 338 K and 345 K) as mentioned in Section 2. Significant changes in aggregation numbers are observed for A1 from 283 K to 330 K without any phase transformation. Fig. 7(a)–(d) show the number of aggregates present in A1 at 283 K, 300 K, 323 K and 330 K respectively. The figure demonstrates the presence of 10 aggregates at 283 K which get fused to 8 aggregates at 300 K. These get fragmented into 11 aggregates at 323 K, which again get fused into 8 aggregates at 330 K without any further change till 345 K (data not shown). Calculation of gðrÞ of A1 at different temperatures show how the sizes of the aggregates change with temperature. Fig. 8(a)–(c)show gðrÞ of M1, M2 and M3 respectively to depict the changes of their sizes from 283 K to 330 K. Since the nature of gðrÞ is symmetric with respect to the peak, micelles did not change their spherical shapes with temperature. The peak of gðrÞ of

108

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111

150

150

(b) 300 K

M1 M2 M3

100

Number of molecules

Number of moleclues

(a) 283 K

50

0

100

50

0

0 1 2 3 4 5 6 7 8 9 10 11 12

0 1 2 3 4 5 6 7 8 9 10 11 12

Aggregation number

Aggregation number 150

150

(d) 330 K Number of molecules

Number of molecules

(c) 323 K 100

50

0

100

50

0

0 1 2 3 4 5 6 7 8 9 10 11 12

0 1 2 3 4 5 6 7 8 9 10 11 12

Aggregation number

Aggregation number

Fig. 7. Number of aggregations present in A1 at (a) 283 K, (b) 300 K, (c) 323 K and (d) 330 K. Different aggregates get fused together or fragmented with different temperatures.

tail

150

- head

(a)

100

g(r)COM

g(r)COM

tail

- head

200

50 0

0

1

2

3

4

5

6

70 60 50 40 30 20 10 0

(b)

0

1

r (nm)

283 K 300 K 323 K 330 K

2

3

4

5

6

r (nm)

g(r)COM

tail

- head

150

(c) 100 50 0

0 0.5 1 1.5 2 2.5 3 3.5 4

r (nm) Fig. 8. Radial distribution function (gðrÞ) of the head beads with respect to the center of mass of the last tail beads of the surfactants for (a) M1, (b) M2 and (c) M3. The size of M1 increases from 283 K to 330 K to accommodate more surfactants from M3.

M1 increases from 283 K to 330 K to accommodate more surfactants from M3. Similar effects have been observed for M2 in Fig. 8(b). However, the size of M3 decreases from 283 K to 323 K (shown in Fig. 8(c)) which remains constant at higher temperatures (data not shown).

A closer look at Fig. 2(a) and 7(a)–(d) reveal that the population of M2 increases with temperature by incorporating more surfactants from M3. This results in increasing the volume of few micelles at the cost of reducing the volume of others. Increasing the volume leads to increase in Boltzmann energy which has a V dependence

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111

and decrease in surface energy which has a V 2=3 dependence, where V is the volume of the aggregate. Thus, a balance between two competing forces arising from Boltzmann penalty and surface energy leads to a lowering in population of the smallest micelle with increase in temperature. As water content decreases from A1 to A2, a micellar phase is transformed into a lamellar phase at 283 K. From the snapshot of A2 at 283 K in Fig. 1(d), presence of two bilayers are evident where each bilayer is strongly interdigitated. The height modulations of both leaflets are calculated using Eq. 3 and superimposed on Voronoi diagram of area per head group of upper and lower surfaces. Fig. 9(a) and (b) clearly show a height modulation along y axis which is the characteristics of a ripple phase. The area per head group of the upper surface is higher than that of the lower surface. This indicates a larger density of head groups in the lower surface than that in the upper surface consistent with the density profiles of each layers in Fig. 3(b). As one increases the temperature, height modulations get inverted at 338 K for both the surfaces where the one dimensional rippling along y-axis is still preserved. At 345 K, the ripple phase is disrupted where height fluctuations are not observed any more (Fig. 9(e) and (f)). The thickness of both the bilayers shows a little modulation along y axis (see Fig. S2 in Supplementary Materials A) due to a periodicity in the degree of interdigitation at 283 K which slowly becomes uniform as temperature increases (data not shown). However, a bilayer thinning or increase in area per head group is not observed at 345 K. This depicts that the bilayer does not melt at 345 K, although the one dimensional ripple phase is transformed into a LbI phase where surfactants from two layers are still interdigitated.

3.4. Influence of water content on geometric packing It is well known that amphiphilic molecules can self-assemble into a wide variety of structures in presence of water. Different solution conditions such as the pH, ionic strength or temperature can drive a transformation of one structure into another. Transformation of one well defined structure to another is induced by two opposing forces: hydrophobic association and head-head repulsion acting at the hydrocarbon–water interface. These two competing

Fig. 9. Voronoi diagram of the area per head groups in the upper leaflet in (a), (c) and (e) and lower leaflet in (b), (d) and (f) on the contour of their height at 283 K, 338 K and 345 K respectively. The unit of all color bars is in nm. The plot shows the presence of rippling along y axis till 338 K which disappears at 345 K for both leaflets.

109

interactions are responsible to create the optimal area per head group exposed to the water interface which is the driving parameter for the geometric packing consideration. In our simulations, the water content plays the major role in constructing the optimal area per head group required for a specific macro-structure. At higher water content at 283 K, BTMAC and SA molecules selfassemble into 10 spherical micelles of three different sizes which get transformed into two bilayers of similar thickness as the water content decreases. The rdf of head beads of BTMAC and SA with respect to the COM of last tail beads along the chain confirms the presence of spherical symmetry as well as the sizes of three classes of micelles. An approximate radius (R) of the spherical micelle is estimated from the location of the peak of the rdf which is used to calculate its surface area. Using the aggregation number in each spherical micelle (M), the area per head group (a0 ) of BTMAC and SA is calculated using the following expression, a0 ¼ 4pR2 =M. The value of a0 after averaging over three classes of micelles are found to be 1:16  0:31 nm2 and presented in Table 2. The average location of the peaks of head-to-tail distributions of surfactants (see Fig. S2 in Supplementary Materials A) in M1, M2 and M3 can be an approximate critical chain length (lc ) for the micelles which is found to be 2:62  0:12 nm (shown in Table 2). At lower water content, the surface area of the bilayers are easily obtained from their box size which is used to calculate the respective a0 (0:38  0:06 nm2) and shown in Table 2 after averaging over 4 surfaces from two bilayers. The value of a0 decreases once the water content decreases due to the less availability of water molecules surrounding the head groups. The smaller a0 due to the less hydrated head groups couples with an increment in lc . The value of lc for A2 is estimated from the thickness of the bilayer (shown in Fig. S3 in Supplementary Materials A), dðx; yÞ ¼ zu ðx; yÞ  zl ðx; yÞ and found to be 3:14  0:09 nm in A2. Table 2 presents that the value of a0 decreases from A1 to A2 although the value of lc increases as the water content decreases. Fig. 10 represents a schematic diagram to show how different water contents can result into different a0 . The interplay between the coupled area per head group (a0 ) and the critical chain length (lc ) and the hydrophobic core volume of the surfactants can be responsible for the required geometric packing parameter associated with the transformation between the micellar to the lamellar phase. Thus a decrease in a0 due to less water molecules surrounding the head groups of the

Fig. 10. A schematic representation to show how change in area per head group from A1 to A2 leads to different macro-structures.

110

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111

surfactants leads into a transformation in macro-structure from a spherical micellar phase to a lamellar phase. 4. Conclusions In the present work, an all atom molecular dynamics simulation is carried out to study how water concentration influences a phase transformation of a model surfactant/co-surfactant/water system. With higher water concentrations, BTMAC and SA molecules selfassemble into a micellar phase which is transformed into a lamellar phase once the water concentration decreases. Aggregation number of the self-assembled structures are calculated to obtain different classes of aggregates present in both cases. The radial distribution function of the head groups of BTMAC and SA in the micellar aggregate is calculated with respect to the COM of the last tail beads. The symmetric nature of the peak confirms the presence of spherical symmetry of the micelles where the size of the micelles are obtained from the location of the peaks. A density distribution of the surfactant molecules of the lamellar phase shows the presence of overlapping layers at lower water concentrations. The radius of gyration, eccentricity and the order of the micellar phase are lower than that of the lamellar phase due to the spherical nature of the micelles. Different micelles get fragmented or fused to change the size of the micelles once their temperature is increased keeping the water content same. Construction of a Voronoi diagram of BTMAC and SA superimposed on their height fluctuations for the lamellar phase confirms the presence of one dimensional rippling at lower temperature. As the temperature increases, the ripple phase transforms into a LbI phase with an interdigitation present between two layers of a bilayer. Our simulations show that the decrease in water content leads to a lowering in head group areas (a0 ) which is coupled with an increase in critical chain length (lc ) of the surfactant molecules. The complex interplay between the coupled a0 and lc and the hydrophobic core volume drives different packing parameters for different macrostructures. Thus, our simulations reveal that water concentration can be an effective parameter to control the phases of surfactant molecules which can be crucial to obtain their desired functionalities. Acknowledgement This work is carried out with a research grant from SERB, DST (SB/FT/CS-134/2013). We are thankful to IIT Jodhpur high performance cluster facility for computational work. We would also like to thank Satyajit Sahu for careful reading of this manuscript. Appendix A. Supplementary data Supplementary data associated with this article can be found, in the online version, at http://dx.doi.org/10.1016/j.chemphys.2016. 11.014. References [1] A.D. MacKerell, Molecular dynamics simulation analysis of a sodium dodecyl sulfate micelle in aqueous solution: Decreased fluidity of the micelle hydrocarbon interior, J. Phys. Chem. 99 (7) (1995) 1846–1855, http://dx.doi. org/10.1021/j100007a011. [2] D.P. Tieleman, D. van der Spoel, H.J.C. Berendsen, Molecular dynamicssimulations of dodecylphosphocholine micelles at three different aggregate sizes: Micellar structure and chain relaxation, J. Phys. Chem. B 104 (27) (2000) 6380–6388, http://dx.doi.org/10.1021/jp001268f. [3] S.J. Marrink, D.P. Tieleman, A.E. Mark, Molecular dynamics simulation of the kinetics of spontaneous micelle formation, J. Phys. Chem. B 104 (51) (2000) 12165–12173, http://dx.doi.org/10.1021/jp001898h.

[4] L. Maibaum, A.R. Dinner, D. Chandler, Micelle formation and the hydrophobic effect, J. Phys. Chem. B 108 (21) (2004) 6778–6781, http://dx.doi.org/10.1021/ jp037487t. [5] Z. Ma, H. Yu, W. Jiang, Bump-surface multicompartment micelles from a linear abc triblock copolymer: a combination study by experiment and computer simulation, J. Phys. Chem. B 113 (11) (2009) 3333–3338. [6] S. Bandyopadhyay, M.L. Klein, G.J. Martyna, M. Tarek, Mol. Phys. 95 (1998) 377–384, http://dx.doi.org/10.1080/00268979809483170. [7] A.R. Rakitin, G.R. Pack, Molecular dynamics simulations of ionic interactions with dodecyl sulfate micelles, J. Phys. Chem. B 108 (8) (2004) 2712–2716, http://dx.doi.org/10.1021/jp030914i. [8] J.F. Nagle, S. Tristram-Nagle, Structure of lipid bilayers, Biochimica et Biophysica Acta 1469 (2000) 159–195. [9] D.F. Evans, H. Wennerstrom, The colloidal domain where physics, Chemistry, Biology and Technology Meet, VCH, New York, NY, 1994, 10010. [10] J.N. Israelachvili, Intermolecular and surface forces, Academic press, San diego, CA, 1992, 92101. [11] C. Tanford, Theory of micelle formation in aqueous solutions, J. Phys. Chem. 78 (24) (1974) 2469–2479, http://dx.doi.org/10.1021/j100617a012. [12] C. Tanford, Thermodynamics of micelle formation: prediction of micelle size and size distribution, Proc. Nat. Acad. Sci. 71 (5) (1974) 1811–1815. [13] P.A. Pearce, H.L. Scott, J. Chem. Phys. 77 (1982) 951–958, http://dx.doi.org/ 10.1063/1.443871. [14] H.-N. Kim, J. Lee, H.-Y. Kim, Y.-R. Kim, Enzymatic synthesis of a drug delivery system based on polyhydroxyalkanoate-proteinblock copolymers, Chem. Commun. (2009) 7104–7106, http://dx.doi.org/10.1039/B912871A. [15] N. Nasongkla, E. Bey, J. Ren, H. Ai, C. Khemtong, J.S. Guthi, S.-F. Chin, A.D. Sherry, D.A. Boothman, J. Gao, Multifunctional polymeric micelles as cancertargeted, MRI-ultrasensitive drug delivery systems, Nano Lett. 6 (11) (2006) 2427–2430, http://dx.doi.org/10.1021/nl061412u. [16] V. Torchilin, Micellar nanocarriers: pharmaceutical perspectives, Pharm. Res. 24 (2007) 1–16, http://dx.doi.org/10.1007/s11095-006-9132-0. [17] Osada Kensuke, Christie R. James, Kataoka Kazunori, Polymeric micelles from poly(ethylene glycol)? poly(amino acid) block copolymer for drug and gene delivery, J. R. Soc. Interface 6 (2009) S325–S339, http://dx.doi.org/10.1098/ rsif.2008.0547.focus. [18] AnneL. Mrkbak, Peter Degn, Wolfgang Zimmermann, Deinking of soy bean oil based ink printed paper with lipases and a neutral surfactant, J. Biotechnol. 67 (1999) 229–236, http://dx.doi.org/10.1016/S0168-1656(98)00179-5. [19] C. Mulligan, R. Yong, B. Gibbs, Eng. Geol. 60 (2001) 371–380, http://dx.doi.org/ 10.1016/S0013-7952(00)00117-4. [20] A.M. Johannessen, K. Spildo, Enhanced oil recovery (eor) by combining surfactant with low salinity injection, Energy Fuels 27 (10) (2013) 5738– 5749, http://dx.doi.org/10.1021/ef400596b. [21] T. Tadros, Future developments in cosmetic formulations, Int. J. Cosmet. Sci. 14 (3) (1992) 93–111, http://dx.doi.org/10.1111/j.1467-2494.1992.tb00045.x. [22] M. Ramanathan, L.K. Shrestha, T. Mori, Q. Ji, J.P. Hill, K. Ariga, Amphiphile nanoarchitectonics: from basic physical chemistry to advanced applications, Phys. Chem. Chem. Phys. 15 (2013) 10580–10611, http://dx.doi.org/10.1039/ C3CP50620G. [23] B.C. Stephenson, K. Beers, D. Blankschtein, Complementary use of simulationsand molecular-thermodynamic theory to model micellization, Langmuir 22 (4) (2006) 1500–1513, pMID: 16460068. [24] G. Landazuri, V.V.A. Fernandez, J.F.A. Soltero, Y. Rharbi, Kinetics of the sphereto-rod like micelle transition in a Pluronic Triblock copolymer, J. Phys. Chem. B 116 (2012) 11720–11727, http://dx.doi.org/10.1021/jp3009089. [25] S. Hayashi, S. Ikeda, Micelle size and shape of sodium dodecyl sulfate in concentrated sodium chloride solutions, J. Phys. Chem. 84 (1980) 744–751, http://dx.doi.org/10.1021/j100444a011. [26] C. Müller-Goymann, Physicochemical characterization of colloidal drug delivery systems such as reverse micelles, vesicles, liquid crystals and nanoparticles for topical administration, Eur. J. Pharmaceutics Biopharm. 58 (2) (2004) 343–356, http://dx.doi.org/10.1016/j.ejpb.2004.03.028. [27] G. Martini, L. Ciani, Electron spin resonance spectroscopy in drug delivery, Phys. Chem. Chem. Phys. 11 (2009) 211–254, http://dx.doi.org/10.1039/ B808263D. [28] H.J.D. McManus, Y.S. Kang, L. Kevan, Electron paramagnetic resonance and proton matrix electron nuclear double resonance studies of N, N, N[prime orminute], N[prime or minute]-tetramethylbenzidine photoionization in sodium dodecyl sulfate micelles: structural effects of added alcohols, J. Chem. Soc., Faraday Trans. 89 (1993) 4085–4089, http://dx.doi.org/10.1039/ FT9938904085. [29] N. Basilio, L. García-Río, M. Martín-Pastor, NMR evidence of slow monomer micelle exchange in a Calixarene-based surfactant, J. Phys. Chem. B 114 (14) (2010) 4816–4820, http://dx.doi.org/10.1021/jp910984g. [30] W. Shinoda, R. DeVane, M.L. Klein, Coarse-grained molecular modeling of nonionic surfactant self-assembly, Soft Matter 4 (2008) 2454–2462, http://dx.doi. org/10.1039/B808701F. [31] P. Posocco, M. Fermeglia, S. Pricl, Morphology prediction of block copolymers for drug delivery by mesoscale simulations, J. Mater. Chem. 20 (2010) 7742– 7753, http://dx.doi.org/10.1039/C0JM01301C. [32] A. Amani, P. York, H. de Waard, J. Anwar, Molecular dynamics simulation of a polysorbate 80 micelle in water, Soft Matter 7 (2011) 2900–2908, http://dx. doi.org/10.1039/C0SM00965B.

R. Lunkad et al. / Chemical Physics 483–484 (2017) 103–111 [33] N. Thota, Y. Ma, J. Jiang, Molecular insights into the self-assembly of short amphiphilic peptides FmDn and FmKn, RSC Adv. 4 (2014) 60741–60748, http://dx.doi.org/10.1039/C4RA10571K. [34] W. Shinoda, R. DeVane, M.L. Klein, Coarse-grained force field for ionic surfactants, Soft Matter 7 (2011) 6178–6186, http://dx.doi.org/10.1039/ C1SM05173C. [35] J.A. Marqusee, K.A. Dill, Solute partitioning into chain molecule interphases: Monolayers, bilayer membranes, and micelles, J. Chem. Phys. 85 (1) (1986) 434–444, http://dx.doi.org/10.1063/1.451621. [36] A.V. Verde, D. Frenkel, Simulation study of micelle formation by bile salts, Soft Matter 6 (2010) 3815–3825, http://dx.doi.org/10.1039/C0SM00011F. [37] M. Sammalkorpi, M. Karttunen, M. Haataja, Structural properties of ionic detergent aggregates: a large-scale molecular dynamics study of sodium dodecyl sulfate, J. Phys. Chem. B 111 (2007) 11722–11733, http://dx.doi.org/ 10.1021/jp072587a. [38] D. Yordanova, I. Smirnova, S. Jakobtorweihen, Molecular modeling of Triton XMicelles: force field parameters, and partition equilibria, J. Chem. Theory Comput. 11 (2015) 2329–2340, http://dx.doi.org/10.1021/acs.jctc.5b00026. [39] B. Aoun, V.K. Sharma, E. Pellegrini, S. Mitra, M. Johnson, R. Mukhopadhyay, Structure and Dynamics of Ionic Micelles: MD Simulation and Neutron Scattering Study, J. Phys. Chem. B 119 (2015) 5079–5086, http://dx.doi.org/ 10.1021/acs.jpcb.5b00020. [40] C. Anzo, A.H. de Vries, H.-D. Hltje, D.P. Tieleman, S. Marrink, J. Phys. Chem. B 107 (35) (2003) 9424–9433, http://dx.doi.org/10.1021/jp0348981. [41] A.H. de Vries, S. Yefimov, A.E. Mark, S.J. Marrink, Proc. Nat. Acad. Sci. U.S.A. 102 (2005) 5392–5396, http://dx.doi.org/10.1073/pnas.0408249102.

111

[42] A. Debnath, K.G. Ayappa, V. Kumaran, P.K. Maiti, The influence of bilayer composition on the gel to liquid crystalline transition, J. Phys. Chem. B 113 (27) (2009) 9183–9196, http://dx.doi.org/10.1021/jp901551d. [43] P. Mark, L. Nilsson, Structure and dynamics of the tip3p, spc, and spc/e water models at 298 k, J. Chem. Phys. A 105 (2001) 9954–9960, http://dx.doi.org/ 10.1021/jp003020w. [44] H.J.C. Berendsen, D. van der Spoel, R. van Drunen, Gromacs: A message-passing parallel molecular dynamics implementation, Comput. Phys. Commun. 91 (1995) 43–56, http://dx.doi.org/10.1016/0010-4655(95)00042-E. [45] E. Lindahl, B. Hess, D. van der Spoel, Gromacs 3.0: a package for molecular simulation and trajectory analysis, J. Mol. Model. 7 (2001) 306–317, http://dx. doi.org/10.1007/s008940100045. [46] B. Hess, C. Kutzner, D. van der Spoel, E. Lindahl, Gromacs 4: algorithms for highly efficient, load-balanced, and scalable molecular simulation, J. Chem. Theory Comput. 4 (3) (2008) 435–447, http://dx.doi.org/10.1021/ct700301q. [47] W. Humphrey, A. Dalke, K. Schulten, VMD – Visual molecular dynamics, J. Mol. Graphics 14 (1996) 33–38. [48] A.R. Leach, Molecular Modelling: Principles and Applications, Pearson Education Limited, Essex, England, 2001. [49] H.J.C. Berendsen, J.P.M. Postma, W.F. van Gunsteren, A. DiNola, J.R. Haak, Molecular dynamics with coupling to an external bath, J. Chem. Phys. 81 (8) (1984) 3684–3690, http://dx.doi.org/10.1063/1.448118. [50] P.P. Ewald, The calculation of optical and electrostatic grid potential, Annalen der Physik (Leipzig) 64 (1921) 253–287, http://dx.doi.org/10.1002/ andp.19213690304.