Comparative Computational Study of Interaction of C60-Fullerene and Tris-Malonyl-C60-Fullerene Isomers with Lipid Bilayer: Relation to Their Antioxidant Effect

Oxidative stress induced by excessive production of reactive oxygen species (ROS) has been implicated in the etiology of many human diseases. It has been reported that fullerenes and some of their derivatives–carboxyfullerenes–exhibits a strong free radical scavenging capacity. The permeation of C60-fullerene and its amphiphilic derivatives–C3-tris-malonic-C60-fullerene (C3) and D3-tris-malonyl-C60-fullerene (D3)–through a lipid bilayer mimicking the eukaryotic cell membrane was studied using molecular dynamics (MD) simulations. The free energy profiles along the normal to the bilayer composed of 1,2-dipalmitoyl-sn-glycero-3-phosphocholine (DPPC) for C60, C3 and D3 were calculated. We found that C60 molecules alone or in clusters spontaneously translocate to the hydrophobic core of the membrane and stay inside the bilayer during the whole period of simulation time. The incorporation of cluster of fullerenes inside the bilayer changes properties of the bilayer and leads to its deformation. In simulations of the tris-malonic fullerenes we discovered that both isomers, C3 and D3, adsorb at the surface of the bilayer but only C3 tends to be buried in the area of the lipid headgroups forming hydrophobic contacts with the lipid tails. We hypothesize that such position has implications for ROS scavenging mechanism in the specific cell compartments.


Introduction
Fullerenes and their derivatives are considered perspective for various applications in medicine and pharmacology [1]. Particularly, C 60 -fullerene and tris-malonic C 60 -fullerene ( Figure 1) are known as free radical scavengers, that have been shown to protect cells from free radicals (including ROS) that can induce apoptotic injuries in vitro [2,3] as well as in different cell types: neuronal cells [4], hepatoma cells [5] and epithelial cells [6]. At the same time, some of fullerene derivatives (for example dendrofullerenes) [7] inhibit HIV-protease, what gives another prospect for their biomedical use. Harhaji et al. demonstrated [8] the antitumor effect of the water suspension of C 60 on glioma cell cultures when irradiated by light: high concentrations of fullerene caused necrosis, while low concentrations stopped proliferation of the cells and eventually lead to autophagy.
Fumelli et al. [9] have demonstrated that carboxyfullerenes protected human keratinocytes from apoptosis induced by ultraviolet-B (UVB) irradiation. Later the possibilities of preventing the neurotoxicity generated by levodopa (L-3, 4-dihydroxyphenylalanine) have been demonstrated as chemical properties of a water-soluble fullerene derivative and ascorbic acid [10]. Wang et al. compared [11] antioxidant activity of the two regioisomeric forms of tris-malonic fullerene ( Figure 1)-C 3 (with the malonic groups localized on one side of the fullerene molecule) and D 3 (with the malonic groups distributed evenly). It was found that C 3 has a higher activity against ROS than D 3 , presumably due to its better interaction with biomembranes. However, the molecular basis of these interactions remains elusive. Particularly, the process of passive transport through cell membranes of C 60 and its derivatives must be elucidated.
Most of the fullerenes get inside the cell by the energydependent mechanism of endocytosis [12]. On the other hand, hydrophobic C 60 molecules can directly interact with the biomembrane and penetrate it in an energy-independent manner. One of the theoretical methods, which allows to tackle this problem, is computational simulation by molecular dynamics. These can be used to study the penetration of small molecules [13] as well as nanoparticles [14] through biomembranes. Several atomistic molecular dynamics simulation studies have been performed that investigate the interactions of C 60 and its derivatives with model biological membranes [15,16,17]. It was shown that water-soluble fullerene derivatives do not penetrate into the membranes as the C 60 does [16,17,18]. In [15] a coarsegrain model was applied to study penetration abilities of multiple fullerenes into different membranes and the authors showed that multiple fullerenes inside the bilayer cause no deformation of the membrane.
Here we report a comparative study of the interaction of C 60fullerene and its amphiphilic derivatives, C 3 and D 3 , with eukaryotic membranes composed of DPPC lipids using molecular dynamics simulations. We performed series of unrestrained MD simulations and calculated free energy profiles as a function of the position of the studied molecules along the bilayer normal. The aim of our investigation is to explain the mechanism of ROS scavenging for fullerene species with different activities. We found out that molecules have different penetration depths into the bilayer what has implications regarding the possible antioxidant effects of C 60 carboxyl derivatives.

Materials and Methods
Four different systems were used for the simulations: a lipid bilayer with one molecule C 60 , with ten molecules of C 60 , with one C 3 molecule and with one D 3 molecule. The model of the eukaryotic membrane consisted of 128 solvated DPPC lipids. The C 60 -, C 3 -or D 3 -molecules were placed in water and did not contact the lipid surface in the beginning of the simulations. In all simulations we used the Berger model for lipids [19] and the SPC/ E water model [20]. The force field parameters for C 60 , C 3 and D 3 were set as follows. The lengths of the two types of C-C bonds in fullerenes were set to 0,139 nm and 0,144 nm according to [21]. Parameters of the van der Waals interactions were taken from [16,17]. The atomic partial charges of carbon atoms in a fullerene were set 0 for the all of them due to the symmetry of the molecule, while the Mulliken partial atomic charges for C 3 and D 3 derivatives were calculated in GAMESS suite (see SI, Table S1) using ab initio Hartree-Fock method with 6-31G** basis set (prefaced by the geometry optimization of the molecules at the same level of theory).
The all of MD simulations were carried out using Gromacs 4.5 [22]. Simulations were carried out in the NPT ensemble: semiisotropic pressure coupling (Parrinnelo-Rahman barostat [23], time constant 2 ps) and constant temperature 323 K (Nose-Hoover thermostat [24,25], time constant 2 ps). Lipids and water were coupled independently to the heat bath. Periodic boundary conditions were applied in all three dimensions. All bond lengths were kept constant using the LINCS [26] algorithm. The time step was 2 fs, long-range electrostatic interactions were treated with the PME algorithm [27] (real space cutoff 1 nm, FFT grid spacing 0.18 nm). The Lennard-Jones potentials were computed by using a cutoff length of 1.2 nm.
Free energy profiles for fullerenes were obtained using the metadynamics approach [28] and the plugin PLUMED [29]. Each of the metadynamics calculations was preceded with an equilibration run. A total of 1.5 ms of molecular dynamics was carried out. We chose a z-component (perpendicular to the membrane plane) of a vector connecting centers of masses of the membrane and the fullerene as a reaction coordinate (called collective variable) for the metadynamics simulations. The following parameters were used for metadynamics: the integration time step was decreased to 1 fs to improve stability of the systems, the time interval between the addition of two Gaussian functions was t = 500 fs, the Gaussian height w = 0.5 kJ/mol, and the Gaussian width d = 0.5 Å . Potential walls were applied to keep the values of the chosen collective variables in the area of a single membrane leaflet with the adjacent water layer, i.e. between 210 and 35 Å . Simulations were carried on until the convergence of the free energy surfaces was achieved after 45 ns. Additional 4 independent trajectories were obtained for each of the studied systems to calculate the average free energy profiles along with standard deviation. Sufficient sampling of the collective variables in metadynamics simulations was checked by plotting their progress along the simulation time ( Figure S1) to confirm adequate sampling of the chosen reaction coordinate. All of the performed MD runs are listed in Table S2. Analysis of the obtained trajectories was done with standard Gromacs 4.5 utilities. Membrane properties calculations were made using GridMAT-MD [30]. Plots of molecular hydrophobic potential (MHP) for the membrane surfaces were built using PLATINUM software [31].

C 60 fullerene rapidly penetrates into membrane
The equilibrium MD simulation of a single C 60 performed in our study showed that during the first nanosecond of the MD trajectory the fullerene was adsorbed in the area of the lipid headgroups. However, in the third nanosecond a spontaneous dive into the tail region was observed (Figure 2A, C), which is consistent with [15,16,17]. Over the following almost 100 ns of the simulation the fullerene remains inside the membrane. Table 1 represents a summary of previous studies of interactions of C 60 and different membranes. According to this table we can conclude that the shape of the PMF profiles for the fullerene interacting with the membrane are quite similar. All methods, which were used to describe interactions between C 60 and membranes, show a low minimum inside the membrane (near its center). Yet the details are different: some simulations show a small barrier for fullerene before entering the tail region. Also it's important to notice that coarse-grain (CG) models underestimate free energy values for C 60 . Different methods of free energy calculations compromise energies: due to the small simulation time it is almost impossible to reach the thermodynamic equilibrium in such inhomogeneous systems as membranes.
We have used the metadynamics approach to calculate the free energy profiles of the fullerenes interacting with the biomembrane, which, to the best of our knowledge, has not been applied to the fullerene-membrane systems before. So we have validated the method against the well-studied case of the C 60 fullerene interacting with the membrane and came up with the conclusion that the metadynamics can be efficiently exploited in case of the C 3 and D 3 derivatives as well. In the figure 2B one can see that the fullerene overcomes the energy barrier of approximately 15 kJ/ mol to pass into the hydrophobic region of the lipid tails. The global energy minimum for the C 60 is located 0.7 nm behind the center of the membrane: if the C 60 is moving closer to the center of the membrane its free energy is increasing by almost 30 kJ/mol. So, the C 60 remains at a certain distance from the precise membrane center in the region of the hydrocarbon chains. To double-check this observation C 60 was placed in the center of the membrane (between C-terminal atoms of lipid tails) and a 50 ns simulation was run. The fullerene in this simulation also migrated to the area located approximately 0,8 nm from the center of the membrane and remained there for the rest of the simulation. Absorption of a single C 60 on the surface of membrane and its ''jump'' to the tails region cause no damage to the membrane (see Table 2). The above-mentioned results appear to be in a good agreement with the previously reported data serving as a good validation of the chosen model.

Cluster of fullerenes causes changes in the physical properties of the membrane upon penetration
In an aqueous medium fullerenes exist only in aggregates because of their hydrophobic nature. In our study we provide calculations of ten fullerenes and a DPPC membrane with heavyatom resolution. To study the permeability of the membrane for a fullerene cluster which exists in real solutions, we assembled a system consisting of 10 fullerenes above the surface of the membrane. At the start of the simulation fullerenes did not contact with each other or with the membrane, but after a few nanoseconds, they aggregated and absorbed in the area of the lipid headgroups. The first fullerene from the cluster ''jumped'' inside the membrane after 3 ns. The rest of the fullerenes penetrated the membrane afterwards, so that eventually, after 100 ns, 9 of 10 fullerenes had dived into the membrane ( Figure 3C). All of them stayed in the membrane during the next 400 ns of the simulation time. We did not observe their disaggregation during simulation time but according to [15] this happens in the microsecond scale. The most probable localization of fullerenes inside the membrane is around 0.5-1 nm from the center. Contrarily to the previously reported coarse-grained MD simulations [15], the membrane is deformed significantly by the addition of a large number of fullerenes: it curves ( Figure 3C, average curvature radius calculated as in [32], Figure S2), the area per lipid head reduces from 69 to 56 Å 2 , on the contrary, thickness increases to 38,7 Å and up to 42 Å at the end of simulation. Such changes are connected to the formation of a fullerene cluster inside the membrane ( Figure 3C) consisting of 8 fullerenes. MHP plots were made for both membrane sides of the two membrane conformations: at the beginning and at the end of the simulation (see SI: Figure S3). We observed that a large hydrophobic area had appeared on the membrane surface from the side where the fullerenes penetrated the membrane ( Figure 3D). It means that when the fullerenes get inside the membrane a small pore appears with the hydrophobic lipid tails being exposed to the solvent. Lipids in this conformation are more vulnerable to ROS. The diffusion coefficients for DPPC molecules do not change upon addition of the fullerenes in the membrane and they are in a good agreement with the experimental data [33]. Pore or micelle formation is not observed though possibly due to insufficient simulation time. The order parameters for acyl chains S CH ( Figure 3A and 3B) for lipid chains sn-1 and sn-2 were calculated using the following equation: From these plots one can see that the order parameters do not change significantly when only one fullerene is inside the membrane while they modestly increase in the simulation with 9 fullerenes inside the membrane.
Our observations indicate that the addition of multiple C 60 molecules does not cause a phase transition of the lipid bilayer but leads to its loosening and bending as well as to the formation of hydrophobic pores, what altogether can increase membrane susceptibility to ROS and explains toxicity of high concentrations of the C 60 , which was reported previously [34].   We have carried out equilibrium MD simulations of the C 3 and D 3 fullerenes interacting with the membrane along with metadynamics simulations in order to obtain a complete picture of molecular and energetical details of the membrane insertion processes for the C 3 and D 3 variants (Figure 4, 5).
The fully deprotonated forms of both studied fullerene derivatives should be the most probable species in aqueous  environment (pK a 1 = 2.3 and pK a 2 = 5.7 for malonic acid in water). This seems to be no longer the case in the membrane hydrophobic region, where pK a s of acidic groups shift towards higher values. Higher pK a s values should favor protonated species and reduce thus the free energy cost for the translocation of the C 3 /D 3 to the region of the lipid tails. On the other hand in the region of the zwitterionic lipid head groups pK a of malonyl groups depends on the interplay of the local interactions of the fullerene derivatives with phosphate or choline groups. However, we assumed that these effects play secondary role at the stage of the C 3 /D 3 absorption to the bilayer and performed simulations of the fully ionized C 3 /D 3 species (charged 26) only.
In the equilibrium simulation with C 3 , the C 3 molecule adsorbed to the membrane surface in two steps: in the beginning, it was localized at a distance of 2 nm from the COM of the membrane and retained its solvent shell (which consisted of 40 water molecules) as seen in Figure 4 and, later, it moved closer to the COM of the membrane by approximately 0.5 nm and established a hydrophobic contact with the lipid hydrocarbon tails as well as ionic bonds with the lipid headgroups ( Figure 4). This notable behavior of the C 3 derivative takes root in its stereochemistry-the molecule has two hemispheres: one is hydrophobic and the opposite carries three malonic acid residues and thus is hydrophilic ( Figure 1C). We calculated the orientation of C 3 molecule in respect to the membrane normal along the equilibrium trajectory of the C 3 absorption to the membrane surface. Since the absorption process was completed after 25 ns, the C 3 upheld a specific orientation with its hydrophobic  hemisphere pointing to the middle of the membrane and the malonic groups interacting with the head groups of the lipids remained accessible to the solvent ( Figure S4). Energetically, this orientation of the C 3 corresponds to the global minimum of the free energy profile.
For both the malonic derivatives of C 60 we also calculated their solvent accessible surface area (SAS). Upon the absorption of C 3 to the membrane the SAS of C 3 decreased significantly ( Figure S5). On the other hand, the SAS of D 3 remained almost the same during the whole simulation time, what indicates the lack of a close hydrophobic contact with the lipid tails in this case.
The observed orientation of C 3 upon its absorption to the membrane surface appears to play a crucial role in efficient scavenging of ROS and radical protection: It resides not inside the membrane but near the lipid tail region, preventing chain peroxidation of the lipids. The isomeric form of C 3 -D 3 -does not have [35] such high cytoprotective properties since its three malonic acid residues cover the fullerene surface symmetrically ( Figure 1B) and D 3 only adsorbs efficiently to the region of the lipid headgroups where the distance from the COM of the membrane equals 2 nm and retains its solvent shell (see Figure 5C).
Neither the C 3 nor the D 3 stereoisomers penetrate the membrane because of the large free energy penalty of being in the middle of the membrane for these molecules ( Figure 4A and 5A).

Conclusions
The processes of penetration and accumulation of C 60 fullerene, its cluster and the two derivatives of C 60 in lipid bilayers were studied using equilibrium MD and metadynamics approach. Absorption of all of the studied fullerene species on the DPPC membrane occurs within the first nanoseconds. C 60 fullerene (a single molecule and a cluster of ten molecules) spontaneously penetrates the membrane and remains inside of it throughout the whole simulation time. After penetration into the eukaryotic membrane, C 60 stays at a distance of 0.7-0.8 nm from the center of the membrane. A cluster of fullerenes deforms the membrane causing its curvature and formation of hydrophobic regions, which presumably make the membrane more susceptible to the ROS attacks.
The studied tris-malonic derivatives of the C 60 fullerene do not penetrate membrane but rather are accumulated at its surface. However, during adsorption to the membrane the C 3 , in contrast to the D 3 , appears to be deeply buried in the area of the lipid headgroups of the membrane forming specific hydrophobic contacts with the hydrocarbon tails of the lipids along with the ionic interactions with the head region and hydrophilic contacts with the solvent. This energetically favorable orientation of the C 3 derivative, as we assume, accounts for its high ROS scavenging activity comparing to the D 3 . Figure S1 Progress of the collective variable during the C 60 , C 3   Text S1 Creating Models of C 3 and D 3. (DOCX)