Peptide interactions with zigzag edges in graphene

Graphene is a one-atom-thick planar sheet of carbon atoms densely packed into a honeycomb lattice. Graphene has four types of experimentally confirmed edges called zigzag, armchair, reconstructed zigzag, and extended Klein edges.^1,2 Below 400 °C, zigzag edges are dominant.³ The superior electronic, thermal, and mechanical properties of graphene have attracted significant attention in the development of high performance nanoelectronics and other applications.^4–10 It has been demonstrated that various graphene properties, including ferromagnetism, superconductivity, and an anomalous quantum Hall effect, are edge-dependent.^11–14 Furthermore, edge functionalization by selectively attaching chemical moieties at the edge of graphene sheets has the advantage of minimal damage of the carbon basal plane, and thus largely retaining the physicochemical properties of the pristine graphene. Therefore, the recognition and functionalization of graphene edges have significant impact on its various properties and applications.^15–18 

Peptides have been extensively explored to noncovalently functionalize single-wall carbon nanotubes and graphene.^19–22 To enhance the selectivity of graphene-based sensors using peptides as biorecognition elements, a graphene binding peptide [GBP: -Glu¹-Pro²-Leu³-Gln⁴-Leu⁵-Lys⁶-Met⁷-] has been recently identified from a combinatorial phage display library and found to preferentially bind to graphene edges without adversely affecting the properties of graphene.^19,23 It was shown that the glutamic acid at the N-terminal plays an important role in recognizing graphene edges through electrostatic interactions.²³ To further investigate the role of the glutamic acid, we have conceived two new peptides by replacing the glutamic acid (Glu) with glycine (Gly) (GBP-gly: -Gly¹-Pro²-Leu³-Gln⁴-Leu⁵-Lys⁶-Met⁷-) and randomly shuffling the whole amino acid sequence (GBP-x: -Gln¹-Leu²-Pro³-Met⁴-Glu⁵-Lys⁶-Leu⁷-). First, atomic force microscopy (AFM) image techniques were used to examine the assembly behavior of the peptides on graphene sheets. Then, extensive all-atom molecular dynamics simulations (MD) were carried out to understand the binding mechanism by constructing the potential of mean force (PMF) along interested reaction coordinates. A well-known challenge in free energy calculations using MD is to enable complete Boltzmann sampling of reaction coordinates.²⁴ Traditional unbiased MD simulations have difficulty bringing a system out of a local free energy minimum and over a high energy barrier to explore other minima and transition states. The binding of a peptide to a surface is even worse and prevents the molecular system from being adequately sampled within limited simulation time. To overcome the sampling challenge, several effective algorithms have been recently developed.^25–31 The biased-energy replica exchange MD method (biased REMD) proposed by Wang et al. is able to predict peptide adsorption free energy comparable to the corresponding experimental measurements.^25–27 The computational expense of the temperature-based REMD (T-REMD) can be greatly reduced by those schemes based on replica exchange with solvent tempering.²⁸ The steered molecular dynamics simulations (SMD) and metadynamics simulations are also proved to be powerful in achieving good estimates of the peptide-surface adsorption free energy.^29–31 Another appealing algorithm is the adaptive biasing force (ABF) method.^32–35 In ABF, the average force $F$ is estimated by sampling the instantaneous force after a threshold simulation. Then, the adaptive biasing force $−F$ is applied to the system on the fly to enable a uniform sampling along the reaction coordinates. An advantage of the ABF method over others is that it requires no prior unknown free energy landscape related parameter tuning such as the bias potential tuning in the biased REMD and the pulling velocity tuning in SMD. The effectiveness of the ABF method has been demonstrated in tackling several challenging problems in protein simulations.^36–39 Here, we demonstrate that the ABF method is powerful in calculating the free energy of peptide adsorption on graphene sheets. The computational results validate the experimental observations and provide insights into the binding mechanism.

Freshly cleaved graphene layers were prepared from highly ordered pyrolytic graphite [(HOPG) Ted Pella, Redding, CA] using the Scotch^™ tape method. AFM was used to ensure that each fresh graphene face was atomically smooth and debris-free. Typical graphene layer-to-layer edge height is ≤3 graphene layers (<1 nm thickness). The HOPG samples were then immersed in a peptide solution (0.2 mg/ml) for 15 min. Thorough rinsing with deionized (DI) water was followed by air-drying (50% relative humidity) for 2 h. AFM topographic images were taken using a Digital Instruments Nanoscope IV-multimode AFM in the noncontact mode. Peptides were purchased from Genescript, Inc. (99% purity).

AFM topographies were obtained using noncontact mode imaging with a 40 N/m spring constant, 300 kHz resonance frequency, and Al-coated asymmetric tip. All AFM images were obtained under 50% relative humidity lab air. AFM under aqueous conditions diminish the appearance of peptide topographic features due to the viscous damping due to water, and samples were imaged after removing unbound peptides and dried prior to imaging.

Since a minimum of three layers of graphene is necessary to avoid underlying substrate and water effect on peptide adsorption as demonstrated by Kim et al.,⁴⁰ a triple-layer graphite model was constructed using Materials Studio 8.0 (©2014 Dassault Systèmes) of dimensions 6.395 × 9.798 nm². When periodic conditions are imposed in MD simulation, they represent infinite graphene sheets as shown in Fig. 1. A smaller sized graphene (6.395 × 6.248 nm²) was placed above the triple-layer graphite. The top graphene layer periodically extends in the armchair direction (X direction) leaving the zigzag edge terminated with hydrogen. The C–C bond length 0.142 nm, the C–C–C bond angle 120°, and the interlayer equilibration distance 0.34 nm between graphene sheets are determined from experimental values.⁴¹ The system setup is shown in Fig. 1.

The i-tasser program⁴² (http://zhanglab.ccmb.med.umich.edu/I-TASSER/) is used for protein structure prediction for lengths of 10–1500 amino acid residues.⁴³ To use i-tasser for the initial structure prediction of 7-mer peptides, we appended one glycine residue at the N-terminus and two glycine residues at the C-terminus to meet the minimal requirement of ten residues. The predicted top models were used as the initial peptide models after deleting the appended glycine residues. The N- and C-terminus was acetylated and amidated, respectively, by patching the acetylated N-terminus residue ACE and amidated C-terminus residue CT2 in the Chemistry at HARvard Macromolecular Mechanics (CHARMM, https://www.charmm.org/) convention to avoid strong terminus–terminus electrostatic interactions. While peptide N- and C-terminus in the experiments were not acetylated and amidated, we capped the termini to exclusively study the sequence effect on binding and avoid the terminus–terminus (intrapeptide) interactions.

Currently, both polarizable and traditional additive force fields are used to describe graphene.^44–50 Since the CHARMM additive force fields yield good results in protein adsorption on graphene^40,49 and graphene water contact angles,⁴⁸ we have generated all the initial model systems using visual molecular dynamics (Ref. 51) and CHARMM27 force field.^52,53 All molecular dynamics simulations were carried out using namd software package⁵⁴ with CHARMM27 force field. The bonded and unbounded parameters for graphite were taken from CHARMM atom type CA as in our previous studies.⁴⁰ Electrostatic partial charges 0.115 e and −0.115 e were, respectively, assigned to the atom H and C on the zigzag edge as imported from the CA-HP group in the CHARMM residue PHE benzene ring, while other carbon atoms in graphene are neutral.

The adsorption of peptides on graphene surface is a very slow process in water.⁴⁵ A crude procedure was first employed to assemble the starting simulation systems, and then extensive molecular dynamics simulations were carried out to equilibrate the systems. First, a peptide was placed 1 nm (center of peptide mass to plane distance) above the top graphene sheet. REMD simulations⁵⁵ were carried out for 10 ns using eight replicas in vacuum with temperature ranging from 300 to 600 K to enhance peptide–graphene interactions. The cutoff, switching, and pair list distance were chosen as 2.4, 2.0, and 2.8 nm, respectively. The configuration with the lowest potential energy at 300 K was selected. A pre-equilibrated TIP3P water box of 4486 water molecules with dimension of 6.395 × 9.798 × 2.5 nm³ was added above the peptide-graphite systems. Any overlapping water molecules with the peptide were randomly relocated on the top surface of the water box. The system was first minimized using a sophisticated conjugate gradient and line search algorithm as described in the NAMD user's guide. Then, the system was gradually heated up to 300 from 0 K by increasing the temperature to 30 K for every 1000 steps. Finally, the system was equilibrated at a constant temperature of 300 K and a constant pressure of 1 atm (NPT ensemble) for 10 ns. The obtained end-point configuration was served as the starting configuration for further simulations. In all the simulations in water, periodic boundary conditions and a time step of 1 fs were adopted in all simulations. The particle mesh Ewald method was used to calculate the long-range electrostatic interactions. The cut-off, switching, and pair-list distance were chosen as 1.2, 1.0, and 1.4 nm, respectively.

Our goal is to identify the most favorable binding region and binding residue of the peptide on graphene. The goal can be achieved by calculating the potential of mean force using the ABF algorithm when the peptide is assumed to move from one location to another on the graphene surface. Therefore, the projection $DY$ on axis Y of a distance vector $D$ is an appropriate reaction coordinate to measure the movement, where $D$ is the distance vector defined from the center of the top graphene layer to the mass center of alpha carbons of the peptide. To enhance the sampling, we have divided the region of interest into nine equally spaced windows of width 4 Å. The window centers are at −36, −32, −28, −24, −20, −16, −12, −8, and −4 Å. To obtain smooth free energy profile, we have divided each window into 40 bins so that the bin size $δDY$ is 0.1 Å. Therefore, the free energy difference can be written as³⁵ 

where $FDY$ is an instance force acting on the reaction coordinate in the ith bin.

To implement the formulae (1), we have first positioned the peptide inside each window using a biased harmonic potential with an actual force constant 30 kcal/mol/Ų. Then, half-harmonic potentials imposed at the boundary of each window with sufficiently strong harmonic restraints 100 kcal/mol/Ų were used to constrain the mass center of the peptide inside each window. A threshold of 1 000 000 force samples (full samples parameter per window) was used to obtain the estimate of the force distribution prior to applying the adaptive biasing force. ABF simulations were conducted for 60 ns in each window so that the total samples are 60 000 000 per window and the total simulation time is 540 ns for each system. The convergence is checked by plotting sequential PMFs during the simulation.

Graphene sheets obtained from freshly cleaved HOPG were exposed to GBP, GBP-gly, or GBP-x peptide solution to allow for binding to the graphene surface. The peptide-coated graphene sheets were then analyzed by AFM in the noncontact mode. Figure 2 shows the topographical images and height histograms of the three different peptides adsorbed onto graphene surfaces.

The peptides exhibit significantly different assembly behavior on the graphene/graphite surface. The parent peptide GBP dominantly forms aggregates at graphene edges. The height of these peptide aggregates is ∼1 nm, suggesting that monomer peptide units indeed accumulate at the edges and aggregate to form short fibrillary nanostructures. Substitution of the N-terminal Glu to Gly (GBP-gly) results in a different binding behavior and formation of distinct globularlike structures across the graphene surface (height = ∼2 nm). In contrast, when the parent peptide is scrambled (GBP-x), a reticular network of peptide structures are formed on the graphene surface. These formed structures are in height (∼1 nm) similar to those of GBP but stretch across the entire planar graphene surface. Taken together, the results suggest that the glutamic acid residue and its position in the context of the peptide sequence are critical for its graphene edge binding property. Next, we resort to molecular dynamics simulation to understand the mechanism for the assembly behavior.

It is demonstrated that zigzag edges are dominant below 400 °C.³ The high-resolution scanning tunneling microscopy and first-principles calculations have shown that the graphene edges are terminated by hydrogen atoms with no rehybridization of the outermost carbon edge atoms.⁵⁶ Therefore, here we investigate the effect of zigzag edges terminated by hydrogen atoms on the adsorption of the peptides. When we take advantage of the periodicity to construct the infinite graphene model from the simulation unit cell as shown in Fig. 1, we have carefully determined the right size of the simulation unit cell. If the size is too big, increased computational resources will be required. On the other hand, if the size is too small, artifacts due to the interaction between a molecule and its image may occur. Since we want to compare the binding difference of our peptides to graphene with and without edge effect, it is necessary to have about 20 Å margin from the edge to avoid the edge effect when a peptide is at the center of the unit cell. Since the maximum length of a 7-mer peptide is about 24.5 Å, we constructed a stripe of width 62.48 Å with zigzag edges. The infinite stripe with zigzag edges is put on the middle of trilayer graphene to model the edge effect. The width of the unit cell of the trilayer graphene is 97.98 Å with a margin 35.5 Å to the zigzag edges to avoid edge–edge interactions. The setting of the model system is shown in Fig. 1.

By construction, the graphite model possesses two different regions for us to investigate the edge effect on the peptide adsorption—at the edge and on basal plane of graphene. The PMF may reveal the binding preference of a peptide on different locations. The adaptive biasing forces method is employed to determine the PMF profile as described in Sec. III. We have chosen the projection $DY$ of the distance vector $D$ on the axis Y as the reaction coordinates, as explained in Secs. III A and III D. The distance vector $D$ is defined from the center of the top graphene layer to the mass center of alpha carbons in a peptide. Due to the symmetry of the simulation system, we need to consider only the range of $DY$ from −38 to −2 Å. For effective sampling, the range was divided into nine equally spaced windows of width 4 Å. In each window, 60 ns ABF molecular dynamics simulations at NPT ensemble and thus total 540 ns ABF simulations for each system were carried out. Since a continuously updated biasing force is added to flatten the free-energy surface as the simulation progresses, a broad sampling is achieved in 60 ns as shown in Fig. 3(a). During our simulation, we monitored the convergence by plotting the sequential PMFs every 10 ns and found out the PMFs converge after 50 ns, as shown in Fig. 3(b). The PMF profile in Fig. 3(b) shows that the peptide has a uniform binding affinity on the graphene plane with D_Y ranging from −24 to −2 Å, but the stronger binding affinity appears near the edge (⁠$DY ≈−36 Å$⁠). The energy barrier from the graphene planar to the edge region is less than 3 kcal/mol, but the energy barrier from the edge to the planar region is >6 kcal/mol. Therefore, the preferable binding site of GBP is on the edge where the binding free energy is ∼4 kcal/mol lower than the plane binding free energy. The modeling results are in agreement with AFM observation shown in Fig. 2(a).

To address the reliability of the above individual free-energy calculation $ΔA⌢$⁠, ideally we should repeat the calculations N times by varying starting atomic velocities or positions and then obtain the average $E(ΔA⌢)$ to see whether it is comparable to the experimental value $ΔA$ by measuring the mean-square error as³⁵ 

where the first term is the measure of the statistical error, and the second term is the measure of the accuracy in comparison with the experimental value. The estimation of the first term requires many independent simulations which are very computationally expensive since each needs 540 ns simulation. The estimation of the second term requires experimental data, which is not available. Therefore, the characterization of errors associated with free-energy calculations is a very challenging problem.^35,57 In recent years, progress has been made to estimate the statistical error.^35,57 Within the adaptive biasing force free energy calculation framework based on Eq. (1), we can directly estimate the statistical error by computing the variance as

A MatLab code is written to calculate $Var(ΔA)$ using Eq. (3), which gives $Var(ΔA)=0.026256$⁠. Thus, the statistical error of the free energy calculation in our case is estimated as $Var(ΔA)=0.16$ kcal/mol.

From the end-point conformation corresponding to the lowest free energy near the edge (⁠$DY ≈−36 Å$⁠) shown in Fig. 3(c), it is clearly seen that the carboxylate group of the N-terminus Glu residue interacts with the zigzag edges. The important role of the Glu amino acid residue binding to graphene edges through electrostatic interaction can be verified by varying the solution pH value. The side chain carboxylic acid functional group has a pKa of 4.1, and thus at neutral pH, it exists in its negatively charged deprotonated carboxylate form. At pH 3, the carboxylate ion is neutralized so that it is expected to weaken the binding to graphene edges. In our earlier studies, we demonstrated that the protonation of the glutamic acid disrupts graphene edge binding.²³ We have recalculated the potential of mean force assuming that the glutamic acid is protonated (the peptide is named GBP-h). From the obtained PMF shown in Fig. 4(a), it is seen that the global free energy minimum has moved toward the graphene plane (⁠$DY≈−36→−28 Å$⁠), whereas the global minimum at −36 Å for GBP has become the global maximum for GBP-h. The difference between the global minimum free energy (⁠$DY≈−28 Å$⁠) and the free energy for GBP-h positioned on plane (⁠$DY≈−24∼−2 Å$⁠) is <1 kcal/mol, whereas the energy barrier is <2 kcal/mol. This implies that the edge-preferred binding property vanishes after the glutamic acid is protonated, which is in agreement with our previous experimental observation.²³ 

In experiments, the variation of pH may also change graphene's properties. To make sure that the loss of peptide edge-preferred binding property is solely due to the protonation state change of the glutamic acid, we have replaced glutamate reside with glycine (GBP-gly) to further test the role of the glutamate residue in neutral pH solution. The AFM topographic image in Fig. 2(b) shows that the peptide GBP-gly appears distributed across the graphene surface. The shape of the recalculated free energy profile shown in Fig. 4(b) for GBP-gly is similar to that for GBP shown in Fig. 3(b). However, the free energy difference between binding to the edge and to the plane is ∼1 kcal/mol. This implies that the edge-preferred binding indeed disappears after the glutamic acid is modified to glycine. The end-point conformation corresponding to the lowest potential energy for GBP-h and GBP-gly is shown in Figs. 4(d) and 4(e), respectively. It is seen that the protonated glutamate and glycine interact with graphene plane and is away from the edge. The study of both GBP-h and GBP-gly reveals the significant role of the carboxylic group of glutamate in interaction with the graphene edges.

One of the main questions that has intrigued researchers studying peptide-surface binding using peptides isolated from phage peptide display libraries is whether amino acid positional context is important for binding. To this end, we have randomly shuffled the amino acid order to produce the GBP-x peptide. The AFM topographic image shown in Fig. 2(c) of GBP-x on graphene/graphite exhibits no preference for the graphene edge. The PMF profile shown in Fig. 4(c) implies that the peptide GBP-x prefers to bind to the graphene plane (⁠$DY≈−24∼−2 Å$⁠). From the end-point conformation shown in Fig. 4(f), it is seen that the carboxylic group sticks out to water, and so does the side chain of its neighbor Lys residue. A possible salt bridge between them may form when the glutamate side chain is deprotonated, whereas a hydrogen bond may form when the side chain is protonated. These two hydrophilic residues together pull the peptide toward water, thereby weakening its interaction with graphene. This explains the broad distribution seen in AFM height histogram with a broader range shown in Fig. 2(c).

The physical properties of graphene-based devices are influenced by graphene edges. To make controllable biofunctionalized graphene-based biosensors, it is necessary to understand biotic–abiotic interactions. An edge-binding peptide has recently been identified from phage display experiments. It is suggested that the glutamic acid at the N-terminus plays the key role in the edge-recognition. We have further investigated the role of the glutamic acid by substituting it with glycine, randomly shuffling the amino acid order and changing its protonation state. In each case, different peptide adsorption features have been observed by atomic force spectroscopy. It is found that if the glutamic acid residue is replaced with glycine or shifted to another position within the peptide sequence, the specific edge-binding ability of the peptide is lost. Free energy calculations were carried out to explain the differences in binding mechanism using the ABF method. It is demonstrated the obtained free energy profiles and end-point conformations are very useful in interpreting the experimental data and providing insights into the different binding mechanism. This study is useful in designing peptides to recognize and functionalize graphene edges.

The authors thank J. C. Gumbart and T. R. Walsh for providing valuable technical comments. The authors are also thankful to Rajiv Berry and Frank Duan for the DoD supercomputer support. This work was funded by AFOSR (R.R.N.). M.C.M. also acknowledges the support from the Air Force Office of Scientific Research (No. FA9550-12-1-0367).