Computational analysis of amino acids’ adhesion to the graphene surface

. The mechanisms of cellular growth have attracted scientists’ attention for a long time, leading to recent eﬀorts in establishing cellular growth on speciﬁc functionalized substrates. In order to fully understand the supported cellular growth mechanisms, one needs ﬁrst to comprehend how individual amino acids interact with the substrate material as cells are known to attach to surfaces through speciﬁc proteins designed to improve adhesion. In this study, we have considered graphene as a candidate material for support-assisted cellular growth and simulated the interaction of all 20 naturally occurring amino acids deposited on graphene. Investigations utilized classical molecular dynamics (MD) for amino acids in aqueous solution and in vacuo, in tandem with quantum chemical calculations. The MD simulations were carried out for classical and polarizable CHARMM force ﬁelds. The simulations performed with the polarizable force ﬁeld conﬁrmed that adhesion of amino acids to the graphene surface may be signiﬁcantly enhanced due to the polarization forces, which was further supported by quantum chemical calculations. The performed analysis thus revealed the role of polarization on amino acids’ adhesion to the graphene surface.


Introduction
Living organisms are connected through the extracellular matrix (ECM), which apart from providing structural support is also believed to be important in initial growth and differentiation of stem cells [1][2][3]. ECM consists of a variety of fiber-like shaped structures, such as e.g. fibrils of the protein fibronectin schematically shown in Figure 1 (as green bars), which are responsible for binding of the cells together [4]. The ECM consists of the plasma membrane with integrin domains, which in turn connect polysaccharide fibers; a schematics of this arrangement is illustrated in Figure 1.
Since there is an increasing interest in artificially growing biological tissue, even entire organs, from stem cells, it is crucial to identify materials that can be used to mimic the ECM and guide stem cell growth and differentiation. Graphene is among the materials of interest for this purpose [5][6][7][8] due to its useful properties, as it is on one side very strong and flexible, while at the same time, its physicochemical properties can be tuned through coating of the surface with oxides or organic molecules [9]. Previous attempts of measuring survivability of cells or growing cells atop pristine graphene layers or coated graphene layers have delivered promising results, with most of the experiments showing graphene to be as good a growth material as the body while in some cases even accelerating cell growth [10][11][12][13][14]. The present investigation aims to populate the bio-engineering toolbox with ingredients for intelligent functionalization of graphene through patterned coating with amino acids. Specifically, we have employed a palette of computational tools to model the interactions between graphene and amino acids which in turn determine the binding and unbinding rates to the graphene surface.
Specifically, we have employed molecular dynamics (MD) jointly with the quantum chemistry (QC) investigations to systematically investigate all the 20 naturally occurring amino acids deposited atop graphene surface. The analysis of the simulations revealed diverse binding configurations in some cases, while it also demonstrated a clear indication of non-binding for some of the amino acids. The combined MD and QC study revealed the degree of polarization of the graphene surface that is expected to arise once charged amino acids are put atop, significantly enhancing the binding affinity of the corresponding amino acid. The results of the present investigation are seen important for the further analysis of graphene/protein interaction as it may suggest specific coatings, which will significantly enhance the binding of the protein to the pristine graphene surface.

Classical molecular dynamics simulations
The structures for the 20 different amino acids were taken from the Protein Data Bank (rcsb.org [15]). Classical MD simulations were carried out for the molecules distributed atop a layered graphene structure, where three 99.49Å×202.77Å graphene layers were built in VMD [16]. The multiple graphene layers were used to make the amino acids/graphene interface as realistic as possible; such a multilayer structure ensures little artifacts in graphene distortion due to amino acid deposition compared to if only one sheet was used. The atoms of the lower graphene layer were kept frozen to further ensure little artifacts in graphene distortion.
Every simulation included 18 identical amino acids distributed atom of the upper graphene layer, and the procedure was repeated for all of the 20 natural amino acids. The graphene layers with deposited amino acids were neutralized in water using NaCl by assuming the height of the simulation box equal to 45Å. The protocol was largely adopted from earlier MD simulation on similar systems [17][18][19]. Particle mesh Ewald method was used to treat the long-range electrostatics [20,21] with a cutoff distance of 12Å. CHARMM36 force fields for proteins and water [22,23] was used to model the interatomic interactions in the amino acids, the graphene layer and in between the two subsystems as well as the water interactions. The carbon atoms of the graphene layers were treated as aromatic carbon atoms, described in the CHARMM force field [24], having a partial charge of 0.0e. Two different types, (i) and (ii) of simulations were carried out.
(i) Soluted amino acids atop graphene were placed in a periodic boundary box of 100.691Å×204.184Å×45Å to ensure periodic continuation of the system box, thus eliminating boundary artifacts. The system was equilibrated by first restricting all molecules of the system but the water molecules by a harmonic potential. As the energy of the system decreased through the first phase of the equilibration process more atoms in the system were gradually released from the harmonic restraints with the sidechains of the amino acids being released at the end of the equilibration procedure. The simulation assumed a pressure control with 1 bar and the temperature control at 300 K [25,26]. A temperature of 300 K was used to simulate a natural environment while preventing the water from freezing. (ii) Amino acids were put in the simulation box above the graphene surface in vacuum and simulated at 10 K with temperature control. A small temperature of 10 K was used to avoid clustering of the amino acids in the course of a simulation. Periodic boundary conditions were also employed in these simulations with the simulation box set to 100.691Å×204.184Å×30Å.
The production simulations for all the 20 amino acid types in the two different simulation setups were carried out for 10 ns using a 2 fs integrator time step.
In the case of soluted amino acids, the average binding energy between the different amino acids and the graphene surface was calculated as where E total is the total energy of the system as derived from the simulation, N = 18 is the total number of amino acids in a simulation, E aa (i), E gra and E solv are the internal energies of the i'th amino acids, the graphene layers and the solvent (water + ions), respectively. E aa/solv (i) and E gra/solv are the respective interaction energies between the i'th amino acids and solvent molecules and the graphene layers and water molecules.
In the case when deposited amino acids are studied in vacuum, the average binding energy was defined as

Polarizable molecular dynamics simulations
To study effects of induced graphene polarization by amino acids, we have also performed polarizable molecular dynamics (PMD) simulations. Note that this type of simulations for deposited amino acids atop graphene have not been carried out before. For the PMD simulations we relied on the Drude model [27][28][29][30][31][32][33][34], which is schematically illustrated in Figure 2 for the case of a water molecule. The Drude model is implemented in the NAMD software [27], and requires the introduction of two virtual particles. Both of these new particles are illustrated in Fig. 2. A: Example of a water molecule as described by the Drude model includes a lone pair particle and a Drude particle. The Drude particle is added at the same coordinate as the oxygen atom with oxygen shown as a transparent sphere around the Drude particle. B: Classic model of a water molecule used in conventional (classical) MD simulations. Figure 2 in the case of a water molecule. The Drude particle is close to the host particle, to which it is bound harmonically through a large force constant k D . In the presence of an electric field E the Drude particle will oscillate around a displaced position d = Eq D /k D with q D being the partial charge of the Drude particle. With the defined partial charge of the Drude particle and the force constant, the polarizability of the host atom can be defined as.
In the simulation the Drude particle shares the mass of the host particle and in this paper it is assumed to have a mass of 0.4 of a hydrogen mass and has a charge of the opposite sign to the host atom. The Drude particle and its host atom make a pair, which in the presence of an electric field E defines an induced dipole, µ = q 2 D E/k D , making the molecule, by definition, polarizable. In the dynamics, the Drude particle-host atom pair are governed by a special thermostat, which is set to a small temperature T * such that the pair only has a low amount of kinetic energy. The performed simulations assume T * equal to 1 K.
The Drude model also introduces the off-centered massless interaction site called the lone pair (LP) particle. This LP particle is a massless charged particle introduced to avoid the limitations of the centrosymmetric-based Coulomb interactions. The analogy for more complex molecules is straight forward as the Drude and lone pair particles are simply introduced for every atom heavier than hydrogen. The LP particle has a total of three guide atoms which will guide the LP particle through the simulation to ensure the charge displacement will be consistent over the simulation and are all defined for all oxygen in this study.
In the present investigation all amino acids were converted to the Drude approximation using the Drude prepper script [35]. The polarizabilities and electron smearing factor (Thole) values were assumed from the given Drude prepper values. The Thole smearing factor is introduced to modify the electrostatic term of the non-bonded potential, in order to allow screened dipole interactions by smearing of the charge using a Slater distribution [36]. The water molecules assumed the same parametrization as used by Lamoureux et al. [28]. The graphene sheets were approximated by us, assuming the aromatic bonds of the graphene to be of the same kind as in cyclic amino acids, such as phenylalanine. The polarizabilities and the thole factor of the carbon atoms in graphene were adopted from the known values for phenylalanine, meaning that the polarizability was set to be −1.66 C 2 /N and a thole smearing factor of 1.30 (dimensionless). The partial charges were set to ±0.4191e with the positive charges corresponding to the real particles, while the negative ones were set for the Drude particles. The box size, system temperature, energy cut-off, water box and pressure were all taken identical to those values in the classical MD simulations. The simulation time for PMD simulations was set to 10 ns. The simulation protocol for the PMD simulations followed that of the classical MD simulations, starting with harmonically restraining the whole system, but water, and then a production run with the amino acids released. The amino acids binding energy was calculated by using equation (1).

Quantum chemical simulations
MD simulations allow essentially for the classical description of a molecular system and permit studying its dynamics on the millisecond time scale [37,38]. On the other hand, QC simulations are important to account for specific polarization effects that are otherwise omitted in the classical MD description and could only be accounted for phenomenologically, by e.g. PMD simulations. The main problem with QC calculations is the poor scalability of the computations, making it difficult to run simulations on systems that exceed 1000 atoms [39]. For the purpose of the QC calculations, the system was therefore reduced to a hexagonal graphene patch with a single amino acid on top as depicted in Figure 3. The QC simulations were carried out using the WB97XD/6-31G(d) method [40,41] employing Gaussian09 program [42], where single configuration geometry optimization was performed for each of the 20 naturally-occurring amino acids. The total energy of the system from these optimizations for each amino acid was obtained, which allowed computing the binding energy of the individual amino acids as where E total is the total energy of the system obtained from the optimization procedure, E aa is the internal energy of the given amino acid and E graphene is the internal energy of the graphene layer patch.

Results
The MD simulations reveal that different amino acids, such as e.g. arginine or tryptophan, exhibit strong binding configurations to the graphene surface while some other amino acids, such as e.g. glycine or valine, bind rather weakly. The established MD interaction energies are compared to those obtained from QC calculations, which are expected to be lower due to quantum effects, such as polarization. The degree of graphene polarization by amino acids is thus established for the 20 naturally-occurring amino acids. are presented for each amino acid. Those probability density distributions represent the sampled energies defined in equations (1) and (2) for the solvated and in vacuo simulations. Every subplot in these figures features a single amino acid type where the main peak corresponds to a binding configuration of an amino acid on graphene surfaces, while the increase in the interaction energy binding probability at around 0 kcal/mol indicates non-bound configurations. The energies corresponding to the peaks are compiled in Table 1.
Note, that for some amino acids several bound configurations were observed, which is reflected in the energies of multiple peaks in the interaction energy distribution. Each peak in the distribution should follow a Gaussian profile as dictated by the rules of statistical physics [43]. Multiple peaks that feature several bound configurations can, therefore, be described as a sum of independent Gaussian distributions given by: where E i is the mean interaction energy of an amino acid with graphene in a specific binding mode i, σ 2 i is the variance of the corresponding probability density distribution, A is the total number of observed bound configurations of an amino acid on top of graphene. The results of fitting the computed interaction energy distributions with equation (5) are shown in Figure 4 with thin grey lines. Figures 4 and 5 demonstrate that all amino acids showed an affinity for binding to graphene in aqueous solution and in vacuum. The strongest binding in water is seen for tryptophan while in vacuum it is for arginine. It is worth mentioning that earlier studies have suggested tryptophan to be the strongest binder to graphene in water [44], while other researchers have suggested arginine [45][46][47]. It is, however, conclusive [46] that the 3 best binders to graphene are arginine, tryptophan, and tyrosine, although the exact value of their adhesion is at present not consistent as the values reported in the present investigation differ from similar earlier studies [44][45][46][47]. Pandey et al. [44] suggests glycine being the weakest binder, which agrees well with the present results summarized in Table 1 and Figures 4 and 5.
The interaction energies of amino acid adhesion obtained from the PMD simulations have also been computed using equation (1) and the fitting of the corresponding probability density distributions was carried out as a sum of gaussian functions defined in equation (5). Note that such simulations were not feasible earlier as no Drude parameters for graphene were available in the literature. The results of the PMD simulations show overall the same trend as Figures 4 and 5, albeit the binding energy of the individual amino acids are overall more negative, and closer to the QC energy, discussed below. This lowering of the binding energy is due to the polarization forces taken into account semi empirically in the PMD simulations and shows that these forces are important for a system containing amino acids deposited on the graphene interface.
It could be suspected that amino acids that feature multiple peaks in the probability density distribution of the interaction energy (in the case of e.g. alanine and cysteine) are caused by the insufficient sampling time. Taking a closer look at the interaction of such amino acids with the graphene surface reveals this not to be the case. Figure 7 shows, for example, that alanine does not have a preferred binding configuration to graphene which leads to several possible binding energies as compiled in Table 1. Similar behavior is also observed for some other amino acids (asparagine, cysteine, and methionine).
The QC calculations were performed for all the 20 naturally occurring amino acids to reveal the degree of polarization that amino acids have on the graphene surface. Figure 8 depicts the Coulomb potential calculated for 5 selected amino acids and demonstrates that the charged amino acids give rise to a large Coulomb potential on graphene surfaces due to the polarization effect. Important to stress here is that the positively charged amino , equation (1), obtained for each amino acid deposited on graphene (colored curves). Each peak represents an energy associated with binding, while the increase in the probability density at 0 kcal/mol represents the non-bound configurations for the respective amino acid. The fittings (grey lines) were done using a sum of Gaussian distributions as introduced in equation (5).  (Fig. 8E) does not lead to any significant polarization, and, therefore, the binding, in this case, is dominated by the van der Waals forces.
The induced surface potential gives rise to a Coulomb interaction energy between an amino acid and the surface, which can be calculated as where q i is the charge of the ith atom in the amino acid while Q j is the charge of an atom j in the graphene surface. N aa is the total number of atoms in a given amino acid, while N g is the total number of atoms in the graphene sheet, r ij is the distance between atom i in an amino acid and atom j in the graphene sheet. The partial charges q i and Q j that follow the Mulliken and Electrostatic Potential Fitting (ESP) definition were obtained from the QC calculations and give rise to the Coulomb potential distribution in Figure 8. The values of all computed energies are provided in Table 1. Table 1 shows that amino acids' binding energies obtained from the QC calculations are generally lower than those from MD simulations. This effect is partially ascribed to the additional Coulomb interactions gained from the shift of charge distribution of the graphene layer. Furthermore the change in energy can be ascribed to QC effects such as introduction of the exchange-correlation term and localization of wavefunctions on the graphene surface around the amino acid binding site. The latter is especially observed for the charged amino acids. The shift in the case of the charged amino acids (arginine, aspartic acid, glutamic acid, and lysine) lowers the interaction energies by several times but appears to create a rather local induced potential on the graphene surface which is not affected by the terminal hydrogen atoms, as revealed in Figure 8. The induced charge distribution shown in Figure 8 is computed based on the calculated electronic wavefunction for the systems. Here the partial charges described to individual atoms could be slightly different, if computed following the Mulliken or ESP definition. Thus, the Coulomb energies computed using equation (6) for the ESP partial charges is overall lower than the energy obtained from Mulliken charges.
The binding energy of an amino acid to the graphene surface, E QC , has been introduced in equation (4) and it must include the different binding terms such that . Energies obtained for QC calculations EQC are also indicated. The values in the parenthesis provided in the, EQC, column lists the Coulomb energy contribution, see equation (6) to the amino acids' binding energies calculated from the Mulliken charges, while the ESP values [49][50][51] are given in the square brackets. The Coulomb energy contribution is arising due to the induced charges on the graphene surface. Note that the energies from MD simulations may have more than one expectation value due to several binding configurations being possible, while QC calculations only reveal one binding configuration. Uncertainty values correspond to the standard deviation obtained from fitting equation (5) to the simulation data. The PMD data for GLY and PRO are missing due to computational difficulties in converting the amino acids into the Drude approximation.   , E Coulomb and E QC are compiled in Table 1. The energy arising due to the polarization effects, E Polarization , can be estimated as the difference E , which turns out to be 4.99 kcal/mol for glutamine, as an example. This estimated difference corresponds to a polarization energy contribution in solution and is expected to be somewhat different in vacuo as water will impact the polarizational interactions for hydrophilic and hydrophobic amino acids. The binding energy found from the quantum chemical simulations can be taken to be representative of the value which should be closest to an experiment. However equation (7) indicates only a proportionality, as the energy contribution E , E Polarization and E Coulomb rely on the force field parameters that may not be exactly matching the resulting energies from the pure QC calculation. Furthermore, the limited basis set effects has contribution to some inconsistencies, making equation (7) an approximation.
From the performed investigation follows that the partial charge analysis in amino acids is dependent on the system rather than the calculation method employed. Fig. 8. Coulomb potentials at the graphene surface induced by five selected amino acids. Panels A and B show the Coulomb potential for the positively charged amino acids lysine and arginine, respectively. Panels C and D show the Coulomb potential for the negatively charged amino acids aspartic acid and glutamic acid, respectively. Panel E shows the Coulomb potential for the neutral amino acid tryptophan. The color code applies to all five panels. Analysis has been carried out using the program UCSF Chimera [48].
Earlier it was already speculated [51], that it is difficult to conclude whether the ESP or the Mulliken approaches deliver the most physically correct method for partial charge analysis. The presented calculations show that e.g. for LEU and ALA the Coulomb energy found using ESP charge analysis combined with the other two terms gives approximately the same value as the QC energy, whereas the same holds true for SER energy computed through the Mulliken approach. This reflects on the fact that it is somewhat difficult to quantify the contribution of Coulomb and van der Waals energies to the total binding energy, even from the rather accurate QC simulations. In general, the observed effect is reduced to a rather well-known problem [52][53][54] of a charged body interaction with a polarizable surface [55]. In some specific cases, such interaction can in principle be modeled accurately through the inclusion of image charges into the system or through polarizable force field models such as the PMD approach. One notes that the overall trend of the PMD data correlates with the QC results, while some discrepancies arise due to the model nature of the PMD approach. The observed discrepancies are also expected due to the fact that the QC simulations produce static optimized configurations of amino acids atop graphene, which are expected to have binding energy somewhat lower as compared to the MD configurations, which are dynamic. In MD simulations the configurations are not necessarily the most energetically favorable, which further explains the energy difference between E

Conclusion
The present investigation has utilized MD simulations to examine interactions of the 20 naturally occurring amino acids with layered pristine graphene in aqueous solution and in vacuo. It was found that some amino acids, such as e.g. arginine and tryptophan, exhibit strong binding to the graphene surface in both environments, while other amino acids, such as e.g. glycine and alanine, showed rather poor adhesion. QC calculations were further utilized in order to investigate whether quantum effects are capable to further enhance the interaction energies. These calculations revealed that a decrease in energy was systematically found for all 20 naturally occurring amino acids, while the decrease in binding energy was ascribed to the polarization interactions created by the amino acids on the graphene sheet, an effect that has also been observed through setting up and performing MD simulations with the polarizable force field. The analysis revealed that the resulting induced Coulomb interactions between amino acids and the graphene surface turned out to be several times larger than the van der Waals forces for induced Mulliken charges, however the ESP charges brought this interaction energy down indicating that polarization interactions could be more important for the interaction energy decrease since the van der Waals forces are the only type of interactions present in the classical MD approach. The binding energies found in MD simulations and the QC calculations are seen as crucial foundation for the further investigation that could demonstrate which amino acid would be most suitable to coat graphene surfaces in order to make it more compatible with proteins and other macromolecules, where fibronectin ( Fig. 1) is among the front line candidates for a follow-up investigation [56].
Open access funding provided by Projekt Deal. The authors would like to thank Peter Reinholdt for providing the geometry of the graphene patch used in the QC simulations. Financial support by the Lundbeck Foundation, the Danish Councils for Independent Research, and Volkswagen Foundation (Lichtenberg Professorship to IAS) and the Deutsche Forschungsgemeinschaft (GRK 1885) is greatly acknowledged. Computational resources for the simulations were provided by the DeiC National HPC Center, University of Southern Denmark.
Open Access This is an open access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.