Molecular modelling of fullerene C60 and its amino acid derivatives in aqueous medium

The hydration free energy ΔGs of fullerene C60 has been determined using the quantum mechanical methods and continuum solvent models, microscopic free energy molecular dynamics techniques, and compilation of known experimental data. Theoretical calculations using standard parameterization schemes of solvent models have predicted negative values of ΔGs(C60) whereas application of the revised parameters, that take into account special features of carbon structures, shifts the computed ΔGs(C60) to positive values. Reasonably accurate estimates of the hydration free energies of the neutral and anionic states of functionalized fullerenes were obtained from the generalized Born method. The described data have shown that the calculations of hydration thermodynamics of fullerenes apparently require separate parameterization for the carbon atoms met in these nanomolecules.


Introduction
Fullerenes are a widely studied carbon nanomolecules that find numerous chemical and biological applications as separate species and as functional fragments in complex molecular constructs [1−3]. The numerous practical applications of fullerenes determine significant interest in the rational understanding and theoretical prediction of the properties of these molecules in condensed phase and particularly in the aqueous biomolecular medium. The carbon atoms in fullerenes are interconnected in the networks of fused rings and form closed strained spheroid cages lined with the π-electrons of quasi-aromatic centers. Such electron structure and large molecular volume determine specific properties of these molecules in solution. Pristine fullerene have good solubility in the aromatic hydrocarbon, much lower solubility in the aliphatic hydrocarbons and polar organic solvents, and zero solubility in water under usual conditions [4−6]. However, the functionalized derivatives of fullerenes can be soluble in water [3]. In view of the active studies of the biologically active compounds containing pharmacophore of fullerene C60, it is important to have good understanding of the mechanisms of fullerene−water interactions.
Modern theoretical chemistry methods can provide useful and accurate information about different physico-chemical characteristics of molecular systems. The computational techniques can especially be helpful in the studies of molecular properties that are difficult to measure experimentally. Exactly this situation exists for fullerenes where experimental data on hydration thermodynamics are insufficient or not available. The computational modelling of the molecules in solution generally involve either quantum mechanical (QM) calculations with the polarizable continuum model (PCM) of solvent or molecular dynamics (MD) force field (FF) free-energy simulations of microscopic systems. In this work we consider the calculations of the thermodynamic potentials of fullerene C60 in water using the quantum mechanical PCM and free energy MD FF computational approaches. The influence of the special separate parameterization for the fullerene nanomolecules on the results of the modelling is reviewed. In addition, the hydration free energies of selected pentasubstituted amino acid derivatives of fullerene in neutral and multi-charged states are examined using the density functional theory (DFT) calculations with the polarizable continuum and generalized Born models of solvent.

Free-energy molecular dynamics simulations of fullerene C60 in water
The MD simulations of hydrated fullerene were conducted using the microscopic models with the minimal image periodic boundary conditions (PBC) and spherical boundary conditions (SBC) (Figure1). The potential energies of the system were computed using the OPLS force field with the OPLS-AA [7] and Bojan-Steele (B-S) [8] parameters for carbons of C60 and the TIP3P FF parameters for explicit water molecules. The OPLS-AA FF parameters are focused on reproduction of the properties of liquids while the B-S carbon-water FF parameters are based on the data for graphiteoxygen interactions. The fullerene atoms were not restrained in the PBC MD simulations whereas the treatment of solute−solvent system in the SBC simulations involved using the artificial spherical cavity for fullerene. The harmonic restraining potential applied to water molecules in order to create empty spherical cavity was derived using the radial distribution function of water molecules around fullerene [9]. The more detailed description of the FF for fullerene and simulation setup can be found in [9,10]. The hydration free energy of fullerene ΔGs(C60) was determined from the alchemical transformation of fullerene into a dummy cage molecule C60 → d60 in water. The end-point d60 particle does not interact with the surrounding medium. The resulting free energy of mutation was calculated using the statistical mechanical free energy perturbation (FEP) method [11]   (1) Here, ΔU(r) is the perturbation potential, r denote the vector of atom coordinates, T is the absolute temperature, R is the gas constant , ‹...›A is the ensemble average over configurational space of the reference state A. Only parameters for the non-bonded interactions in the solute topology were changed in the perturbed state. In case of the spherical model, the calculation of ΔGs(C60) involved a correction ΔGcav for the free energy of creation of the spherical cavity for C60 The cavitation free energy ΔGcav of C60 was calculated using scaled particle theory [12]. For calculation of ΔGs(C60) we used the cavity radius Rcav=5.414 Å which renders ΔGcav =144.6 kJ mol -1 . Rcav was determined as the sum of the distance of fullerene carbons to the center of C60 and the calibrated atom radius of atoms C [10].
The FEP/MD calculations using the OPLS-AA and B-S FF parameters for carbons predict sufficiently differing values of ΔGs(C60) [10]. The ΔGs(C60) free energies calculated using OPLS-AA have negative values of -50.2±2.6 kJ mol -1 for the SBC and -50.3±1.1 kJ mol -1 − for the PBC model. In case of B-S FF, the computed free energies ΔGs(C60) are positive and equal 33.6±1.7 and 20.1±0.7 kJ mol -1 for the for SBC and PBC models, respectively. We can see in this connection that the 7.5-ns PBC multi-step FEP and 3-ns SBC single-step FEP calculations yield quite coherent results. The calculated with the B-S FF free energies well agree with the data from other MD simulations [13,14]. The FEP/MD simulations have also shown that ΔGs(C60) is determined by the balance between large negative enthalpy and entropy terms [10].

QM PCM calculations of the hydration free energy of fullerene
The quantum mechanical calculations of hydration free energies ΔGs were carried out using the 'united atom' (UA) [15] and 'standard model density' [16] versions of PCM in Gaussian. The computed ΔGs for the solute in water cavity include electrostatic ΔGel and non-polar ΔGnp contributions The electrostatic contribution ΔGel was obtained from the quantum mechanical calculation using integral equation formalism with the solvent reaction field included in solute Hamiltonian. The nonpolar term ΔGnp was determined using empirical schemes. In the PCM(UA) method, ΔGnp includes the solvent cavity ΔGcav, electron dispersion ΔGdisp and repulsion ΔGrep contributions [15]. In the PCM(SMD) method, calculation of ΔGnp is based on the solvent-accessible surface area and surface tension parameters [16]. The calculated ΔGs depend on the parameters of PCM for the solute−solvent interactions and on the level of theory employed. The cavity polarizable dielectric is assembled from the overlapping spheres around the atomic centers of the solute. The atomic radii, that define the solvent cavity in PCM, are of crucial importance for obtaining accurate ΔGs values. For instance, the computed ΔGs(C60) changes from negative to positive values upon increase of the atomic radius R(C) (Fig. 2a). It should be mentioned in this respect that the original UA and SMD models do not contain separate parameters for the atom types met in these nanomolecules.
The rational approach to estimate atomic cavity radius of the fullerene carbons R(C) in PCM considered several possible ways to derive R(C) for fullerene C60 [17]. This involved quantitative description of pyramidal hybridization of carbons in fullerene C60 and derivation of atomic radii from the constant electron density surfaces (Figure 2b) around the solute. The computed isodensity surfaces qualitatively support using spherical cavity in the MD SBC simulations of C60 (Figure 2b). The atomic radii were also calibrated using the thermodynamic cycle of the pKa changes for dihydrofullerene and saccharin in DMSO and in a toluene−DMSO mixture. Overall, all three applied criteria have shown that the R(C) value in C60 exceeds the standard value of 1.725 Å for the aromatic carbon in the UA model [17]. The derived R(C) values from the calibration for pKa are 1.842 Å for the HF method and 1.864 Å for DFT B3LYP [16].
The calculated values of ΔGs(C60) are displayed in Table 1. The DFT B3LYP PCM(UA) calculations predict negative ΔGs of C60 of -3.7 kJ mol -1 for standard value of R(C) in the united atom model and positive free energy ΔGs = 14.8 kJ mol -1 for the recalibrated for fullerene R(C). The PCM(SMD) method produces negative free energy ΔGs = -31.0 kJ mol -1 . The largest contribution to ΔGs(C60) comes from the non-polar term although there is also present a small polar contribution (Table 1).

Generalized Born calculations of hydration of fullerenes
The free energies of solvation were also determined from the calculations of the solute interactions with a continuum water dielectric using generalized Born (GB) treatment. The GB approach finds broad applications in the FF and QM studies of solvation of small molecules and proteins [18,19]. In this work, the calculations of ΔGs of fullerenes were performed by means of the GB method with the fixed atomic charges and with linear dependence of the non-polar ΔGnp contrinution to the solventaccessible surface area (SA) method [18]. We consider Eq. (3) where the polar ΔGel and non-polar ΔGnp terms are expressed as the sums over solute atoms In these equations, qi are the partial atomic charges, fij GB is the function of interatomic distances and the Born radii of atoms, SAi is the solvent-accessible surface area of atom i, σi is an atomic solvation parameter, ε is the dielectric constant of water [18]. In the calculations of ΔGs with Eq. (3−5), we have used the Amber-95 atomic radii and σi =30 J mol -1 (for all atom types). The employed Mulliken charges qi were determined in gas phase at the B3LYP/6−31G(d)//B3LYP/6−31+G(d,p) level. The side groups in the studied compounds 1,2 ( Figure 3a) were considered in a single extended conformation optimized in gas phase for uncharged state. The anionic states correspond to deprotonation of the carboxylic acid of side groups. In 2, the oxidation of the more distant from the fullerene core −COOH group was considered. In computations of hydration free energies we used the partial atomic charges from the QM calculations, but the partial charges from the standard FFs can also be used in GB−SA to assess ΔGs of the model structures. The computed GB−SA value of ΔGs(C60) is 17.3 kJ mol -1 ( Table 1) that is close to the result from the PCM(UA) model with the re-calibrated carbon radius. The hydration free energy of C60 in the GB−SA method totally comes from non-polar contribution (Table 1). We have also explored energetics in water for pentasubstituted amino acid fullerene derivatives 1, 2 (Figure 3a). Hydration thermodynamics of these compounds are determined by the solute−water interactions with the central hydrophobic cage of fullerene and the peripheral hydrophilic groups.
Step by step deprotonation of the side groups is accompanied by a steep reduction of ΔGs (Figure 3b). Good liner correlations (with the slopes of the plot and correlation coefficients about unity) are observed between the data from the GB−SA model and more accurate PCM(SMD) method (Figure 3b). These results show that despite the simplified treatment of the solute−solvent system the GB−SA method can be useful in the studies of hydration of complex nanomolecules containing fullerene fragments.

Compilation of experimental data
The hydration free energy of fullerene C60 has not been measured experimentally (to our knowledge) due to extremely low solubility of this compound in water [5]. However, there are estimates of the aqueous solubility Saq(C60) of fullerene C60 based on the Hildebrand−Scatchard model and parabolic extrapolation of the solubilities in normal alcohols [5]. The reported value of Saq(C60) is 1.310 −11 ng ml −1 [5,6]. There are also experimental data for the vapour pressure Pvap(C60) of pure fullerene in equilibrium between the pure solid and gas phase (logPvap(C60)=−22.6) [6]. Consideration of the triple system 'solid C60 -gas state of C60 -water solution of C60' allows to determine ΔGs(C60) from the thermodynamic data for transitions solid↔gas and solid↔water. The standard hydration free energy ΔGs° of fullerene was determined [10] In this equation we use standard conditions of T=298 K and 1 M concentration and the pressure P° of an ideal gas of 24.45 atm. The described treatment predicts negative value of ΔGs°(C60)= −24.5 kJ