Interactions between C₆₀ and vesicles: a coarse-grained molecular dynamics simulation

Abstract

The interaction of fullerene with biological systems has been an issue of great research interest for the past decade. Mechanisms of C₆₀ penetrating and disrupting cell membranes have been widely investigated but are not fully understood. Here we report on coarse grained molecular dynamics (CGMD) simulations of the translocation of monomeric C₆₀ and a fullerene pair across a model DPPC (dipalmitoylphosphatidylcholine) vesicle. Our simulations suggest that monomeric C₆₀ tends to dwell between the head groups of the inner leaflet of the vesicle. This characteristic can be verified from the PMF profiles and the Lennard-Jones interaction energy analysis. The fullerene pair can enter the vesicle membrane as a whole, then disaggregate on a nanosecond timescale inside the membrane. To study the toxicity of C₆₀ on the vesicle, the interactions of the fullerene aggregations with the vesicle are also examined in our simulations. Small fullerene aggregations can penetrate into the vesicle membrane, and do not cause significant mechanical damage to the vesicle membrane. However, as the size of the aggregations increases to greater than the thickness of the vesicle membrane, this can induce a change in the structure of the vesicle membrane or even lead to rupture of the vesicle. Our simulations describe the mechanism for the interactions of C₆₀ with the vesicle, and point out the potential toxicity of fullerenes on the vesicle. These results may provide a useful blue print for drug or gene delivery, and improve our understanding of fullerene cytotoxicity.

---

I. Introduction

Carbon-nanoparticles (CNPs), especially fullerenes, carbon nanotubes and graphene, have great potential and are of scientific and commercial merit because of their unique characteristic properties.^1–3 Investigations of the chemical and physical properties of CNPs, including their size, hydrophobicity, dimensionality, electronic configurations and so on, have yielded promising information and make them an appealing subject in medicinal chemistry and biotechnological applications.^1,4–8 Fullerenes, which are one of the first instances of CNPs, have been developed to be one of the most promising nanomaterials and good candidates for applications in the biological and medical fields.^5,9–11 However, with large-scale production and application of fullerenes, the problem of their potential toxicity towards biological systems is becoming more and more serious.^12,13 High concentrations of fullerenes, which can accumulate in water, were observed to be harmful to the health of people and animals.^14–16 Studies on the interaction of fullerenes with cell membranes may provide a key to understanding the basic questions of their pharmacological and toxicological effects.

Recently, in order to gain a better understanding of the interactions and dynamics of nanoparticles (fullerene is mainly considered here) with biological membranes, several computational techniques including all atom molecular dynamics (MD),¹⁷ coarse-grained (CG) methods,¹⁸ dissipative particle dynamics (DPD) simulations¹⁹ and others, have been used.^20,21 The first choice of computational techniques is all atom molecular dynamics. All-atom methods treat every atom explicitly, which enables a description of structural evolution within atomistic resolution. However it is computationally expensive, so it has been limited by the small sizes of the systems, which usually include small fragments of a lipid membrane and a couple of NPs.^22,23 Li et al. studied the stability of a fullerene pair in a POMF bilayer based on all atom MD simulations. They indicated that the bilayer is less able to accommodate the larger aggregated fullerene pair than isolated single fullerenes, likely due to the distortion of the bilayer structure.²² Qiao et al. explored the interactions of pristine and functionalized C₆₀ with a DPPC bilayer using all atom atomistic simulations, and determined that pristine C₆₀ can diffuse easily into the bilayer but hydrophilic C₆₀(OH)₂₀ molecules do not significantly diffuse into the bilayer.²³

Compared to atomistic models, coarse-grained models, which group several atoms into a single site and have near-atomic resolution, allow systems to be studied on longer time and length scales.^18,24,25 Coarse-grained methods have been widely used to investigate nanoparticles interacting with membranes.^27,28 Wong-Ekkabut’s group²⁷ reported fullerene translocation through a lipid membrane based on computer simulations. Their results showed that fullerene aggregates can permeate through a cell membrane in a microsecond time scale. They concluded that high concentrations of fullerene have effects on the structural and elastic properties of a lipid membrane, and suggested that mechanical damage is not likely to be responsible for membrane disruption and fullerene toxicity. D’Rozario et al. used CGMD simulations to explore the interactions of C₆₀ and its derivatives with lipid bilayers. They suggested that pristine C₆₀ partitions into the bilayer core, whilst C₆₀(OH)₂₀ experiences a central energetic barrier to permeation across the bilayer.²⁸

DPD solves hydrodynamics equations for materials by representing the materials as discrete particles that can be much larger than individual molecules.^19,29 If the particle volume is chosen to coincide with atom groups in the CG MD model, its resolution is similar to CG MD, although it typically treats problems that are more generic.²⁹ Recently, the DPD method has been used to study the internalization of single nanoparticles³⁰ and the interaction of membranes with multiple nanoparticles.^31,32 Yang and Ma found that the shape anisotropy and initial orientation of a nanoparticle can play a complicated role in the physical translocation of a particle across a lipid membrane via DPD simulations.³⁰ Chen et al. performed large scale DPD simulations to study the interaction between nanoparticles (NPs) and vesicles.³¹ They demonstrated that the adsorption of NPs often induces NP clustering on the outer and/or inner surfaces of vesicles, leading to different vesicle responses, including a change of vesicle morphology, formation of protuberance, and vesicle rupture. Zhang et al. performed DPD simulations to elucidate the transmembrane penetration mechanisms of multiple NPs.³² Their simulations demonstrated that NP translocation proceeds in a cooperative manner, where the interplay of the quantity and surface chemistry of the NPs regulates the translocation efficiency.

More detailed discussion of different computational techniques for investigating the interactions of nanoparticles with membranes can be found in review articles^33–35 and references therein.

Experimental studies supported the results from computer simulations that fullerenes can enter the lipid bilayer^36,37 and modify the bending and compressional modulus of the host lipid membrane³⁸ or induce the consequent contraction of the cell membrane.³⁹ Ikeda et al. confirmed that the C₆₀ units were located in the hydrophobic core of the lipid bilayer membrane in liposomes using differential scanning calorimetry and ¹³C NMR spectra.⁴⁰ More about the simulation and experimental studies of the interaction of fullerenes with membranes, especially the partitioning and solubility of C₆₀ fullerenes in lipid membranes, can be found in a review article⁴¹ and references therein.

Vesicles, which are also called lipsomes when they are prepared artificially, can be considered as a simplified model of biological membranes that allows the investigation of dynamic properties of living cells. Lipid vesicles have been devised to create effective drug delivery carriers which can provide precise spatial and temporal delivery of therapeutic agents to the target site.⁴² They can also be used to describe extensive biological phenomena and thus have become a subject of intensive research both in experiments and theory. The research of vesicles can deepen our understanding of the characteristics of cell membranes and it can also provide an effective way to study the interaction between CNPs and biological membranes.^43–45 In order to explore the fullerene–lipid interaction, Katsamenis et al. studied the effects caused by C₆₀ on a large unilamellar vesicle.⁴³ Their results provided evidence for the successful incorporation of C₆₀ into liposomal bilayers. Zupanc et al. studied the disruption of a lipid bilayer by C₆₀ suspended in water via an experimental method.⁴⁴ Their results showed that C₆₀ aggregates in aqueous suspension and a C₆₀ water suspension can cause changes to the average mean curvature of the lipid membrane, leading to the rupture of POPC vesicles. Bouropoulos et al. investigated the interaction between lipid bilayers and unmodified C₆₀ fullerenes to understand the fundamental problems via physicochemical and theoretical methods.⁴⁵ Experimental studies revealed that the presence of C₆₀ in the membranes induced a distortion in vesicle morphology and an increase in vesicle size. The theoretical calculations supported the experimental findings that high concentrations of fullerenes induce disruption to the surface of the bilayer, which is not the case for low concentrations of fullerenes as the membrane remains intact.⁴⁵

As shown above, in the past twenty years, many studies have been performed on the interactions between fullerenes and lipid membranes. However, several aspects remain unclear, such as the mechanism of fullerene translocation into the lipid membrane and their possible effects on biological membrane. These aspects are very important not only because of the possible toxicity of fullerenes in biosystems, but also because fullerenes may be chosen as a candidate for drug delivery carriers. In past theoretical studies, especially molecular dynamics simulations, the interaction between fullerenes and planar membranes was mainly considered.^22,23,27,28 However, there are few investigations about the influence induced by the interactions of fullerenes with bending membranes (like those in small vesicles, which have high membrane curvature). In experiments, scientists have begun to study the interaction of fullerenes with vesicles and have obtained some useful results.^43–45 However, little is known about the mechanism of fullerene–vesicle interactions in detail via coarse-grained molecular dynamics simulations. In this paper, based on a computer simulation, we describe the mechanism of permeation of C₆₀ and its aggregates through a DPPC vesicle membrane. We built a coarse-grained (CG) model to investigate the interaction between fullerenes and vesicles. Our work provides insights into the thermodynamics of fullerene and its aggregations permeating through the vesicle and points out the possible mechanism for the potential toxicity of fullerene on the vesicle.

II. Models and simulation method

A. Systems studied

The CGMD simulations can acquire an ideal process in atomic scale, which is hard to observe in experimental studies. Firstly, we studied the translocation of monomeric C₆₀ into a small DPPC vesicle membrane using unstrained CGMD simulations to investigate the mechanism of the interaction between the isolated fullerene and the vesicle. Secondly, penetration of the fullerene pair into the vesicle membrane is simulated in the same box with the same vesicle. Initially, both the fullerene monomer and pair are placed about several nanometers from the centre of the vesicle. Furthermore, to test the fullerene toxic effect on the vesicle, we performed simulations of the vesicle with different concentrations of fullerene molecules. The fullerene molecules are initially randomly placed separately in the water outside the vesicle. Performing these unbiased MD simulations could observe the penetration of C₆₀, in the aggregated state, into the vesicle membrane bilayer core, tracking the behavior of C₆₀ inside the membrane, and analyzing the effect of these aggregations on the structure and shape of the vesicle. All of the simulations lasted from several hundreds of nanoseconds to several thousands of nanoseconds and we present the main observations during at least two unbiased simulations.

B. CG models of DPPC vesicle and fullerene C₆₀

The coarse grained force field was employed in our simulations of the vesicle and fullerene. The Martini force field developed by Marrink et al. is one of the most popular coarse-grained models^18,24 which faithfully represents the intrinsic chemical structure and interactions between different molecules like lipids, water, ions and so on. In this model an average of 4 heavy atoms are represented by one interaction centre, each with an effective size of 0.47 nm and a mass of 72 amu.²⁴ A DPPC lipid is described by 12 beads where the choline (NC₃) group is represented by type Q0, the phosphate group (PO₄) is represented by type Qa, and the glycerol group (GL1, GL2) and the hydrocarbon tail are represented by type Na and C1, respectively.²⁴ A spherical vesicle is used in our simulations, which consists of 1267 DPPC lipids: 410 lipids in the inner leaflet and 857 lipids in the outer one. The spherical vesicle used in the simulation is spontaneously aggregated in aqueous solution, and its diameter is about 16.5 nm.

The force field for C₆₀ is the same as that used in ref. 26, which has been successfully used to study the fullerenes’ properties and their interactions with the lipid membrane.^27,46,47 The fullerene model used is based on an approximate four-to-one mapping of carbon atoms onto CG beads, which is formed by placing 16 evenly spaced C-type particles on the surface of a sphere of diameter 0.68 nm, extremely approaching the real diameter of 0.71 nm for a C₆₀ molecule. The CG C₆₀ particles are restrained so that the overall shape of the molecule is organized firmly.

C. PMF calculations via umbrella sampling^48–51 MD simulations

The free energy profile of the C₆₀–vesicle system along the vesicle radial direction can provide more valuable information about the equilibrium position of C₆₀ in the vesicle membrane and characterize the propensity of C₆₀ penetration into the vesicle bilayer. This is usually referred to as the potential of mean force (PMF).⁴⁸ Deep minima in the potential profile indicate the locations where the C₆₀ molecule prefers to stay, whereas a shallow and uniform PMF corresponds to a relatively easy, barrier-free translocation process. In the present work, PMF is obtained from the standard combination of the umbrella sampling protocol and an implementation of the Weighted Histogram Analysis Method (WHAM).⁵¹ Not only do we do the PMF calculation for C₆₀ penetrating into the vesicle but also for C₆₀ penetrating into a planar membrane for comparison. Firstly, we obtain the PMF profile as a function of the centre of mass (C.O.M.) distance between C₆₀ and the vesicle along the vesicle radial direction via umbrella sampling with CG MD simulations. The total separation distance between C₆₀ and the centre of the vesicle (11.0 nm) is divided into 55 windows that are 0.2 nm apart. The C.O.M. of the C₆₀ particle is restrained relative to the centre of the vesicle with a harmonic biasing potential. For each window, 50 ns CG MD simulations are employed. In total, we accumulate ∼3.8 μs of biased simulation data. Secondly, we also derive a 1D PMF profile as a function of the C.O.M. distance between C₆₀ and the planar membrane in the same way. Fifty-eight sampling windows are used along the bilayer normal z with the same separation (0.2 nm) (z = −5.8 to +5.8 nm; the bilayer centre is at z = 0 nm) to sample from bulk water, across the bilayer and into bulk water. A force constant of 1000 kJ mol⁻¹ nm⁻² is utilized in the positional restraints and biasing potential.

D. Simulation protocol

GROMACS 4.5 (ref. 52) is used to perform our MD simulations. The simulation box has a dimension of 24.75 nm in the x-, y- and z-directions respectively. All bonds are constrained in length according to the LINCS protocol.⁵³ This allows the use of an integration step of 0.02 ps in the simulations, whereas the pair list distances are 16 Å. We employed the Berendsen weak coupling method⁵⁴ with a time constant of 0.3 ps, and the pressure is coupled isotropically using the Berendsen algorithm⁵⁴ at 1 bar with a coupling constant of 1 ps and compressibility of 3 × 10⁻⁵. Electrostatic interactions are treated with the particle mesh Ewald method,⁵⁵ and a cutoff of 12 Å is used in the calculation of van der Waals interactions. All simulations are started from the initial water/vesicle system using different seed numbers for the generation of initial velocities. All of these MD simulations are carried out at a constant temperature (298 K). It has been found that the DPPC bilayer is well stable throughout the simulation with the system dimensions well converged.⁵⁶ Periodic boundary conditions are adopted and the box lengths along the x, y, and z directions are allowed to change independently. The VMD software is used to visualize and analyze results.⁵⁷

III. Results and discussion

A. Unstrained MD simulations for monomeric C₆₀

In the first series of simulations, monomeric C₆₀ is initially positioned beside the vesicle, 2 nm apart from the edge of the outer surface (the C.O.M. distance between C₆₀ and the vesicle is about 9 nm). An unbiased simulation confirms that monomeric C₆₀ placed in bulk water spontaneously enters into the lipid bilayer interior rapidly once it touches the outer surface of the vesicle membrane. Fig. 1 shows one case of the simple progress measurement for the distance from the C.O.M. of C₆₀ to the C.O.M. of the vesicle. We can see that C₆₀ remains in water for 110 ns and then passes the lipid head group region of the outer monolayer of the vesicle membrane and penetrates into the bilayer hydrophobic core within 0.1 ns. After penetrating the outer monolayer of the vesicle membrane, for the first 570 ns, C₆₀ stays close to the head groups of the outer monolayer of the vesicle. But later (at about t = 690 ns), C₆₀ jumps towards the inner monolayer of the vesicle and adsorbs to the head groups of lipids in the inner leaflet within 10 ns. Afterwards it dwells within about 1.2 nm off the core centre of the lipid bilayer. The system keeps this state for subsequent simulation and finally C₆₀ stays in the region of the inner monolayer of the vesicle. After penetration, the translocation of C₆₀ outside the vesicle membrane is not observed, in accordance with the large energy barrier gap for transferring C₆₀ from the centre of the lipid membrane to the bulk water outside of the vesicle. Differing from some small hydrophobic molecules, previous works suggested that C₆₀ shows no preference for the centre of the planar membrane bilayer and is more often found displaced off the centre of the planar membrane bilayer.²⁷ It has also been previously observed that the preference of C₆₀ is for a region at a 1 nm distance from the centre of the planar membrane bilayer in atomistic simulations.²³ This is also similar with the results shown in ref. 28 where C₆₀ just fluctuated ±0.25 nm on either side of the centre of the planar bilayer. But in our study we find a new preference for C₆₀ after penetrating into the bilayer vesicle. Similarly, we find that C₆₀ prefers to stay off of the centre of the vesicle membrane, but close to the head groups of the outer monolayer of the vesicle.

Fig. 1 Trajectory of the C₆₀ molecule in the radial direction of the vesicle during unrestrained MD simulations. The red line depicts the C₆₀ trajectory. The green (purple) line depicts the average positions of the head groups of the lipids in the outer (inner) leaflet of the vesicle. The snapshot was rendered using VMD.

B. Free energy profiles

To estimate the free energy profiles as a function of the position of the C₆₀ molecule along the radius direction of the spherical vesicle, we perform PMF⁵⁸ calculations from bulk water inside the vesicle across the membrane to bulk water outside the vesicle. For comparison, we also performed PMF calculations for the planar membrane along the direction of the bilayer normal from bulk water to the membrane core. In Fig. 2 we show the PMF profiles for C₆₀ across the vesicle and the planar membrane. The profile for the planar membrane performed from bulk water to the core of the bilayer is close to the results of previous works in ref. 26. The differences between the minima may arise from the different force fields used in the studies. Overall, our results are also consistent with previous work concerning the equilibrium distance of C₆₀ from the centre of the membrane: C₆₀ was predicted to remain preferentially off the centre of the membrane (∼±1.230 nm), roughly in the middle between the head group region and the bilayer mid-plane.⁴¹ The PMF of the membrane system has two distinct energy minima, one in each leaflet of the bilayer. So in the planar membrane C₆₀ has two equilibrium places, in both sides off the centre of the membrane and near the head groups. Thus C₆₀ prefers to stay around the equilibrium point and cannot move through the mid-plane of the membrane bilayer unless incurring an energetic cost (about 13.69 kJ mol⁻¹). The profile for C₆₀ across the vesicle membrane is shown as the red line in Fig. 2. It is obviously observed that a minimum locates exactly at the interface between water inside the vesicle and the head groups of the inner leaflet lipid (∼±0.6 nm), which allows C₆₀ to tend to stabilize at this place. Inside the vesicle, a small energy peak (∼±4.93 kJ mol⁻¹) appears near the inner leaflet with a distance of about 1.075 nm from the centre of the membrane, which can prevent C₆₀ from entering the bulk water into the vesicle. When C₆₀ starts to enter into the lipid membrane, a great fall appears (∼±85.93 kJ mol⁻¹), thus C₆₀ cannot easily escape from the vesicle to the outside once it penetrates into the vesicle membrane. It can also be found that there is a pronounced energy fall (∼±21.60 kJ mol⁻¹) in the centre of the vesicle bilayer. It suggests that C₆₀ prefers to jump from the outer leaflet to the inner leaflet of the vesicle. The profile is flat between 0 nm to 0.857 nm, which may explain why once C₆₀ molecules enter into the outer monolayer of the vesicle, some of them would keep this state for a while before jumping towards the inner monolayer of the vesicle. For the small vesicle, the membrane is not flat but curved; as a result, the PMF only has one minimum in the bilayer membrane, which is close to the region of the lipid head groups in the inner leaflet of the vesicle. This is different from the case of the planar membrane where the PMF profile has two symmetry minima. As we know, the fullerene prefers to stay in the regions where PMF gets its minima. Therefore, the fullerene finally tends to stay close to the lipid head group region of the inner leaflet of the vesicle. This is consistent with our unstrained CGMD simulations for the penetration of monomeric C₆₀ into the vesicle membrane.

Fig. 2 Potential of mean force (PMF) for C₆₀ as a fuction of its distance from the center of the planar or vesicle membrane. For comparison, the horizontal axis of zero represents the center of the hydrophobic core of the vesicle membrane and the planar membrane. The red (black) line depicts the PMF for the C₆₀–vesicle (C₆₀–planar membrane) system. The green (purple) dash line depicts the phosphate particles of the outer (inner) leaflet of the vesicle.

C. LJ interaction energy

In terms of energy calculations, we used coarse-grained molecular dynamics simulations to study the non-bonded interaction of the bilayer of the DPPC vesicle with fullerenes and employed the 6–12 Lennard-Jones potential function as the force field. Fig. 3 shows the LJ interaction energy between monomeric C₆₀ and the DPPC vesicle. Initially, C₆₀ is placed artificially at a distance about 2 nm from the surface of the vesicle. It takes several nanosecond MD simulations for the system to reach equilibrium. In order to verify the interaction of different parts of the vesicle system with C₆₀, we further investigate the interaction energy of C₆₀ with the inner and outer leaflet lipids of the vesicle, as well as the water molecules inside and outside of the vesicle, and the whole system except C₆₀ from 0 ns to 2000 ns in Fig. 3. At the beginning, as mentioned above, the distance between C₆₀ and the surface of the vesicle is about 2 nm, which is larger than the distance of the energy cutoff for the LJ potential. As a result, the interaction between C₆₀ and the vesicle is not observed and the C₆₀ molecule is only interacting with the water molecules around it when it is far away from the vesicle. While the C₆₀ molecule is penetrating through the outer leaflet of the vesicle, the interaction energy of C₆₀ with the water outside the vesicle is decreasing, and that of C₆₀ with the lipids on the outer leaflet of the vesicle is increasing. Once the C₆₀ molecule enters across the barrier, the outer leaflet lipids dominate the interactions and the effect of the water outside is almost cut off. As is shown in Fig. 3, after about 110 ns, it is obvious that the lipids of the outer leaflet clench C₆₀ with an interaction energy around 300 kJ mol⁻¹, and the interaction energy of C₆₀ and the water outside the vesicle decreases to a relatively weak value around 30 kJ mol⁻¹. The interaction energy of C₆₀ with the outer leaflet lipids of the vesicle gets up to a value of about 300 kJ mol⁻¹ after penetration, which is stronger than that of C₆₀ with water outside the vesicle (about 250 kJ mol⁻¹) before penetration, so it is possible that it is an established trend for C₆₀ to penetrate into the hydrophobic core of the DPPC vesicle bilayer. After penetration, C₆₀ stays and fluctuates in the region of the outer monolayer of the vesicle for about 600 ns as mentioned above. The state of interaction between C₆₀ and the vesicle is not changed too much until C₆₀ goes towards the inner monolayer of the vesicle. Once C₆₀ moves from the outer leaflet lipids to the inner ones, it mainly interacts with the latter with a negligible effect from water inside the vesicle. The interaction energy of C₆₀ with the inner leaflet lipids is about 270 kJ mol⁻¹, which is much higher than that (about 50 kJ mol⁻¹) between C₆₀ and the outer leaflet lipids. C₆₀ prefers to stay and fluctuate in the region of the inner monolayer of the vesicle due to the high attraction from the inner leaflet lipids at the end of our simulation. An analysis of the interaction energy of C₆₀ with different parts of the vesicle systems could be a good tool for explaining the mechanism of penetration during the whole simulating process. It seems that the fullerene prefers to move along the pathway to cause an interaction energy maximization between the fullerene and the vesicle, which can be obviously observed in Fig. 3. When C₆₀ penetrates through the outer leaflet of the vesicle, the total interaction energy of C₆₀ with the vesicle will increase rapidly. The interaction energy continues to increase while C₆₀ moves from the outer monolayer to the inner monolayer of the vesicle, although the increase is small. The cause of the fullerene penetration into the vesicle membrane can be attributed to the fact that the interaction energy of C₆₀ with the outer leaflet lipids of the vesicle is stronger than that of C₆₀ with the water outside the vesicle. The moving of C₆₀ from the outer monolayer to the inner monolayer of the vesicle can also be due to the bigger interaction energy of C₆₀ with the lipids of the inner leaflet and the water inside the vesicle than that of the lipids of the outer leaflet and the water outside the vesicle (see Fig. 3). Fullerene has the most favorable interaction with its surroundings in the region of the inner monolayer of the vesicle, indicating that C₆₀ finally prefers to fluctuate in this region. The results from the interaction energy analysis are consistent with the PMF profile mentioned above. The minimum in the PMF profile is around the position where the interaction energy of C₆₀ with its surroundings reaches a maximum. According to the interaction energy analysis, we can also obviously explain why the predicted position of fullerene is away from the center of the vesicle membrane. Some similar discussions have also been proposed in previous researches about the interactions between C₆₀ and the membrane.^22,27

Fig. 3 (A) LJ interaction energy of fullerene with the other parts of the system (all of the other atoms: blue line; vesicle: black line; water: red line) versus time. (B) LJ interaction energy of C₆₀ with the outer leaflet of the vesicle (black line) and the water molecules outside the vesicle (red line) versus time. (C) LJ interaction energy of C₆₀ with the inner leaflet of the vesicle (black line) and the water molecules inside the vesicle (red line) versus time.

D. Unrestrained MD simulations for fullerene pair

An aggregated fullerene pair is found to be stable in water, based on our equilibrium CG MD simulations lasting for a 10 ns duration. The fullerene pair is also initially placed at a distance of 2 nm from the vesicle surface. It stays in water for about 12.5 ns, then penetrates into the vesicle membrane. The translocation of a fullerene pair into a DPPC vesicle membrane is shown in Fig. 4. The fullerene pair as a whole can penetrate through the outside membrane of the vesicle within 1.0 ns. It stays beside the head group region of the lipids in the outer monolayer of the vesicle for a while. Then the departure of the fullerene pair is observed at about 23 ns. Next, one fullerene of the pair moves towards the inner layer of the vesicle and locates itself close to the lipid head group region of the inner layer of the vesicle at about 37 ns. The trajectory of the two C₆₀ in pair and their distance during the penetration process are shown in Fig. 5. The result confirms that the fullerene pair remains intact in water, even though it is in the process of penetrating into the outer monolayer of the vesicle membrane. However, after penetration, it disaggregates on the nanosecond timescale (about 10 ns) in our simulation. Our CG MD simulations indicate that fullerene dimerization in the vesicle membrane should be significantly less favorable than in bulk water. These results are consistent with former studies.^22,27

Fig. 4 Penetration of a fullerene pair into the vesicle membrane. C₆₀ is shown in red. The head groups of the lipids in the outer (inner) leaflet are shown in green (purple), lipid tails and water molecules are not shown for clarity. The simulation time is indicated in each snapshot.

Fig. 5 (a) Trajectory of the C₆₀ pair in the radial direction of the vesicle. The black and red lines depict the trajectory of two C₆₀ molecules in pair. The green (purple) line depicts the average head group positions of the outer (inner) leaflet of the vesicle. (b) The center of mass distance between two fullerenes in pair as a function of simulation time.

E. Unstrained MD simulations for fullerene aggregations

To test the effect of the fullerene aggregations on the properties of the vesicle membrane, we perform a second series of simulations. Initially, a number of fullerenes (from several hundred to several thousand) are randomly distributed in the whole box outside of the vesicle.

1. Low concentration of fullerene molecules. When the concentration is low, fullerene aggregates with a small size can be formed in the aqueous environment outside of the vesicle. The aggregates can enter into the vesicle membrane smoothly due to their small size, which is smaller than or comparable to the size of the bilayer (about 50 nm). The solution of fullerene aggregates can be found in the vesicle membrane within our simulation time. Along with their low density in the vesicle membrane, fullerenes do not induce serious deformation of the vesicle membrane, which retains its initial intact spherical structure. The detailed process of 150 C₆₀ molecules permeating from the aqueous phase into the vesicle membrane interior is shown in Fig. 6. The C₆₀ molecules rapidly aggregated in water to form clusters in different sizes from several molecules to several tens of molecules within several tens of nanoseconds. The sizes of these clusters are less than the thickness of the vesicle membrane. Then these clusters permeate into the vesicle membrane interior gradually and disaggregate. By 700 ns, all of the fullerenes enter into the vesicle membrane (Fig. 6). During the whole simulation time, the vesicle mostly keeps its spherical shape. This means that when the sizes of the formed clusters are small and there is a low concentration of fullerene molecules, fullerenes are unable to cause mechanical damage in the vesicle membrane. The result is similar to a former study, in which they performed a CG MD simulation of fullerene translocation through a planar membrane, and they suggested that the presence of fullerene at a low or high concentration, at which only small aggregations formed, is unlikely to result in mechanical damage to the lipid membranes.¹⁰ Thereafter, it is found that after penetration, the C₆₀ clusters disaggregate inside the bilayer on simulation timescales similar to other planar membrane simulations, which indicates that when the fullerenes are placed inside the membrane, they do not form stable aggregations.^22,27

Fig. 6 Penetration of 150 C₆₀ molecules into the vesicle. C₆₀ molecules are shown in red. The head groups of lipids in the outer (inner) leaflet are shown in green (blue), lipid tails and water molecules are not shown for clarity. C₆₀ molecules aggregate in water to form clusters first, then they permeate into the vesicle membrane. Permeation starts with the insertion of a single fullerene in the lipid head group region. The simulation time is indicated in each snapshot.

2. High concentration of fullerene molecules. To test the effect of higher fullerene concentrations on the properties of the vesicle membrane, we performed simulations in the vesicle with 1000 to 2000 fullerenes. The initial position of the fullerene molecules is chosen randomly in the bulk water outside the vesicle. When the concentration of fullerenes is high, besides the small clusters formed in water, the fullerene aggregates include some large clusters, which are larger than the thickness of the vesicle membrane. In this case, not all of the fullerene clusters can penetrate into the vesicle membrane completely. Large clusters may partly enter into the vesicle or adsorb on the surface of the vesicle membrane. As a result, they can induce a change in the vesicle membrane structure due to large interactions between them and the lipids of the bilayer membrane. In addition, the density of fullerenes in the vesicle membrane will be highly increased because the number of small clusters entering into the vesicle membrane is large at a high concentration. The fullerenes with a high density in the vesicle membrane can induce expansion of the vesicle and furthermore cause deformation of the vesicle. Fig. 7 shows the CG MD simulation results of the vesicle in aqueous solution with 1500 fullerenes. As is shown in Fig. 7, after 200 ns, most of the fullerenes have entered into the vesicle membrane, except for several large fullerene aggregates, which only partly penetrate into the membrane. The vesicle is induced to deform from a spherical shape to an ellipsoidal shape due to the effects of the fullerenes. Pore formation is also observed at the membrane where the big fullerene clusters stay. In Fig. 7, we can also find that some lipids of the inner layer displace to the outer layer. This means that a rupture of the vesicle may be occurring. The disruption of the vesicle can be observed obviously when the concentration of the fullerene molecules is large enough (the number of fullerenes is 2000, not shown). The results are also consistent with former experimental studies by Bouropoulos et al.⁴⁵ Their studies revealed that the presence of C₆₀ in the membranes induced distortion in the vesicle morphology, resulting in nonspherical vesicles, and they also suggested that C₆₀ can rupture the liposome membrane.

Fig. 7 (a) Side view of the simulation system (1500 C₆₀ molecules and the vesicle) at 200 ns. The red balls denote the C₆₀ molecules, the green balls denote the lipids. (b) Cross section of the simulation system at 200 ns. The red balls denote the C₆₀ molecules, the green (cyan) balls denote the lipid head groups (tails) of the outer leaflet of the vesicle, the blue (yellow) balls denote the lipid head groups (tails) of the inner leaflet of the vesicle. Water molecules are not shown for clarity.

IV. Conclusions

In summary, we use coarse-grained molecular dynamics simulations to investigate the thermodynamics and kinetics of C₆₀ permeation into a DPPC vesicle membrane and their potential toxic effects on the vesicle. Our simulations show that monomeric C₆₀ can easily break into a vesicle bilayer membrane by means of passive transport. These results are in good agreement with experiments.^40,43,45 Taking the overall thermodynamic balance for the transfer of C₆₀ from the aqueous phase into the vesicle membrane into consideration, a free energy of penetration value of −110.73 kJ mol⁻¹ suggests that the whole process is spontaneous. This is similar to the case of a single fullerene translocation into a planar membrane, for which the relative value is about −87.73 kJ mol⁻¹. For the planar membrane, the PMF is nearly symmetric about the core of the membrane bilayer. But the result is different for the vesicle, a spherical membrane. A large decline of free energy from the outer monolayer to the inner monolayer of the vesicle is found. This distinction may be caused by the bending of the vesicle membrane. As a curved membrane, the density of the lipid head groups in the inner vesicle leaflet is larger than that in the outer vesicle leaflet. Compared with the planar membrane monolayer, the inner leaflet lipids of the vesicle squeeze together due to bending and it would cost more energy to break through the pyknotic head group of the DPPC lipids. Relatively, the outer leaflet lipids of the vesicle do not arrange themselves, and the energy cost for penetration is relatively small. The determination of the location of C₆₀ incorporated in the vesicle membrane is very important not only in terms of the design of effective carriers for drug or gene therapy, but also for identifying any potential toxic effects of C₆₀. Our simulations provide a possible location of fullerene in the vesicle membrane, and may provide a basis for further correlation researches of other nanoparticles.

In simulations of the interaction between fullerenes and membranes, most of the systems used are constructed with small fragments of a flat lipid membrane and a couple of fullerenes. This may give some uncertain or incomplete results because of the size effect of the small system. In contrast to the case with a fragment of flat membranes, the advantage of using vesicles is that a vesicle can serve as a source of lipid molecules, and includes the local membrane curvature and asymmetry between the two leaflets of the membrane.^31,59 In our simulations, the fullerene monomer prefers to stay in the region of the head groups in the inner leaflet of the vesicle after penetration, which is different from simulations of fullerene with a flat membrane. Besides, the model we used can also give physical pictures of the shape changes and ruptures of the vesicles induced by fullerene aggregations. Additionally, the former MD simulations did not consider the large aggregation state of fullerenes, which may cause significant differences in the interactions with membrane. In our simulations, we studied the large fullerene aggregations and found that they can induce shape changes or ruptures of the vesicle.

Fullerenes can aggregate in water to form clusters with different sizes. Small clusters (especially as a fullerene pair) can translocate as a whole into the vesicle bilayer easily within a short time. After translocation, the clusters can disaggregate inside the vesicle bilayer membrane. The observation of the disaggregation of the fullerene pair in our simulations is consistent with the MD simulations showing that the pair prefers to depart from each other in the lipid bilayer membrane.²² Our simulations show that the vesicle lipid bilayers reduce the propensity of fullerene to aggregate. The disaggregation of small clusters in a planar membrane was also observed in CGMD simulations.²⁷ Furthermore, our results are significant in relation to studies of the cytotoxicity of C₆₀ and its aggregates. The most likely mechanisms of C₆₀ toxicity are mediated by lipid membranes, as it needs to cross or break biological membranes to manifest toxicity. When the concentration of fullerenes is low, at which only small clusters are formed (their diameters are less than the thickness of the membrane), they cannot cause significant mechanical damage to the vesicle membrane but disaggregate in the vesicle membrane. However, as the concentration of fullerenes becomes high, at which some of the clusters are larger than the thickness of the membrane, the big aggregations can induce a change in the structure of the vesicle membrane, even leading to the rupture of the vesicle. In addition, pores formed on the vesicle membrane, which are induced by big fullerene clusters, are in fact closely related to cytotoxicity. They can induce a toxic effect by destroying the concentration balance of the ions and other molecules between the inside and outside environments of the vesicle.

Acknowledgements

This work is supported by the China Postdoctoral Science Foundation (No. 2012M511447) and the National Natural Science Foundation of China (No. 11204257).

References

1. L. J. dos Santos, G. P. Rocha, R. B. Alves and R. P. de Freitas, Quim. Nova, 2010, 33, 680–693.
2. R. H. Baughman, A. A. Zakhidov and W. A. de Heer, Science, 2002, 297, 787–792.
3. V. Singh, D. Joung, L. Zhai, S. Das, S. I. Khondaker and S. Seal, Prog. Mater. Sci., 2011, 56, 1178–1271.
4. Y. Li, X. Chen and N. Gu, J. Phys. Chem. B, 2008, 112, 16647–16653.
5. T. Da Ros and M. Prato, Chem. Commun., 1999, 10, 663–669.
6. G. Pagona and N. Tagmatarchis, Curr. Med. Chem., 2006, 13, 1789–1798.
7. J. L. Li, B. Tang, B. Yuan, L. Sun and X. G. Wang, Biomaterials, 2013, 34, 9519–9534.
8. F. Wahid, T. Khan, A. Shehzad, M. Ul-Islam and Y. Y. Kim, J. Nanosci. Nanotechnol., 2014, 14, 744–754.
9. T. Da Ros, G. Spalluto and M. Prato, Croat. Chem. Acta, 2001, 74, 743–755.
10. R. Bakry, R. M. Vallant, M. Najam-ul-Haq, M. Rainer, Z. Szabo, C. W. Huck and G. K. Bonn, Int. J. Nanomed., 2007, 2, 639–649.
11. P. Chawla, V. Chawla, R. Maheshwari, S. A. Saraf and S. K. Saraf, Mini-Rev. Med. Chem., 2010, 10, 662–677.
12. T. M. Benn, P. Westerhoff and P. Herckes, Environ. Pollut., 2011, 159, 1334–1342.
13. V. L. Colvin, Nat. Biotechnol., 2003, 21, 1166–1170.
14. C. M. Sayes, J. D. Fortner, W. Guo, D. Lyon, A. M. Boyd, K. D. Ausman, Y. J. Tao, B. Sitharaman, L. J. Wilson, J. B. Hughes, J. L. West and V. L. Colvin, Nano Lett., 2004, 4, 1881–1887.
15. C. M. Sayes, A. M. Gobin, K. D. Ausman, J. Mendez, J. L. West and V. L. Colvin, Biomaterials, 2005, 26, 7587–7595.
16. E. Oberdörster, Environ. Health Perspect., 2004, 112, 1058–1062.
17. D. Frenkel, and B. Smit, Understanding Molecular Simulations: From Algorithms to Applications, Elsevier, Singapore, 2nd edn, 2000.
18. S. J. Marrink, A. H. de Vries and A. E. Mark, J. Phys. Chem. B, 2004, 108, 750–760.
19. R. D. Groot and P. B. Warren, J. Chem. Phys., 1997, 107, 4423–4435.
20. M. Fosnaric, A. Iglic, D. M. Kroll and S. May, J. Chem. Phys., 2009, 131, 3282–3285.
21. A. H. Bahrami, R. Lipowsky and T. R. Weikl, Phys. Rev. Lett., 2012, 109, 188102.
22. L. Li, H. Davande, D. Bedrov and G. D. Smith, J. Phys. Chem. B, 2007, 111, 4067–4072.
23. R. Qiao, A. P. Roberts, A. S. Mount, S. J. Klaine and P. C. Ke, Nano Lett., 2007, 7, 614–619.
24. S. J. Marrink, H. J. Risselada, S. Yefimov, D. P. Tieleman and A. H. de Vries, J. Phys. Chem. B, 2007, 111, 7812–7824.
25. L. Monticelli, S. K. Kandasamy, X. Periole, R. G. Larson, D. P. Tieleman and S. J. Marrink, J. Chem. Theory Comput., 2008, 4, 819–834.
26. L. Monticelli, J. Chem. Theory Comput., 2012, 8, 1370–1378.
27. J. Wong-Ekkabut, S. Baoukina, W. Triampo, I. M. Tang, D. P. Tieleman and L. Monticelli, Nat. Nanotechnol., 2008, 3, 363–368.
28. R. S. D’Rozario, C. L. Wee, E. J. Wallace and M. S. Sansom, Nanotechnology, 2009, 20, 115102–115108.
29. M. Ramezanpour, S. S. W. Leung, K. H. Delgado-Magnero, B. Y. M. Bashe, J. Thewalt and D. P. Tieleman, Biochim. Biophys. Acta Biomembr., 2016, 1858, 1688–1709.
30. K. Yang and Y. Q. Ma, Nat. Nanotechnol., 2010, 5, 579–583.
31. X. Chen, F. Tian, X. Zhang and W. Wang, Soft Matter, 2013, 9, 7592–7600.
32. H. Zhang, Q. Ji, C. Huang, S. Zhang, B. Yuan, K. Yang and Y. Ma, Sci. Rep., 2015, 5, 10525,  DOI:10.1038/srep10525.
33. H. Ding and Y. Ma, Small, 2015, 11, 1055–1071.
34. F. Tian, T. Yue, Y. Li and X. Zhang, Sci. China: Chem., 2014, 57, 1662–1671.
35. A. J. Makarucha, N. Todorova and I. Yarovsky, Eur. Biophys. J., 2011, 40, 103–115.
36. M. Bortolus, G. Parisio, A. L. Maniero and A. Ferrarini, Langmuir, 2011, 27, 12560–12568.
37. S. M. Santos, A. M. Dinis, F. Peixoto, L. Ferreira, A. S. Jurado and R. A. Videira, Toxicol. Sci., 2014, 138, 117–129.
38. U. Jeng, T. L. Lin, K. Shin, C. H. Hsu, H. Y. Lee, M. H. Wu, Z. A. Chi, M. C. Shih and L. Y. Chiang, Phys. Rev. B: Condens. Matter Mater. Phys., 2003, 336, 204–210.
39. E. Salonen, S. Lin, M. L. Reid, M. Allegood, X. Wang, A. M. Rao, I. Vattulainen and P. C. Ke, Small, 2008, 4, 1986–1992.
40. A. Ikeda, K. Kiguchi, T. Shigematsu, K. Nobusawa, J.-i. Kikuchi and M. Akiyama, Chem. Commun., 2011, 47, 12095–12097.
41. G. Rossi, J. Barnoud and L. Monticelli, Phys. Scr., 2013, 87, 058503–058510.
42. V. P. Torchilin and W. Weissig, Liposomes: A Practical Approach, Oxford University Press, Oxford, 2003.
43. O. L. Katsamenis, N. Bouropoulos and D. G. Fatouros, J. Biomed. Nanotechnol., 2009, 5, 416–420.
44. J. Zupanc, D. Drobne, B. Drasler, J. Valant, A. Iglic, V. Kralj-Iglic, D. Makovec, M. Rappolt, B. Sartori and K. Kogej, Carbon, 2012, 50, 1170–1178.
45. N. Bouropoulos, O. L. Katsamenis, P. A. Cox, S. Norman, P. Kallinteri, M. E. Favretto and S. N. Yannopoulos, J. Phys. Chem. C, 2012, 116, 3867–3874.
46. J. Barnoud, G. Rossi and L. Monticelli, Phys. Rev. Lett., 2014, 112, 068102–068105.
47. N. Nisoh, M. Karttunen, L. Monticellicd and J. Wong-Ekkabut, RSC Adv., 2015, 5, 11676–11685.
48. P. Virnau and M. Muller, J. Chem. Phys., 2004, 120, 10925–10930.
49. M. Souaille and B. Roux, Comput. Phys. Commun., 2001, 135, 40–57.
50. G. M. Torrie and J. P. Valleau, J. Comput. Phys., 1977, 23, 187–199.
51. S. Kumar, J. M. Rosenberg, D. Bouzida, R. H. Swendsen and P. A. Kollman, J. Comput. Chem., 1992, 13, 1011–1021.
52. S. Pronk, S. Pall, R. Schulz, P. Larsson, P. Bjelkmar, R. Apostolov, M. R. Shirts, J. C. Smith, P. M. Kasson, D. van der Spoel, B. Hess and E. Lindahl, Bioinformatics, 2013, 29, 845–854.
53. B. Hess, H. Bekker, H. J. C. Berendsen and J. G. E. M. Fraaije, J. Comput. Chem., 1997, 18, 1463–1472.
54. H. J. C. Berendsen, J. P. M. Postma, W. F. Vangunsteren, A. Dinola and J. R. Haak, J. Chem. Phys., 1984, 81, 3684–3690.
55. T. Darden, D. York and L. Pedersen, J. Chem. Phys., 1993, 98, 10089–10092.
56. T. Husslein, D. M. Newns, P. C. Pattnaik, Q. F. Zhong, P. B. Moore and M. L. Klein, J. Chem. Phys., 1998, 109, 2826–2832.
57. W. Humphrey, A. Dalke and K. Schulten, J. Mol. Graphics, 1996, 14, 33–38.
58. M. Makowski, C. Czaplewski, A. Liwo and H. A. Scheraga, J. Phys. Chem. B, 2010, 114, 993–1003.
59. K. Yang, B. Yuan and Y. Q. Ma, J. Phys. Chem. B, 2012, 116, 7196–7202.

---

Footnote

† PACS numbers: 89.75.Hc, 05.40.Fb, 05.60.Cd, 02.10.Yn.

---