Dependence of fullerene aggregation on lipid saturation due to a balance between entropy and enthalpy

It is well-known that fullerenes aggregate inside lipid membranes and that increasing the concentration may lead to (lethal) membrane rupture. It is not known, however, how aggregation and rupture depend on the lipid type, what physical mechanisms control this behavior and what experimental signatures detect such changes in membranes. In this paper, we attempt to answer these questions with molecular simulations, and we show that aggregation and membrane damage depend critically on the degree of saturation of the lipid acyl chains: unsaturated bonds, or “kinks”, impose a subtle but crucial compartmentalization of the bilayer into core and surface regions leading to three distinct fullerene density maxima. In contrast, when the membrane has only fully saturated lipids, fullerenes prefer to be located close to the surface under the head groups until the concentration becomes too large and the fullerenes begin clustering. No clustering is observed in membranes with unsaturated lipids. The presence of “kinks” reverses the free energy balance; although the overall free energy profiles are similar, entropy is the dominant component in unsaturated bilayers whereas enthalpy controls the fully saturated ones. Fully saturated systems show two unique signatures: 1) membrane thickness behaves non-monotonously while the area per lipid increases monotonously. We propose this as a potential reason for the observations of low fullerene concentrations being effective against bacteria. 2) The fullerene-fullerene radial distribution function (RDF) shows splitting of the second peak indicating the emergence short-range order and the importance of the second-nearest neighbor interactions. Similar second peak splitting has been reported in metal glasses.

to few hundreds of nanometers 27,[30][31][32] . We have previously shown how fullerenes aggregate in water and how such clusters can passively translocate into lipid membranes 14 . It was further shown that once inside a membrane, the aggregates dissociate and fullerenes become dispersed. It was also observed that the dispersed fullerene molecules stay inside the membrane for simulation times of microseconds with no damage to the bilayer. Free energy analyses have indicated that a dimer is energetically unfavorable due to the distortion of the bilayer structure [33][34][35] . The above explanation, however, holds only at low fullerene concentrations 36 . Several experimental and computational studies have shown that fullerene aggregates are harmful to biological membranes 17,[19][20][21][37][38][39] . However, the experimental results of fullerene aggregation in lipid membrane are still controversial especially since the aggregate sizes are not known 20,[40][41][42] . Importantly, the mechanisms of fullerene aggregation inside the lipid membrane remain poorly understood 14,20,28,37,39,43 . Computer simulations offer a complementary method 44,45 , however, fullerene aggregation may depend on the concentration and the accessible time and length scales set some limitations that necessitate the use of coarse-grained (CG) methods 46 .
In this study, we performed coarse-grained molecular dynamics (CGMD) simulations to investigate the physical mechanisms of fullerene aggregation by varying the concentration of fullerenes in model lipid bilayers. Free energy, enthalpy, and entropy profiles were examined which allowed us to explain the thermodynamics of aggregation. The results show why fullerenes favor dispersion and aggregation in lipid bilayers at low and high concentrations, respectively.

Methodology
We use CGMD simulations to study the aggregation of fullerenes in lipid membranes at different fullerene concentrations. Saturated and unsaturated lipid membranes are modeled with 1,2-dipalmitoyl-sn-glycero-3-p hosphocholine (DPPC), 1,2-dioleoyl-sn-glycero-3-phosphocholine (DOPC), and 1-palmitoyl-2-oleoyl-sn-glycero-3 -phosphocholine (POPC) lipids. The simulations were performed by using the GROMACS 4.5.5 package 47 . The Martini force field version 2.1 was used for lipids and water 48 . The coarse grained fullerene molecule (F16) consist of 16 beads and the latest updated fullerene parameters were used to reproduce the atomistic potential of mean force (PMF) profile of transferring a fullerene through a lipid membrane 49 . This fullerene model has been previously validated against experimental solvation free energies and atomistic simulations and it has been used in numerous independent studies with various lipid bilayers (DPPC, POPC, DOPC, DSPC and DUPC) and shown to work well 14,19,21,28,34,39,42,[49][50][51][52][53] . The models and atom numberings are shown in Fig. S1. The simulated systems consisted of 512 lipids and 16,000 water CG molecules. The fullerene to lipid ratio was varied between 0 and 0.5. All simulations were carried out under the constant particle number, pressure and temperature (NPT) ensemble. Fullerenes were initially randomly placed in the bulk water phase. Over the course of the simulation they spontaneously came into contact with the membrane, and after a few microseconds fullerenes had passively penetrated inside the bilayer. The simulations were run for 20-50 μs with a 20 fs time step. The last 10 μs were used for analysis and errors were estimated based on the standard deviation. Temperature was controlled using the Parrinello-Donadio-Bussi velocity rescale algorithm 54,55 with a time constant of 1.0 ps. The Parrinello-Rahman algorithm 56 with a time constant of 5.0 ps and compressibility of 4.5 × 10 −5 bar −1 was used to keep a constant 1 bar pressure in a semi-isotropic fashion. Periodic boundary conditions were applied in all directions. A cut-off of 1.2 nm was used for the long-range neighbor list. The Lennard-Jones interactions were shifted to zero between 0.9 and 1.2 nm, and the Coulomb potential was shifted to zero between 0 and 1.2 nm. The Visual Molecular Dynamics (VMD) software was used for molecular visualization 57 .
The total simulation time exceeds 500 microseconds. The details of all systems are provided in Table 1, and the data presented in the main text is at the temperature of 298 K. We would like to point out that although the experimentally observed main phase transition of DPPC is at about 314 K 58 , the transition temperature for DPPC within the coarse-grained Martini model has been determined to be 295 ± 5 K 59 . This discrepancy is due to the mapping in the process of coarse-graining, but it is important to notice that although the temperature becomes shifted, the qualitative features remain intact. In addition, it has been shown that modeling lipid systems is a particular strength of the Martini model, see e.g. Refs 46,[60][61][62] . Another related question is how good is the Martini model for describing fullerene aggregation? Monticelli 49 compared dimerization in different CG models and atomistic models, and showed that while some CG models overestimate dimerization (about 5.6 kJ/mol in water in the case of Chiu et al. 63 ), the Martini parameterization performs remarkably well, within 2.9 kJ/mol (favorable) in water and 0.5 kJ/mol (unfavorable) in octane. Thus, the current Martini model for fullerenes 49 provides quite a high resolution in terms of free energy differences. Similarly, Hsu et al. used the Martini model to study the interactions of fullerenes and bacterial membranes, and showed very good agreement with both atomistic simulations and experiments 64 . This is different from the case of protein aggregation in membranes where Martini appears to excessively favor dimerization 65 .

Results and Discussion
Membrane dimensions. Figure 1 shows snapshots of structures at the end of the simulations and Fig. S2 shows final structures at additional concentrations for completeness. As the figures show, aggregation was observed to depend on fullerene concentration in accord with previous studies 28,34 . Fig. 2a,b show that both the area (Fig. 2a) and volume per lipid (Fig. 2b) of DPPC, DOPC and POPC bilayers increase with increasing fullerene concentration. The area per lipid has been studied before and our results are in agreement with previous studies 14,21,28,66 . We also analyzed membrane thickness, Fig. 2c. The thicknesses of both the DOPC and POPC bilayers increase monotonously with increasing fullerene concentration in agreement with previous studies 14, 66 .
In contrast, however, the thickness of the DPPC bilayer shows non-monotonous behavior as it decreases at fullerene concentrations below 20% where it reaches its minimum, after which the thickness starts to increase. This non-monotonous change in thickness appears to be characteristic of fully saturated (here DPPC) bilayers only. We return to the origin of this below in the discussion regarding the free energy, but the top views in Fig. 1 show (see the difference between 20% and 30%) that qualitatively, this is linked to the clustering of fullerenes on   As is clear from Fig. 1, systems with double bonds (POPC and DOPC) have qualitatively different behavior from the fully saturated DPPC system: instead of clustering, fullerenes organize themselves in layers. To characterize this behavior, we measured the component-wise density profiles across the bilayers, Fig. 3 (see in Figs. S3, S4). As the figure shows, at concentrations up to 20%, fullerenes stay close to the surface in DPPC whereas in DOPC fullerenes layer into three different regions beginning at concentrations of 5% and becoming pronounced at 20%, subsequently developing in a systematic fashion. This is also evident in the side views shown in Fig. 1. This layering is induced by the double bonds whose density is also shown in Fig. 3. Although the area and volume per lipid are higher in the DOPC bilayer, the double bonds present a subtle constriction thus dividing the membrane interior into three regions. In the case of DPPC the area per lipid is smaller and at increasing concentrations it appears that after saturating the region close to the head groups, the presence of a larger local distortion becomes more likely. Note that the largest bilayer deformation was observed at 50% concentration of fullerene in DPPC bilayer caused by the accumulation of fullerenes inside the bilayer. To examine changes in elastic properties, we computed the isothermal area compressibility in all the cases (Fig. S5). At zero fullerene concentration the values are in excellent agreement with those obtained by Daily et al. 67 . For concentrations up to 20%, DOPC and POPC show qualitatively similar behavior while compressibility for DPPC increases rapidly and exceeds that of both the POPC and DOPC systems at about 10%. This rapid increase is correlated with decreasing thickness as reported in Fig. 2c. At higher concentrations (30% and above), the error bars become large due to lack of data (correlation time increases). Quantitative analysis of fullerene clusters will be shown in the next section.

Fullerene aggregation.
To study the aggregation of fullerenes in lipid bilayer, we used previously validated software from Barnoud et al. 34 . Two fullerene molecules are considered to be clustered when their centers of masses are closer than 1.30 nm 19,34 . Fig. 4 shows that the largest cluster increases with increasing fullerene concentration and Fig. S6 shows the time dependence of cluster formation in two representative cases. In the case of DPPC, the onset of rapid cluster growth starts around 20%, correlating well with the minimum in the non-monotonous behavior of bilayer thickness in Fig. 2c. Aggregation in DPPC and DOPC (POPC is similar to DOPC) is qualitatively different from each other. As Fig. 4 shows, the cluster sizes are significantly larger in the case of DPPC, whereas in DOPC and POPC, only small clusters are observed; the percentage of fullerenes belonging to the largest cluster is up to ~70% in DPPC (the largest cluster did not span all of the fullerenes) and only 20% in DOPC bilayers (Fig. S7).
To further analyze aggregation, we determined the center-of-mass positions for fullerenes as a function of time in each case and computed the (unnormalized due to the thin film geometry) radial pair distribution function (RDF) for them. In the case of the unsaturated DOPC matrix, the RDF is liquid-like (Fig. 5). This is drastically different in the DPPC matrix (Fig. 5d): While below 20% the RDFs are liquid-like, the second peak of the RDF has fully split at 30%. In the case of POPC, the RDF displays a characteristic broadened shoulder at 40% (the inset in Fig. 5c shows the onset). DOPC with double bonds in both of its tails shows no signs of such behavior even at 40%. These results together with the density profiles in Fig. 3 show that the peak splitting is due to lipid saturation (double bond position is indicated in the figure). Peak splitting is often observed in metallic glasses and liquids. The second peak is due to, in general, second nearest neighbors. As the analysis of Ding et al. shows 68 , splitting arises when more than two atoms share parts of the second-nearest neighbor shell. In terms of their three dimensional Voronoi cells, splitting arises when the second shells share more than a single vertex (a single vertex equals to one atom and leads to liquid-like behavior), that is, they have to share either an edge or a face of their Voronoi cell indicating short range order 68 . This order is due to cluster formation in the DPPC system with fully saturated chains.
Free energy. Figure 6 shows the free energy profiles and the decomposition of free energy into its entropic and enthalpic components. The free energy profiles were calculated using where k B is the Boltzmann constant, T is the absolute temperature, and r g( ) is the radial distribution function with r being the distance between two fullerenes 69 . Figs. 6 and S8 show that the first and second minima of the free energy are always located at 1.0 and 1.5 nm, respectively (Fig. 6a-c); at 1.0 nm two neighboring fullerenes are in contact. The positions of first and second minima are in agreement with the work of Barnoud et al. who studied the interactions of fullerenes and POPC membranes and found dispersion of fullerenes up to concentrations of 18% (the highest concentration in that study) 34 . Figs. 6a-c also show that at low concentrations the global minimum is at 1.5 nm (the second minimum), meaning that fullerenes prefer to be dispersed as monomers. The origin of the second minimum at 1.5 nm is qualitatively similar to what was observed by Cao et al. 70 who studied fullerenes in ethanol solutions. They observed a second (albeit shallow) minimum and traced that to be due to the fullerene-fullerene distance at which one of ethanol's functional groups is able to become inserted between a pair of fullerenes. Here, 1.5 nm is the distance at which a lipid bead can be inserted between two fullerenes. This is demonstrated in Fig. 7 which shows the in-plane density between fullerenes at separations of 1.0 nm and 1.5 nm. The free energy decomposition in Fig. 6 also shows that independent of the system and fullerene concentration, the enthalpic component always shows a minimum at that distance. The above suggests that fullerenes do not damage biological membranes at low concentrations via aggregation and consequent buckling of the membrane. We would like to point out, however, that thinning of the saturated DPPC membrane with simultaneous increase in the area per lipid, may expose the membrane to, for example, leakage. This is consistent with the observations At around 20%, the first minimum at 1.0 nm becomes the global one indicating that aggregation becomes dominant. Furthermore, as the concentration increases from 20% to 30%, in the case of DPPC the second maximum breaks into two smaller maxima separated by a shallow local minimum. This is consistent with the non-monotonicity observed in membrane thickness (Fig. 2c) and the growth of clusters (Fig. 4). No such splitting occurs in the case of DOPC. To study the effect of system size on free energy profiles, we simulated larger systems consisting of 2048 DPPC lipids with 10% and 30% fullerene concentrations. The results did not show any significant differences (Fig. S9).
To understand this better, following MacCallum et al., the free energy of lipids bilayer was decomposed into its entropic and enthalpic components as 72 .
Additional simulations were performed with varying temperatures of 283 and 313 K in order to obtain a reliable estimate for the entropic component. The results are shown in Fig. 6d-f. First, Fig. 6a shows that as the fullerene concentration in a DPPC bilayer increases from 10 to 30%, the dimerization free energy (r = 1.0 nm) decreases from −1.75 to −6.76 kJ/mol. Fig. 6d shows that this is due to the entropic component decreasing faster than the enthalpic component increases. This is in contrast to DOPC where the entropic component dominates; this is similar to POPC bilayers studied by Barnoud et al. 34 and also verified here (data not shown). Thus, the above results clearly show that the aggregation behavior is controlled by the balance between entropy and enthalpy, and that this depends on whether the acyl chains are saturated or not. Meanwhile the free energy of the second minimum decreases only slightly (~0.15 kJ/mol) and the changes in the enthalpic and entropic components compensate each other (thermal energy is about 2.5 kJ/mol). The inset in Fig. 6b shows the free energy difference between the first and the second minima as a function of concentration. The difference increases faster for the fully saturated DPPC. Importantly, the thermal energy corresponds to about 2.5 kJ/mol indicating (together with the fact that the enthalpic component dominates the free energy) that in DPPC, aggregation is preferred at  concentrations above 20%. In DOPC, however, the thermal energy is enough to disperse the fullerenes even at the highest concentrations.
Diffusion. Finally, we also computed fullerene diffusion coefficients using mean square displacements both laterally and in the direction normal to the membrane surface. The data is shown in Fig. 8. Slower diffusion was observed at higher fullerene concentrations in agreement with Sastre et al. 39 . Perhaps surprisingly, the differences in lateral fullerene diffusion in the three cases are very small. Diffusion in the direction normal to membrane surface, however, shows clear differences: diffusion in DOPC and POPC is very similar, but up to 20% concentration the DPPC system behaves differently with the diffusion coefficient being practically independent of fullerene concentration. The plausible origin of this appears to be explained by Fig. 2: while the area per lipid increases for all systems, thickness of the DPPC systems decreases up to 20% fullerene concentration.

Conclusions
We have investigated the aggregation of fullerenes in lipid bilayers using coarse-grained molecular dynamics simulations over the total time of more than 500 microseconds. We focused on two lipid matrices: 1) DOPC in which both of the acyl chains have a double bond and 2) the fully saturated DPPC. Third matrix, POPC, was used as a control. We decomposed the free energy into its entropic and enthalpic components at different fullerene concentrations. We show, to our knowledge, for the first time that in the fully saturated bilayers fullerenes are controlled by the enthalpic term which favors aggregation. This is in contrast to DOPC (and POPC) in which the entropy dominates and favors dispersion. This is manifested in the layer-like ordering of fullerenes inside the DOPC bilayer and, the emergence of splitting of the second peak in the fullerene-fullerene center-of-mass pair correlation function in the saturated DPPC matrix as the fullerene concentration increases. These results have several consequences: First, several experimental studies have shown that fullerenes have antibacterial effects. It has been suggested that this is due to the increased area per lipid. However, as our results and those of others show, the area per lipid increases with increasing fullerene concentration independent of the lipid type. Here we have shown, however, that membranes consisting of fully saturated lipids become thinner at small and moderate fullerene concentrations yet their area per lipid increases monotonously with increasing fullerene concentration. This is consistent and can explain the observations of Lyon et al. 71 who reported that the antibacterial effects against B. subtilis were the strongest at small fullerene concentrations. Their data also shows (albeit not discussed) that the area per volume ratio changed non-monotonically. Here, we focused on the interactions of fullerenes with pure bilayers. One of the interesting questions is what happens when membranes have varying compositions. Sastre et al. 39 studied cholesterol-containing membranes (PC-Chol-fullerene) and reported that cholesterols and fullerenes do not prefer to be in contact with each other and that the presence of fullerenes may suppress the inter-leaflet flip-flop rate. They also simulated raft-like membranes and observed that fullerenes prefer the more disordered matrix. There are several open questions including how the asymmetric distribution of cholesterol and different lipid types in real cellular membranes influence the partitioning of fullerenes (and nanoparticles in general) between the leaflets, and how the local variations in the available free volume (e.g., due to presence of cholesterol) influence the distribution. Many of the issues remain unresolved yet they have very important consequences both for applications as well as real biological systems as discussed, e.g., by Chen and Bothun 73 .
These results also open the questions of how the balance between entropy and enthalpy depends on the number or/and location of double bonds, and if by modulating the positions allows for localization of fullerenes in some particular region of a membrane. Finally, we would also like to emphasize that our predictions are experimentally testable using model bilayer systems.