Atomistic Perspective on Biomolecular Adsorption on Functionalized Carbon Nanomaterials under Ambient Conditions.

The use of carbon-based nanomaterials is tremendously increasing in various areas of technological, bioengineering, and biomedical applications. The functionality of carbon-based nanomaterials can be further broadened via chemical functionalization of carbon nanomaterial surfaces. On the other hand, concern is rising on possible adverse effects when nanomaterials are taken up by biological organisms. In order to contribute into understanding of interactions of carbon-based nanomaterials with biological matter, we have investigated adsorption of small biomolecules on nanomaterials using enhanced sampling molecular dynamics. The biomolecules included amino acid side chain analogues, fragments of lipids, and sugar monomers. The adsorption behavior on unstructured amorphous carbon, pristine graphene and its derivatives (such as few-layer graphene, graphene oxide, and reduced graphene oxide) as well as pristine carbon nanotubes, and those functionalized with OH-, COOH-, COO-, NH2-, and NH3+ groups was investigated with respect to surface concentration. An adsorption profile, that is, the free energy as a function of distance from the nanomaterial surfaces, was determined for each molecule and surface using the Metadynamics approach. The results were analyzed in terms of chemical specificity, surface charge, and surface concentration. It was shown that although morphology of the nanomaterial has a limited effect on the adsorption properties, functionalization of the surface by various molecular groups can drastically change the adsorption behavior that can be used in the design of nanosurfaces with highly selective adsorption properties and safe for human health and environment.

Introduction

Nanostructured carbon holds a prominent place in material science because of its potential in a virtually endless line of engineering applications, for example, in nanoelectronics and energy technology.[1] From noncrystalline carbon (tetrahedral amorphous carbon), fullerene, graphene sheet, and ribbon to nanotube, nanoengineered carbon has also been used for biomedical applications in cancer diagnostic and therapeutic,[2] biosensing,[3−5] and targeted drug delivery.[6−8] For example, carbon nanotube (CNT) localization near amyloid beta peptide not only increases the specificity of drug delivery but also prevents peptide self-assembly in Alzheimer’s disease.[9] Properties of carbon nanomaterials can be further modified by various functional groups at the surface which can substantially change the interaction of the material with the surrounding molecules which can be further exploited in various applications. Along with the technological use, the nanotoxic hazard associated with the uptake of carbon nanoparticles by biological organisms is a highly debated subject. For example, confinement of CNTs in lysosome through the intracellular vesicle uptake can damage the internal lipid membrane and induce lysosomal permeability.[10] Toxic effects of graphene can be caused by interrupting the hydrophobic protein–protein interaction inside the cell which can lead to the failure in biological function of the cell and even cell mortality.[11]

In recent years, considerable efforts have been made to correlate nanomaterials’ toxicity to their chemical, structural, and morphological characterizations.[12−14] Computational data-driven modeling such as quantitative structure–activity relationship (QSAR) has been used to predict toxic effects of nanomaterials on the basis of their physiochemical properties[15−17] and theoretical molecular descriptors. The binding strength of biomolecules to nanomaterials expressed in terms of binding affinities and adsorption free energies is an essential descriptor in the prediction of biological and environmental effects of nanomaterials.[18,19] Binding free energies of amino acids determine in a large extent composition of protein corona,[19,20] which in turn determines the further fate of nanoparticles in the organism and their possible toxic effects. Because the toxicity mechanism of nanomaterials is highly dependent on the surface[12] and nanoparticles’ surface is often modified by covalently linked molecules, studying the response of functionalized nanomaterials to biomolecules is essential for nanotoxicity assessment.

Experimental quantification of the interaction strength between biomolecules and nanomaterial surfaces in terms of binding free energies is difficult, and existing studies are rather scarce.[21,22] Atomistic computer simulations suggest another route in characterization of bio-nanointeractions. Previously, binding of peptides to carbon nanomaterials, such as graphene and CNTs, has been studied via all-atom molecular dynamics simulations in a number of publications.[21,23,24] The free energy of adsorption for a range of amino acids on graphene sheets has been reported.[25] Comer et al.(26) studied adsorption of small organic molecules on pristine and OH-functionalized CNTs. The influence of oxidation defects on peptide adsorption onto the CNT surface has also been studied.[27] Furthermore, ab initio computations have been carried out to investigate the stability of CNTs covalently functionalized by COOH, OH, NH (for n = 1,2,3,and 4), and other groups at different sites and in various concentrations.[28−30]

In the present work, we have investigated adsorption of small biomolecules on a variety of carbon nanomaterials which include tetrahedral amorphous carbon (ta-C), graphene and its derivatives such as monolayer graphene, few-layer graphene, graphene oxide (GO), reduced graphene oxide (rGO) as well as CNT, and both pristine (nonfunctionalized) and functionalized with hydroxyl, carboxyl, and amino groups. The studied biomolecules represent the basic components of proteins and lipids, particularly they include amino acid side chain analogues and fragments of essential lipids and sugars. We have carried out Metadynamics molecular dynamics simulations computing the potential of mean force (PMF) between the biomolecules and nanomaterial surfaces and determined the adsorption free energies from that. Binding free energies of a given set of small biomolecules to a particular nanosurface can be useful in developing thermodynamic fingerprints for bio-nanomaterial interfacial interactions by providing a set of descriptors which characterize the nanomaterial in a biological media that can be an important component in statistical approaches (QSAR, neural networks, machine learning) to evaluate possible biological effects of the nanomaterial.[31] Furthermore, the computed PMF can be used in coarse-grained models of the complexation of nanomaterials and biological macromolecules, including modeling of protein corona formation. Thus, Power et al.(32) showed that PMFs of amino acids at gold surfaces allow us to build a coarse-grained model to describe protein adsorption on these surfaces, and predict the composition of protein corona.

The simulations of this work were carried out under ambient (biological) conditions, with salinity and solvent conditions that mimic those found in the blood. The effect of functionalizing groups’ concentration on the biomolecular binding strength was investigated as well.

Materials and Methods

Material Models and Force Field

We have developed models of carbon nanomaterials with functionalizing molecular groups attached to their surfaces with covalent bonds. The general design of our models is based on the generalized AMBER force field (GAFF),[33] with a few modifications described in the Supporting Information. GAFF has previously showed good performance in the description of various solution properties of small molecules.[34,35] GAFF parameters exist for molecular fragments that correspond to different hybridizations of carbon atoms (sp, sp2, and sp3) and for the functionalizing group used in our simulations. Thus, pristine graphene and CNTs were modeled as an infinite material composed of benzene (C6H6) fragments with sp2-hybridized carbon atoms (three-folded coordinated carbon). Previously, benzene parameters for graphene were used in a number of simulation works,[23,36] and their suitability was confirmed by ab initio computations.[37] Amorphous carbon materials containing sp3-, sp2-, and sp-hybridized carbon atoms were modeled with GAFF parameters (version 2.11) for the respective types of carbon. We assigned force field types and the corresponding parameters by running antechamber[38]via acpype[39] (for GO, rGO, and all adsorbent molecules) or by homemade scripts.

The groups considered in this work to functionalize CNTs are OH, COOH, COO–, NH2, and NH3+, while for graphene oxide, they are OH and epoxy (C–O–C). The protonation state of a functionalizing group is pH- and pKa-dependent. For example, CNT–COOH terminations are dominant when the environment is strongly acidic while at neutral pH, this group is dominantly presented as CNT–COO–. We considered both protonated and unprotonated forms of carboxyl- and amino-functionalizing groups taking in mind that under some biologically relevant conditions, pH may strongly differ from the neutral value. For example, in phagosomes, which is a key component of a cell clearance system and which can be affected by the uptake of nanoparticles, pH is becoming acidic upon phagosome maturing.[40]

The partial atom charges were determined according to the following principles. Charges of all atoms of pristine CNTs, graphene, and amorphous carbon were set to zero. The net atomic charges of the functionalizing CNT groups were determined with the AM1 bond charge correction (AM1-BCC) method[41] running antechamber[38]via acpype.[39] In these computations, the functionalizing groups are represented by hydrogen-terminated fragments. For example, the hydroxyl (CNT–OH) group is neutral and represented by a water molecule (H2O). The carboxyl (CNT–COOH) group is represented by formic acid (HCOOH) when neutral and the formate ion (HCOO–) when negatively charged. Amine groups (CNT–NH2) are represented by ammonia (NH3) when neutral and ammonium (NH4+) when positively charged. Partial charges of the functionalizing groups on CNTs were then determined by rescaling the charge on the corresponding atom in the fragment to account for missing hydrogen. For GO and rGO models, the density-derived electrostatic and chemical (DDEC6)[42] approach was used to calculate net atomic charges for all atoms. The computed atom charges are indicated in Figure a–d for GO and rGO by colors, and charges for functionalizing groups of CNTs are shown in Figure e. Topology files with force field parameters in Gromacs format for all simulated systems are available as a part of the Supporting Information.

Figure 1

Graphical representation of graphene oxide [side (a) and top views (b)] and reduced graphene oxide [side (c) and top views (d)]. Partial atomic charges are shown by RGB coloring scheme varying between −0.65 to 0.45 amu. (e) Schematic representation of functionalizing CNT groups and partial atomic charges on the functionalizing groups. (f–h) Structural models for tetrahedral amorphous carbon taken from ref (43) and replicated 3 × 3 in the x–y plane. Atoms are colored according to their coordination number, red (two-folded), green (three-folded), and blue (four-folded) coordinated carbon, using ATOMEYE.[44]

A set of 29 biomolecules (as listed in Table and Figure S1 in the Supporting Information) was selected to investigate their adsorption behavior on carbon nanomaterials. The set consists of 18 amino acid side chain analogues (amino acids with the backbone fragment replaced with hydrogen[45]) for all naturally occurring amino acids except glycine and proline, two full amino acids that include protein backbones (glycine and proline), four molecules representing lipid fragments, and one sugar monomer (d-glucose). For amino acid side chain analogues HIS, CYM, and GLU, having pKa within the range from 4 to 10, we have considered both charged and uncharged forms, which makes 3 additional molecules. For histidine, we considered two isomeric forms, HID and HIE. In total, five molecules are positively charged, four molecules are negatively charged, and 20 molecules are electrically neutral. This choice of molecules was made in order to cover most of the principal fragments of biomolecules present in biological environments, such as proteins and lipids. The set covers hydrophobic, polar, aromatic, cyclic, and charged moieties. A set of binding energies of these molecules to a material could represent an essential set of descriptors to be used in in silico characterization of the biological activity of nanomaterials.

Table 1

List of Biomolecules Whose Adsorption Behavior was Investigated in This Worka

code	description	charge (e)
ALA	SCAb of alanine	0
ARG	SCA of arginine	+1
ASN	SCA of asparagine	0
ASP	SCA of aspartic acid	–1
CYS	SCA of cysteine	0
GLN	SCA of glutamine	0
GLU	SCA of glutamic acid	–1
HID	SCA of histidine (δ)	0
HIE	SCA of histidine (ϵ)	0
ILE	SCA of isoleucine	0
LEU	SCA of leucine	0
LYS	SCA of lysine	+1
MET	SCA of methionine	0
PHE	SCA of phenylalanine	0
SER	SCA of serine	0
THR	SCA of threonine	0
TRP	SCA of tryptophan	0
TYR	SCA of tyrosine	0
VAL	SCA of valine	0
HIP	SCA of histidine ion	+1
CYM	SCA of cysteine ion	–1
GAN	SCA glutamic acid (neutral)	0
GLY	glycine (amino acid)	0
PRO	proline (amino acid)	0
CHL	choline group of lipid	+1
PHO	phosphate group of lipid	–1
ETA	etanolamine group of lipid	+1
EST	ester group of lipid	0
DGL	d-glucose	0

Chemical structures are provided in Figure S1 of the Supporting Information.

Side chain analogues.

GAFF parameters for the biomolecules in the set were determined by running antechamber[38]via acpype.[39] Partial atom charges were determined with the AM1-BCC method.[41] Prepared topologies and coordinates that are ready for simulations with Gromacs[46] are available as part of the Supporting Information.

Simulated Systems

A honeycomb model of pristine graphene containing 416 carbon atoms was constructed with initial dimensions of 3.4 × 3.2 nm. We have simulated a single layer in water, as well as two- and three-layer graphene. The latter models were also intended to describe surfaces of large-diameter multiwalled CNTs (MWCNTs). The distance between graphene planes was set to 0.34 nm[47] with the initial orientation of A-B stacking and A-B-A stacking for bilayer and trilayer graphene, respectively.[48] In the case of two- and three-layer graphene, we put two respective layers of graphene at 8 nm vacuum distance between them to resemble models for NRCWE-062 multiwalled CNTs with an outer diameter of 8 nm (more information on characterization of MW-CNT can be found on http://www.enanomapper.net/data), while the other compartment between the layers in the periodic cell was filled by water.

The pristine graphene surface was also used to construct GO and rGO models. GO and rGO consist of different percentages of oxygen-containing functional groups such as hydroxyl (OH), epoxy (C–O–C), carbonyl (CO), and carboxyl (COOH). Hydroxyl and epoxy groups are found on the graphene basal plane and other groups are mostly attached to the carbon atoms on the edges. Because the reduction of graphene oxide is mainly aiming at eliminating epoxy and hydroxyl groups on the plane, we have modeled a hydrogen-terminated graphene plane with randomly distributed epoxy- and hydroxyl-functionalizing groups onto the plane to resemble the model for GO and rGO. A more elaborated structural model of GO by considering defect formation and out-of-plane distortion[49,50] can be investigated in the future work. The sites for hydroxyl and epoxy groups were chosen on the top of the carbon atom and between two carbons (bridge), respectively, as reported in the Lerf–Klinowski model.[51] In our model, a required number of hydroxyl and epoxy groups were randomly distributed on the basal plane (with the same amount in both sides) to have a C/O ratio (the fraction of carbon atoms to oxygen containing groups) in agreement with the experiment (see Table for chemical composition). The range of C/O for GO is reported as 1.3–2.7,[52−55] while the reduction of GO increases the ratio to 2.8–10.3 in the rGO surface. We have simulated finite fragments of GO and rGO, with carbon atoms on the edges terminated by hydrogen atoms. The structures were optimized using the Hartree–Fock (HF) method and 3–21 G* basis set with polarization functions on heavy atoms. The optimized structures of GO and rGO surfaces used in the simulations are shown in Figure a–d.

Table 2

Chemical Characterization of all Simulated Systems Representing Different Carbon Nanomaterials

material	functionalization	chemical composition
tetrahedral amorphous carbon (ta-C)		C1944
graphene		C416
bilayer graphene		C1664
trilayer graphene		C2496
graphene oxide (GO)	–OH, –C–O–C–	C336 (OH)86 (O)39 H54
reduced graphene oxide (rGO)	–OH, –C–O–C–	C336 (OH)20 (O)14 H57
CNT		C660
CNT–OH–low	–OH	C660 (OH)19
CNT–OH-high	–OH	C660 (OH)76
CNT–COOH-low	–COOH	C660 (COOH)5
CNT–COOH-high	–COOH	C660 (COOH)77
CNT–COO–-low	–COO–	C660 (COO–)5
CNT–COO–-high	–COO–	C660 (COO–)20
CNT–NH2-low	–NH2	C660 (NH2)10
CNT–NH2-high	–NH2	C660 (NH2)79
CNT–NH3+-low	–NH3+	C660 (NH3+)10
CNT–NH3+-high	–NH3+	C660 (NH3+)20

Atomistic models for tetrahedral amorphous carbon (ta-C) were taken from Deringer et al.(43) The authors of that work generated and characterized a representative library of optimized tetrahedral amorphous carbon slabs of different sizes by applying machine learning combined with DFT modeling.[43,56] They showed that 216-atom ta-C slab models (with 1.13 nm length size) are the smallest system size that can represent the ta-C local structure correctly in which the slab interior is structurally similar to diamond (the prototype for sp3 carbon) and the slab surface is similar to graphene (the prototype for sp2 carbon). Because ta-C does not have an ordered unit cell, three randomly selected slab models composed of 216 atoms were chosen for this study (Figure f–h). We replicated the system 3 × 3 in the x–y plane to have a comparable system size with other studied surfaces. Larger-sized amorphous carbon models reported recently[57] can be used for more realistic surface structures in the future work.

A single-walled CNT (SWCNT) with armchair symmetry expressed by the indices (11, 11) was generated using scikit-nanoa. This pristine CNT was used as a scaffold for functionalization with −OH, −COOH, −COO–, −NH2, and −NH3+ groups. Adsorption sites were chosen on top of carbon atoms of the CNT unit cell (Figure e), as ab initio calculations have shown this to be the preferred adsorption site for all functionalizing molecules considered in the present work.[28] We programmed an automated procedure to functionalize the CNTs based on the open-source network x Python library.[58] Except for the pristine CNT (which corresponds to zero concentration of the functionalizing groups), two concentrations were considered for each functionalized CNT surface: low and high. The low value corresponds to typical experimental conditions (a few wt %) while the high value is the maximum concentration allowed while keeping an intact CNT. This threshold has been previously determined with ab initio calculations.[28] Our numerical results indicate maximum concentration slightly below 25% (in fraction of six-membered rings) for all functionalizing groups. We take this value as the high concentration (see Table for chemical composition of each simulated system).

For all simulated systems, a suitable number of counterions were added (Cl– or Na+ depending on the charge) to neutralize the system. Once neutral, an additional number of the same ions were added to yield salt concentrations at 0.15 mol dm–3, which corresponds to physiological conditions. At start, all ions were added at random positions outside the nanomaterial. The adsorbing molecule was added 0.8 nm away from the surface. The systems were then solvated by adding an appropriate number of water molecules outside the nanomaterial surfaces and not overlapping with the inserted molecule and ions. The Packmol program[59] was used to generate the starting configurations by putting all parts together as described above. A summary of the studied nanosystems and their chemical composition is given in Table .

Free Energy Computations

We have used classical molecular dynamics with enhanced sampling metadynamics algorithm[60] to calculate adsorption profiles and adsorption free energy for small biomolecules on carbon nanomaterial surfaces. The surface separation distance (SSD) between the adsorbing molecule and material surface was used as a collective variable. In all flat nanomaterials (graphene, few-layer graphene, GO, and rGO), SSD is measured as the z-component of the distance between the center of mass of the surface atoms (outermost layer in the case of few-layer graphene) and center of mass of the adsorbate. For CNTs, we modified the definition of s to be adapted to its cylindrical geometry, as followswhere is the difference between the distance from the center of mass of the adsorbing molecule to the CNT axis and R0 is the CNT radius, for which (11, 11) armchair symmetry CNT is 0.75 nm. For amorphous carbon, because of the existence of the surface roughness (typically ∼0.12 nm for ultrathin 2 nm-thick ta-C films[61]), we considered SSD as the minimum distance between all atoms of the ta-C slab and the center of mass of the adsorbate, modified to be a function with a continuous derivative using ZDISTANCES directive in PLUMED plug-in v 2.5.[62]

The potential of mean force is the free energy of the molecule as a function of the SSD and can be written as a histogram over the ensemblewhich, in principle, is computed by counting the number of times the molecule is found at different s values in unbiased MD simulations. This approach is hampered by sampling, particularly, for strong binding surfaces when the molecule stays at the surface most of the simulation. Enhanced sampling is thus needed to observe enough binding/unbinding events and obtain a converged adsorption profile. In metadynamics, a history-dependent bias composed of the Gaussian function is added to the potential which enforces the system to explore regions that were not visited previously.[60] The potential of mean force can be constructed either from the accumulated bias potentials[60] or by force integration.[63]

In this work, we have employed two versions of metadynamics: the adaptive Gaussian well-tempered metadynamics[64] (AWT-MetaD) with the histogram estimator of the free energy, as it was used in the previous work of our group for simulations of adsorption of amino acids at the TiO2 surface,[45] and constant Gaussian height metadynamics with force integration (MetaDF) to calculate the adsorption profile. It was shown previously[63,65] that combining metadynamics and thermodynamic integration can accelerate the convergence of PMF calculation. In MetaDF, the force acting on SSD (i.e., derivative of the energy over the SSD without including the bias potentials) is used to calculate the potential of mean forcewhere s runs over the relevant SSD values and ⟨F(s)⟩ is evaluated from averaging of the force acting on SSD from configurations having SSD in a small range around s, not including contribution from the bias potential. In simulations of this work, AWT-MetaD was used for CNTs while MetaDF for graphene, graphene oxide, and amorphous carbon. The relevant metadynamics parameters in the case of AWT-MetaD were bias factor 15 and an initial height of Gaussians 2.5 kJ/mol, while in MetaDF, a constant Gaussian height of 0.001 kJ/mol was used. In MetaDF, a low value of the Gaussian height was selected in order to provide close-to-equlibrium condition of the metadynamics simulation, while allowing us to accumulate bias in reasonable time to provide good sampling over the collective variable. We have also carried out comparative studies of AWT-MetaD and MetaDF in the case of adsorption on the graphene oxide surface which demonstrated equivalence of the two methods and also better performance of MetaDF compared to AWT-MetaD.

Additionally, in order to reduce convergence time of the simulations, a fictitious wall potential was introduced to prevent the adsorbate to visit states in the bulk region far from the surfacewhere the wall was set at a = 1.5 nm with a force constant of κ = 40 kJ mol–1 Å–4. For finite-size GO and rGO surfaces, an additional cylindrical fictitious wall with a radius of 1 nm from the center of mass of GO and rGO and axis perpendicular to the surface was introduced to prevent the adsorbing molecule from visiting surface edges of GO and rGO. These restrained potential walls limit the region of the space accessible for the adsorbing molecule during the simulation while leaving essential parts of the material surface and bulk solvent region unaffected and thus not affecting the result of simulation. Such an approach, with limiting available space for the molecule of interest to only relevant region, is often used for free energy and PMF calculation with PLUMED,[62] for example, in ligand, binding free energy calculations with funnel metadynamics.[66]

The binding free energy of a molecule is calculated from a converged free-energy profile aswhere δ is the adsorption layer thickness. This parameter must be specified when calculating binding free energies from experimental data[67] as well as by converting the area concentration of bound molecules to a volume concentration that can be compared to the volume concentration of the unbound molecules that are free to diffuse in the surrounding bulk. We used δ = 0.8 nm in our calculations. ΔFads depends weakly (logarithmically) on δ and its exact value is not of high importance for the binding free energy.

Note that the binding free energy defined in eq is that of a single molecule binding to the surface. This situation corresponds to the dilute limit of the Langmuir isotherm, which is the standard model when extracting binding free energies with experimental techniques. Binding free energies are only defined with respect to standard states. Conversions between standard states are carried out by adding the term kBT ln(c/cst), where c and cst are the concentration of the substance in the simulation and in the other standard state, respectively.[45]

All metadynamics simulations in this work were carried out with the PLUMED plug-in v 2.5 62 to Gromacs v.2019.[46] Relevant MD parameters were time step 2 fs, all bonds to hydrogen atoms are constrained, v-rescale thermostate[68] with relaxation time 1 ps with separated temperature control for the nanoparticle and for the rest of the system, cutoff 1 nm, and particle-mesh Ewald summation[69] of electrostatic and Lennard-Jones interactions.

After the start, each system was first pre-equilibrated in semianisotropic NPT ensemble at pressure 1 bar with separate pressure control in the XY plane and in the Z-direction. After that, Metadynamics simulations were run in NVT ensemble with the box sizes equal to the average box sizes of pre-equilibration simulations. In MetaDF simulations with a constant Gaussian height, the first 50 ns of simulations were excluded from average force computations. Further details on convergence of the simulations and evaluation of uncertainty of the computed data are provided in the corresponding section of the Supporting Information.

Results

Adsorption of Biomolecules on Pure Carbon Surfaces

We begin our analysis from nonfunctionalized carbon surfaces. Potentials of mean force for biomolecules on graphene and pristine CNT are shown in Figure . Adsorbents with aromatic groups (TRP, TYR, PHE, and HIE) or cyclic structure (PRO and DGL) show stronger binding than other molecules, and they are getting closer to the surface without facing any positive energy barrier observed for some other, mostly nonpolar molecules at 0.6 nm SSD. Representative configurations of aromatic biomolecules adsorbed on the graphene surface show π–π stacking interaction with the graphene surface (Figure ). Previously, the π–π stacking interactions have been reported to play an important role in the interaction between proteins and carbon nanomaterials.[70] Also, arginine shows strong binding because of its flat sp2-hybridized moiety joining amino groups. A similar picture of biomolecular binding is observed for CNTs.

Figure 2

Adsorption profile (PMF) of different biomolecules on graphene (black) and pristine CNT (red).

Figure 3

Molecular orientation of strongly absorbed biomolecules on the graphene surface. Pictures were taken from the top view at the distance between the center of mass of the adsorbing molecule and graphene corresponding to the free energy minimum, which is found at 0.34 nm for TRP, TYR, PHE, and HIP, 0.35 nm for ARG, and 0.44 nm for DGL.

The calculated binding free energies shows that TRP is the strongest binder to graphene (−24.3 ± 0.5 kJ mol–1), followed by TYR (−21.4 ± 0.3 kJ mol–1), ARG (−17.2 ± 0.4 kJ mol–1), DGL (−16.1 ± 0.3 kJ mol–1), PHE (−15.6 ± 0.1 kJ mol–1), HIP (−14.4 ± 0.3 kJ mol–1 and similar values for HIE and HID), and PRO (−14.0 ± 0.3 kJ mol–1). The weakest binding free energies are for ALA (−1.2 ± 0.1 kJ mol–1) and ASP (−2.3 ± 0.1 kJ mol–1) molecules. All binding free energies values can be found in Figure , and numerical values are provided as a part of the Supporting Information. Similar trends in interaction of amino acid side chains and peptides with graphitic interfaces have been reported.[21,23,25,71]

Figure 4

Binding free energies of different biomolecules on different carbon nanomaterials.

To analyze how the surface curvature could affect the adsorption of biomolecules on carbon nanomaterials, a comparison between the adsorption profile of biomolecules on graphene and pristine CNT has been made. As can be seen in Figure , the potentials of mean force along the SSD for biomolecules on graphene and pristine CNT are similar. Aromatic and cyclic fragments show deeper minima on graphene compared to the CNT. The binding free energy difference between graphene and CNT is TRP (−10.8 kJ mol–1), TYR (−8.0 kJ mol–1), PHE (−6.8 kJ mol–1), DGL (−6.1 kJ mol–1), ARG (−6.0 kJ mol–1), HID (−4.7 kJ mol–1 and similar values for HIE and HIP), and PRO (−4.0 kJ mol–1) while for other fragments, the difference is less than −2.0 kJ mol–1. The reason of this is that these biomolecules can lay flat on the graphene surface and form stronger π–π stacking interaction to the graphene surface compared to CNTs with a relatively high curvature. Surface curvature relation on bio-nanomaterial adsorption is further discussed in the Discussion part.

To understand the effect of layers in multiwalled CNTs on the interaction with biomaterials, we have investigated the interaction of biomolecules with bilayer and trilayer graphene models. As shown in Figure S2 of the Supporting Information, only a small increase (about −2 kJ mol–1) in the binding free energy was observed for aromatic and cyclic fragments on the three-layer graphene compared to single-layer.

Computations of PMF and binding free energies have been done also for three different samples of amorphous carbon. Comparison of PMF of biomolecules on amorphous carbon slabs with a graphene surface (Figure ) shows similar adsorption profiles which differ mostly by a shift over SSD. Note that displacement of the PMFs to larger distances for ta-C slabs compared to graphene surface can be due to different ways of SSD measurement, while the shape of the profile is similar. Deringer et al.(43) showed that the surface of ta-C models has 70% similarity to the graphene surface. An almost similar adsorption profile has been observed for biomolecules on three different ta-C slabs, but shallower minima have been observed for PHE, TRP, and TYP biomolecules on one of the ta-C slabs compared to the other ones. This difference causes 5–12 kJ mol–1 changes in the binding free energy of PHE, TRP, and TYP biomolecules for one ta-C slab compared to the other ones (Figures and S3).

Figure 5

Adsorption profile of different biomolecules on graphene and different ta-C slabs. Surfaces are colored as follows. Black, graphene; red, green, and gray are three different ta-C slabs.

Adsorption of Biomolecules on Oxidized Graphene Surfaces

Because GO and rGO surfaces are functionalized by randomly distributed hydroxyl and epoxy groups, computations of the potential of mean force require longer simulation time to sample different binding regions and average interaction of the biomolecule over the surface. Therefore, we extended metadynamics simulation on GO and rGO to 600 ns and for some stronger binders to 1000 ns. As Figure shows, there is a noticeable difference between the adsorption profile of biomolecules on GO and graphene surfaces while PMF of biomolecules on rGO is similar to that on the graphene surface. Our analysis shows that hydroxyl and epoxy groups form a hydrogen bond with water molecules which hinder biomolecules to be closer to the GO surface. For most of the biomolecules, the potential of mean force is shifted to the larger distances on GO compared to the graphene surface, and the free energy profile showed weaker binding than for pristine graphene. There are however several molecules, particularly GLU, ASP, CYM, CHL, and ETA which show stronger binding to GO than to pristine graphene and they can be found on closer distance to the surface compared to nonoxidized graphene. These molecules are able to form hydrogen bonds to epoxy and/or hydroxyl groups which are stronger than hydrogen bonds formed by water. Binding of these molecules at the GO surface with the formation of hydrogen bonds is shown in Figure .

Figure 6

Adsorption profiles of different biomolecules on various carbon nanomaterial surfaces. Graphene is colored in black, graphene oxide in red, reduced graphene oxide in gray, pristine CNTs in blue, functionalized CNTs with low concentration (0.39 wt %) −OH in green, and functionalized CNTs with high concentration (14 wt %) −OH in orange.

Figure 7

Molecular orientation of some charged biomolecules on the graphene oxide surface. Pictures were taken from the top view at 0.34 nm SSD for ASP, 0.35 nm SSD for GLU, 0.38 nm SSD for GLU, CYM, and ETA, and 0.46 nm SSD for CHL. Hydrogen bonds between charged biomolecules and surface are shown using the purple-dashed line.

Analyzing further how binding free energy is changed by functionalizing groups, we note that the presence of hydroxyl and epoxy groups on the GO surface leads to increase in the binding free energy of negatively charged ASP and GLU biomolecules but also of positively charged ETA and CHL lipid fragments. We can conclude that the presence of hydroxyl- and epoxy-functionalizing groups on the graphene surface makes it selective to groups with thr carboxylate group (GLU and ASP) or fragments with positively charged nitrogen (CHL and ETA).

Reduced graphene oxide has much lower density of epoxy and hydroxyl groups and larger areas of unmodified graphene. That is why most of the biomolecules bind to rGO in the same manner and with the same strength as to pristine graphene, as there is no hydrogen bound water in the patches of rGO without epoxy and hydroxyl groups. The molecules which showed stronger binding to GO also show somewhat stronger binding to rGO compared to graphene. The obvious reason for this is a smaller fraction of epoxy and hydroxyl groups to form hydrogen bonds with.

Adsorption of Biomolecules on Functionalized CNT Surfaces

Binding free energies of biomolecules to functionalized CNTs are collected in Figure . The corresponding PMF on hydroxyl-functionalized CNTs is shown in Figure in comparison with other surfaces, and all other PMFs at functionalized CNTs are shown in Figure S4 of the Supporting Information. As can be seen, the adsorption of biomolecules is highly dependent on functionalizing groups, as well as on their concentration. By comparing the result of OH-functionalized CNTs with oxidized graphene surfaces (Figure ), we can see that by increasing the OH concentration in both graphene and CNT surfaces, the PMF is shifted to the larger distances and cause weaker interaction with nanomaterial surfaces for most of the biomolecules (except for GLU and ASP), in a similar fashion as it was observed for GO in comparison with graphene. One can also see in many cases similarities in behavior of PMF at GO and high-density OH-functionalized CNTs.

Figure 8

Binding free energies of different biomolecules on functionalized CNTs with different groups and concentrations.

Functionalizing CNTs with charged amino or carboxyl groups also affects binding to these surfaces. Charged molecules bind particularly strong and show a strong response to concentration. For example, PHO, GLU, and ASP with negative charges show stronger response to the concentration when the CNT is functionalized by −NH3+ compared to −NH2 or to pristine CNTs. Similarly, positively charged LYS, HIP (protonated histidine), and ARG bind stronger to the negatively charged −COO–-functionalized CNT. Thus, we can say that electrostatic interaction plays an important role in response of biomolecules to nanomaterial surfaces. Still, there are strong responses on CNT functionalization observed for neutral molecules or/and neutral functional groups. For example, d-glucose binds very strongly to both positive −NH3+- and negative −COO–-functionalized CNTs, while GAN (neutral form of glutamine acid) binds strongly to COOH–-functionalized CNTs. These specific features of binding are determined by the possibility of forming hydrogen bonds between adsorbents and functional groups of CNTs. Concerning the effect of concentration of the functional group, in almost all cases, low-density functionalization shows response in the same direction as high-density functionalization but of lower magnitude. Values of adsorption free energies obtained at high and low density of functionalizing groups, as well as for pristine CNTs, can thus be used as a reference for interpolation for other values of density.

Discussion

Our results on the potentials of mean force for biomolecules at different carbon nanomaterials reveal that both the nanomaterial surface and chemical nature of biomolecules play key roles in the selective adsorption of biomolecules on the surface. To get more insight on how the chemical nature of molecules affect their binding abilities, we replotted the data on PMFs grouping amino acids into categories (polar, aromatic, hydrophobic, and charged) which is shown in Figure S5 of the Supporting Information. We can see that with increasing hydrophobicity in the nonpolar group, LEU and ILE bind at larger distances compared to ALA and they bind stronger to the graphene while weaker to the GO surface. With increasing polarity of biomolecules in the polar group, GLN and ASN bind stronger to the surface compared to the SER and THR. This selectivity is more noticeable on the graphene and ta-C compared to the GO surface. Biomolecules in the polar group are getting closer to the GO easier, without positive energy barrier, compared to graphene and ta-C surfaces. A larger aromatic amino acid, TRP, binds stronger compared to the other aromatic biomolecules on the graphene and ta-C while it does not show any distinguishable selectivity toward the GO surface. All aromatic amino acids have similar PMF profiles at the GO surface. Selective adsorption of ETA, ASP, and GLU on GO compared to the other members in the same basic or acidic groups has been observed. Water molecules can play profound roles in mediating biomolecular interactions. Figure illustrates the role of water in the interaction of biomolecules with the carbon nanomaterials by showing binding configurations of amino acids in each category (polar, aromatic, etc.) to specific nanomaterials and hydrogen bonds formed between adsorbing molecules, surface, and water. Water molecules are repelled from the graphene surface because of impossibility to form hydrogen bonds to it, thus forming a dewetting layer above graphene that makes it easier for the aromatic (TRP) or nonpolar (ALA) molecules to get closer to the surface. For molecules with polar and charged groups, favorable hydrophilic and hydrogen-bond interactions between interfacial water molecules and biomolecule stabilize the latter at a larger distance to the graphene surface compared to the molecules representing nonpolar and aromatic categories.

Figure 9

Representative configurations of absorbed biomolecules of each category (polar, aromatic, etc.) at the free-energy minimum on graphene (left) and graphene oxide (right). Distance between the center of mass of the adsorbing molecule and surface is 0.35 nm and 0.44 nm for ALA; 0.35 nm and 0.39 nm for SER; 0.34 nm and 0.42 nm for TRP; 0.35 nm and 0.42 nm for ARG; and 0.44 nm and 0.35 nm for GLU on graphene and graphene oxide, respectively. Hydrogen bonds between the absorbed biomolecules and carbon nanomaterial or water are shown in blue and the rest are in orange. Snapshots were taken by the funnel fast script developed for metadynamics simulation.[66]

When the surface is oxidized by hydroxyl groups, they form hydrogen bonds with water molecules which prevent biomolecules to get closer to the surface and dictate special orientation and conformation of biomolecules toward the surface (see, e.g., different orientations of SER and GLU on graphene and graphene oxide). The competing behavior of water with biomolecules for adsorption to the surface has been also reported by previous simulation works.[72,73] For the acidic GLU biomolecule, water mediates the interaction of GLU to the graphene oxide surface by increasing the number of hydrogen bonds and stabilizes it at a closer distance to the surface compared to noncharged molecules. The effect of surface curvature on adsorption of biomolecules on nanomaterials is an often discussed question.[74,75] Jana et al.(75) showed the strong dependence of polypeptide adsorption on the surface curvature of the carbon nanomaterial, with the inner (concave) surface of the CNT adsorbing the peptide most strongly, followed by the planar (graphene) surface, and the outer (convex) surface showing weaker binding. In our work, we have not studied a concave surface, but we addressed to the curvature effect by comparing biomolecule adsorption on graphene and CNT(11,11). Our results for PMF of aromatic molecules show deeper (up to 10 kJ/mol) minima on graphene compared to CNTs which is in line with the interaction trend discussed by Jana et al. Because we have studied the adsorption of individual amino acid fragments on carbon nanomaterials, the binding free energy difference between the graphene and CNT is much less than what Jana et al. have reported for the adsorption of peptides. The authors of ref (75) have also partitioned the interaction energy of the peptide with the nanosurface into contributions from different amino acids and showed that aliphatic and aromatic amino acid side chains react more distinctively to the surface curvature compared to other amino acid side chains. In addition, their QM calculations for benzene and toluene (aromatic representative) revealed 2–3 kcal/mol interaction energy difference on plane compared to convex carbon nanosurfaces. The corresponding value for ethane or methane (aliphatic representative) was reported about 1 kcal/mol and even less than 1 kcal/mol for other amino acid representatives. These conclusions in the energy difference for small molecules at planar and convex surfaces are well in agreement with our results.

Conclusions

In this work, we have used Metadynamics molecular dynamics simulations to study the adsorption behavior of a variety of biomolecules on carbon nanomaterial surfaces. We have computed potentials of mean force and binding free energies for 29 small molecules representing basic building blocks of biomolecules to each of 17 considered carbon-based nanosurfaces differing by morphology and functionalization. We have found that morphology of pure carbon materials affects adsorption behavior only slightly because potentials of mean force of the considered molecules in most of the cases were similar on graphene, multilayer graphene, pristine CNT, and amorphous carbon. Pure carbon nanomaterials were also found to be highly selective to molecules containing aromatic and cyclic moieties, which also show some dependence of binding on the curvature, while polar and charged molecules have generally low affinity to these surfaces. On the other hand, functionalization of surfaces by different groups affects the adsorption behavior of biomolecules very strongly. Particularly, charged biomolecules bind strongly via electrostatic interaction and show a strong response to density of the functionalizing groups on the surface. Possibility of the formation of a hydrogen bond to water and to adsorbing molecules as well as their competition is another factor affecting adsorption. For example, surfaces with a high fraction of polar groups, such as graphene oxide, show generally weaker binding of biomolecules than pristine carbon surfaces because of the water layer attached to these functional groups by hydrogen bonds. However, in some specific cases, polar or charges amino acids or lipid fragments show strong binding to such surfaces. Thus, the modification of carbon nanomaterial surfaces by different functionalizing groups can be used to tune highly selective binding of specific biomolecules to the nanomaterial. Data collected in this work can be used further as guidelines (and in the perspective as a data source for data-driven approaches) to, for example, design surfaces selective for adsorption of certain types of proteins, or to predict biocomplexation around nanoparticles taken up by living organisms in in silico schemes of toxicity assessment.