Free-Energy Simulations of Hydrogen Bonding versus Stacking of

ARTICLE pubs.acs.org/JPCC

Free-Energy Simulations of Hydrogen Bonding versus Stacking of Nucleobases on a Graphene Surface
eza
c*,‡ Vojt
ech Spiwok,† Pavel Hobza,‡,§ and Jan R †

Department of Biochemistry and Microbiology, Institute of Chemical Technology, Prague, Technicka 3, 166 28 Prague 6, Czech Republic ‡ Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic and Center for Biomolecules and Complex Molecular Systems, Flemingovo nam. 2, 166 10, Prague 6, Czech Republic § Regional Centre of Advanced Technologies and Materials, Department of Physical Chemistry, Palacky University, 771 46 Olomouc, Czech Republic

bS Supporting Information ABSTRACT: It has been demonstrated by molecular modeling and experiments that free nucleic acid bases form hydrogen-bonded complexes in vacuum but prefer π
π stacking in partially and fully solvated systems. Here we show using molecular dynamics simulations and metadynamics that the addition of a surface (in this case a nanographene monolayer) reverts the situation from stacking back to hydrogen bonding. Watson
Crick as well as several non-Watson
Crick base pairs lying on a graphene surface are significantly more stable in a water environment than a π
π
π-stacked graphene
base
base assembly. It illustrates that the thermodynamics of nucleobase interactions results from a fine balance among hydrogen bonding, stacking, and solvation, and that these effects must be considered in molecular design.

’ INTRODUCTION The importance of hydrogen bonding in the stabilization of DNA duplex has been recognized simultaneously with the solution of the DNA structure.1 Much later, it was shown that the typical double-helical shape of DNA is also stabilized by inter- and intrastrand π
π stacking interactions.2,3 An interplay between these two types of noncovalent interactions is the key to understanding genetic and epigenetic information storage, DNA thermo- and photostability, its damage and repair, RNA secondary-structure formation, ribozyme function, and so on. For these reasons, isolated nucleic acid bases and their derivatives have been intensively studied in various contexts differing in their propensity for the formation of either hydrogen-bonded or stacked pairs. In the gas phase, the guanine
cytosine (G
C) Watson
Crick (WC) base pair featuring three strong H bonds is considerably more stable than the stacked configuration. In the adenine
thymine (A
T) pair, the difference between the H-bonded and stacked configuration is much smaller. This difference is further lowered in methylated analogs. Although the methyl-adenine
methyl-thymine (mA
mT) H-bonded pair is still lower in potential energy, transition to free energy tips the scales, as the stacked structure is preferred in terms of entropy.4,5 The situation changes dramatically when passing from gas phase to a solvent. A molecular dynamics simulation of the A
T pair in vacuum at 300 K revealed that hydrogen bonding is favored over stacking.2,6 The addition of a single water molecule does not r 2011 American Chemical Society

change the preference. However, hydrogen-bonded and stacked pairs become almost equivalently populated in the system containing two water molecules. Furthermore, microhydration by additional water molecules resulted in a significant preference for stacked structures. Hydrogen-bonded and stacked pairs of mA and mT are approximately equivalently populated in the absence of water, and addition of water molecules stabilizes stacked pairs. Similar results were obtained for G
C and mG
mC, where six water molecules change the preference from hydrogen-bonded to stacked pairs.5,6 The preference of nucleobases to form stacked complexes in water has also been verified experimentally.7,8 The above-described systems containing unsolvated, microsolvated, and solvated nucleobases can freely “choose” between either hydrogen-bonded or stacked complexes depending on the conditions. However, in DNA, both types of interaction operate simultaneously, and the final helical structure is determined mainly by the interplay between the interaction of bases under the constraints imposed by the backbone. Similarly, the movement of free bases can be limited by the presence of a surface. Here interaction with the surface competes with the interaction between the bases, and the adsorption of the bases on a hydophobic surface also affects the solvation of the bases. Received: March 16, 2011 Revised: August 30, 2011 Published: September 15, 2011 19455

dx.doi.org/10.1021/jp202491j | J. Phys. Chem. C 2011, 115, 19455–19462

The Journal of Physical Chemistry C

Figure 1. Example of one of the studied systems (graphene, mA and mT in the periodic box; water molecules not shown).

The smallest suitable model of such a system contains a base pair, water molecules, and a continuous sheet of benzene rings such as those in graphite or graphene (Figure 1). In this model system, bases can form different interbase pairings while both being stacked on the graphene/graphite surface. Alternatively, one base can lie on the graphene/graphite surface, whereas the second is stacked on the first one in a sandwich assembly. Understanding the factors affecting the equilibrium between these two forms is interesting from the point of view of the roles of noncovalent interactions in the stabilization of nucleic acids. This system is therefore a suitable extension in the series of systems studied in our laboratory, starting from the isolated base pairs in vacuum, through microsolvated and fully solvated base pairs to the situation in nucleic acids. Graphene, a monomolecular layer of graphite, is a “2010 Nobel Prize winning” material with unique electronic, optical, and other physical properties.9 The adsorption of nucleobases and nucleic acids and the subsequent formation of base pairs, self-assembled monolayers, and nucleic acid assemblies on graphene, graphite, and carbon nanotubes has been demonstrated experimentally10
17 and studied theoretically.18 The assembly of nucleic acid components on graphite and other water
solid interfaces is also interesting in the context of prebiotic evolution. Moreover, understanding these phenomena is also necessary for the design of future nanodevices exploiting the unique properties of graphene. Whereas most studies referenced above focus on the formation of an ordered 2D network of molecules adsorbed on the surface, we have focused on the interaction of a single base pair with the surface to understand the process of such self-assembly and its energetics and dynamics in detail. The binding of a pair of nucleobases on a graphene surface in a water environment is a very complex process. Its thermodynamics is determined by base
base, base
graphene, base
water, graphene
water, and water
water interactions. Some of these interactions are likely to favor the formation of hydrogen-bonded base pairs, whereas others disfavor them. A system containing microsolvated nucleic acid base pairs is already relatively complex; it is therefore not possible to rely on static structures and single-point energies, which means that dynamic calculations are required. However, the prediction of the equilibria from unbiased molecular dynamics trajectories

ARTICLE

can be successful only if the studied process is fast enough to be sampled on reasonable time scales. A proper sampling can be obtained in unbiased simulations of microhydrated base pairs, but a fully solvated system containing graphene is already too complex and requires the application of free energy modeling via biased simulation methods. Several methods have been developed to improve the sampling in molecular dynamics simulations or the Monte Carlo method. In metadynamics,19
21 the system is simulated by a molecular dynamics simulation with an additional bias potential. This potential accumulates during the simulation and energetically disfavors the previously visited configurations. This helps the system escape from the traps of the free energy minima. Moreover, there is a direct relationship between the free energy of a certain minimum and the bias potential that must be added to allow the system to escape this minimum. It must be kept in mind that metadynamics does not sample the studied system with an equilibrium probabilities. A free-energy surface is not calculated from the probabilities of different states along the simulation. Instead, the bias potential accumulated during the whole simulation (or more precisely its negative value) can be used as an estimate of the free-energy surface of the studied process. We studied free-energy surfaces of systems containing one molecule of nanographene (coronene 61, C150H30), two molecules of methylated nucleic acid bases (mA and mT, mG and mC), and 800 or 817 explicitly modeled water molecules. An example of the studied molecular system is illustrated in Figure 1. The metadynamics bias potential was acting to favor/disfavor the formation of stacked structure and WC base pairing.

’ COMPUTATIONAL DETAILS The metadynamics simulation was performed in Gromacs 4.0.722 with a Plumed 1.1.023 plug-in. The methyl derivatives of the nucleobases were modeled using the AMBER03 force field.24 Additional parameters of methyl groups were taken from the 5-methyl moiety of thymine. It is possible to model graphene as a continuous sheet by means of periodic boundary conditions, but using a finite sized nanographene has its advantages. We modeled graphene as a single molecule of nanographene (coronene 61) molecule to allow for a direct comparison with quantum mechanical (QM) calculations. It was modeled using generalized AMBER force field (GAFF)25,26 with the Gasteiger charges.27 The TIP3P model28 was used for water. The nanographene, both nucleobase derivatives, and water molecules were placed in a periodic 3  3  3 nm box (Figure 1). No molecules were restrained/constrained to their initial positions, and the box was used solely to confine water molecules. The system was geometry-optimized and then pre-equilibrated by one nanosecond unbiased molecular dynamics simulation. This was followed by a 100 ns metadynamics simulation, during which the hills of the bias potential were added every picosecond (100 000 hills in total). The bias potential was acting on two collective variables, namely, on stacking and WC hydrogen bonding. It was defined as coordination numbers Sstack ¼

Npurine Npyrimidine

∑i ∑j

1
ðrij =r0 Þ6 1
ðrij =r0 Þ12

for stacking and SWC ¼ 19456

Npurine Npyrimidine

∑i ∑j

1
ðrij =r0 Þ10 1
ðrij =r0 Þ12

dx.doi.org/10.1021/jp202491j |J. Phys. Chem. C 2011, 115, 19455–19462

The Journal of Physical Chemistry C

ARTICLE

Figure 2. Free-energy surfaces of the association of mA
mT and mG
mC pairs on graphene, calculated as the negative of the metadynamics bias potential. Minima A
N are highlighted. The axes are the collective variables describing base stacking and Watson
Crick pair formation, as described in the text.

Figure 3. Development of the Sstack (A), plane orientation (B), and free-energy difference between the stacked and WC structures (WC to stack free energy (C)) during the metadynamics simulation of the mA
mT pair association on graphene.

for hydrogen bonding. Sstack was calculated from distances rij between nine ring atoms of purine (i) and six ring atoms of pyrimidine (j). SW
C acts on donors/acceptors of WC hydrogen bonding (two pairs for mA
mT and three pairs for mG
mC). The reference distances r0 were 0.4 nm for Sstack and 0.3 nm for SW
C. The minima obtained from the metadynamics simulations of both systems were further studied using QM methods to verify the MM results independently. First of all, the geometry of the minima was optimized using the self-consistent-charges densityfunctional tight-binding29,30 with dispersion correction31 (SCCDFTB-D). This method has been shown to provide reliable relative values of the stabilization energies for H-bonded and stacked DNA base pairs,31,32 and our recent experience indicates that it can be successfully applied to graphene, where an efficient

method is needed to treat large systems. The free energy of these structures was estimated by evaluating the separate contributions at the best level practically applicable to systems of this size: The electronic energy is calculated using the density functional theory (DFT) with the TPSS functional in the TZVP basis set, supplemented by empirical dispersion correction.33 Both of these methods rely on the empirical dispersion correction, which makes it possible to handle such large systems efficiently with an accuracy comparable to the advanced correlated post Hartree
Fock (HF) methods. The solvation free energy of the aggregates is calculated using the solvation model D (SMD)34 based on the HF calculation in 6-31G* basis set. Finally, the thermodynamics at room temperature were evaluated using a rigid rotor/harmonic oscillator approximation based on a vibrational analysis at the 19457

dx.doi.org/10.1021/jp202491j |J. Phys. Chem. C 2011, 115, 19455–19462

The Journal of Physical Chemistry C

ARTICLE

Figure 4. Development of the Sstack (A), plane orientation (B), and then free-energy difference between the stacked and WC structures (WC to stack free energy (C)) during the metadynamics simulation of mG
mC pair association on graphene.

SCC-DFTB-D level. The SCC-DFTB-D calculations were carried out in the DFTB+ program,35 the DFT energies were calculated using Turbomole 5.10,36 and the solvation free energies were calculated using Gaussian 03.37

’ RESULTS The coordination number turned out to be a versatile collective variable for metadynamics studies because it efficiently recognizes the bound and dissociated states. A drawback of the coordination number is the fact that it approaches zero and hardly changes when molecules are fully dissociated (r . r0). Therefore, when both molecules are fully dissociated, their reassociation cannot be efficiently facilitated by a bias potential because a force acting in the space of collective variables becomes weak when converted to forces acting in Cartesian space on individual atoms. Therefore, the metadynamics simulation of dissociation/association by biasing a coordination number is efficient only for a partial dissociation/association, when the process is reasonably fast or in dissociation-only simulations. Fortunately, the graphene-assisted assembly of nucleobases is reasonably fast because the problem of association is converted from a 3D search in solution to a 2D search on a graphene surface. The addition of the bias potential during the metadynamics simulation promoted transitions between the state where both nucleobases were stuck on the nanographene as well as the graphene
base-1
base-2 π
π
π-stacked assembly (Figures 2, 3 and 4). Graphene
purine
pyrimidine as well as graphene
pyrimidine
purine assemblies have been observed, with the former being more common. Approximately 20 transitions between hydrogen-bonded pairs and π
π
π-stacked assemblies were observed during the 100 ns metadynamics simulations of mA
mT and mG
mC pairs. Besides this, also transitions between WC and non-WC base pairing were observed. Non-WC

base pairings often contain an H-bond of the WC pairing or at least a close proximity of WC H-bond donors/acceptors. As a result, the SWC as well as Sstack of non-WC base pairing can be distinguished from a fully dissociated form. Figure 2 shows the free-energy surface in the space of coordination numbers SWC and Sstack. It can clearly distinguish the dissociated, WC, and stacked forms. Structures corresponding to detected free-energy minima are depicted in Figure 5. A visual inspection of the trajectories revealed that all of the association processes took place in the 2D space on the graphene surface. Interestingly, the edges of the nanographene molecule did not cause any significant disturbance of the process of association. Several events of base flipping has been observed during both simulations. One of bases can spontaneously flip during a simulation by 180 degrees and stack again onto the graphene. As a result, the orientation of the bases changes from the cisorientation, suitable for WC pairing, to the trans-orientation, suitable for a reversed WC hydrogen bonding, or vice versa. This process has been monitored as the angle between the axes of both bases. The profile of this value (referred to as the plane orientation) is illustrated in Figures 3 and 4. Values of ∼0° correspond to the cis orientation, whereas values of ∼180° correspond to the trans orientation. Two different kinds of base-flipping transitions have been observed (Figure 6). The first pathway occurs on a graphene
base-1
base-2 π
π
πstacked complex. Base 2 is highly movable and can occasionally flip. In the second flipping pathway, the base flips directly on the graphene surface, and the second base acts as a “ski jump”. The process of base flipping complicates the calculation of the free energy profiles because it represents a slowly evolving degree of freedom that is not addressed by the metadynamics bias potential. The most straightforward way to eliminate or at least reduce this problem is to increase the length of the simulation to observe numerous flipping transitions. We observed 10 and 3 flipping 19458

dx.doi.org/10.1021/jp202491j |J. Phys. Chem. C 2011, 115, 19455–19462

The Journal of Physical Chemistry C

ARTICLE

Figure 5. Snapshots of the metadynamics simulations corresponding to minima A
N. The structures of the purines are superimposed within a group. These minima are Watson
Crick pairs (B,K), stacked pairs (F,N), dissociated bases (A,G), and other hydrogen-bonded arrangements (others).

Figure 6. Schematic view of two pathways of the base flipping between Watson
Crick and reverse Watson
Crick orientations.

events for mA
mT and mG
mC pairs, respectively. Nevertheless, base flipping remains the major source of error in the determination of the free-energy surface by metadynamics. Although the setup of the metadynamics favors conditions suitable for base flipping, namely, a formation of a π
π
πstacking, these events are rather rare and occur on a time scale of tens of nanoseconds. In principle, it is possible to include base-flipping as another collective variable in metadynamics, but this approach has its advantages and disadvantages. Besides (possible) improvement of sampling of base flipping, it can cause some problems such as release of bases from the graphene surface. The addition of another CV also reduces the numerical accuracy of metadynamics. Selected snapshots of simulations representing free-energy minima are depicted in Figure 5. It must be kept in mind that each free energy minimum calculated by metadynamics represents a family of configurations having similar values of CVs rather than a single unique configuration. Selection of snapshots of metadynamics with values of CVs close to the predicted minimum is therefore the most straightforward way how to structurally interpret the results of metadynamics. The purine base was superimposed to better illustrate the mutual orientation of both bases. Six minima were observed in the mA
mT pair. Besides the WC pairing (B), stacking (F), and separated complexes (A), the reversed WC pair (C) in the trans orientation is

also relatively stable. Another two structures (D,E) were characterized by interaction via the N7 of the purine as an acceptor (i.e., Hoogsteen and Hoogsteen-like pairing). Similarly to minima A and B, the structure D is also in cis orientation, whereas the structure E is in trans orientation. The situation was more complicated in the mG
mC pair. Besides the WC pair (K), stacked structure (N), and separated pair (G), there were five other energy minima. Some of them (H
J) were relatively heterogeneous. For example, minimum I was characterized by the presence of a single hydrogen bond with N1
H, N2
H, or both of mG and O2 of mC. The fact that there is only a single H-bond acceptor means that the mC may adopt different positions. These heterogeneous minima can usually be formed by both cis- and trans-orientated pairs. However, minima H
J and L are relatively shallow and are therefore not very important. The relatively stable minimum M is mostly composed of trans-oriented pairs similar to sugar edge pairs found in some complex nucleic acids. The quantum-mechanical SCC-DFTB-D optimization was applied to all of the minima obtained from the metadynamics shown in Figure 2; however, only some of them (A, B, and D
F for mA
mT and G and J
N for mG
mC) retained the original binding mode. Apparently, the partially dissociated structures are supported by water bridges missing in the gas-phase optimization. Because all of the important minima optimized well, 19459

dx.doi.org/10.1021/jp202491j |J. Phys. Chem. C 2011, 115, 19455–19462

The Journal of Physical Chemistry C

ARTICLE

Table 1. Relative Values of the Free Energies from QM Calculations (ΔGQM,total) and Their Components (electronic energy E; solvation free energy ΔGsolv; free energy of heating the system to room temperature ΔH0f300
TS), and from Metadynamics (ΔGmetadynamics, mean ( s.d.)a ΔGsolv ΔH0fT
TS ΔGQM,total ΔGmetadynamics

E

structure mA
mT A B C

0

0


58.5

35.0

8.8

0
14.7

0 6.8 ( 3.8 3.9 ( 3.7

D


64.8

37.7

19.1


8.0

4.9 ( 4.7

E


52.4

38.2

9.9


4.3

8.4 ( 4.0

F

33.3

22.5

9.7

65.5

0

0

mG
mC G

0

0

34.6 ( 3.4 0

H

2.8 ( 4.6

I

2.4 ( 3.0

J
63.0 K
129.1

a

0

43.5 78.2

11.8 24.5


7.7
26.4

1.1 ( 6.1
14.3 ( 2.2

L


32.9

33.1

11.9

12.1

M


69.0

35.9

22.4


10.7

8.3 ( 2.7 1.0 ( 19.2

N


3.4

32.3

11.9

40.8

27.2 ( 3.4

All of the values are in kilojoules per mole, relative to a dissociated state.

we have used only these (listed above) for comparison with results from the metadynamics. The separate terms of the QM free energy are listed in Table 1. Coordinates of complexes after geometry optimization are included in the Supporting Information. The results confirm that the absorption of both bases onto the surface is strongly preferred compared with the π
π
π arrangement. A detailed comparison of the multiple minima in a planar arrangement is not possible because of the inaccuracies in both metadynamics results (the estimated errors are listed in Table 1) and limitations of the approximations used in the QM calculations.

’ DISCUSSION Our results show a high preference of the nucleobases to bind on a graphene layer via stacking interactions. This is in agreement with the experiments10
17 showing that nucleobases tend to adsorb to polyaromatic surfaces, and it therefore implies a good performance of the molecular mechanics force field used in modeling of π
π stacking between a base and graphene/ graphite. It has already been shown that conventional empirical potentials provide a relatively accurate modeling of hydrogenbonding and stacking interactions between bases,32,38 despite the different physical nature of the types of interactions and relatively simple force-field formulation. This study indicates that π
π interactions can be accurately modeled by empirical potentials also in the context of nucleobase
graphene/graphite interactions. An interesting feature of graphene is the presence of low absolute atomic charges on the ring atoms when compared with benzene, small polyaromatics, or nucleobases. The charges on the carbon atoms tend to average toward zero with the increasing size of the (nano)graphene and with a distance from its edges (in this study formed by partially polar C
H bonds). In a hypothetical graphene with an infinite number of atoms, the atomic charges would be zero. Here we have used the easy-to-obtain Gasteiger charges for the nanographene molecule. The partial charges on the carbon atoms were ranging form
0.1261 to

+0.0414 with 64% of atoms in range of
0.01 to +0.01. This shows that electrostatics does not play an important role in nucleobase binding on graphene unless there is an important contribution from the induction (polarization) effects not modeled by a (nonpolarizable) force field. The studied system is well-suited for the application of a force field because the effect of polarization is expected to be negligible. The charge transfer that can be important in some graphene complexes39,40 is also very small because both graphene and nucleobases are electron donors. An interesting result is the existence of multiple minima for each pair. The collective variables used in this metadynamics calculation can distinguish not only separated, hydrogen-bonded, and stacked bases but also different WC and non-WC pairings. The existence of multiple minima in mG
mC showed the limit of the collective variables. For example, the minima H
J (Figure 5) are heterogeneous. Nevertheless, these minima are relatively unimportant and do not significantly change the overall picture of the pairing of mG
mC on graphene. Relative free energies of a dissociated state (A,G) might be influenced by peculiarities of modeling dissociation/association equilibria.41 Hydrogen-bonded WC pairing is relatively strong in terms of the stabilizing potential energy, especially in the mG
mC pair. Ab initio calculations have shown that WC pairing is stronger than stacking by 13.4 and 44.0 kJ/mol in mA
mT and mG
mC pairs, respectively, as calculated at the CCSD(T)/CBS level of theory in vacuum.2 This highly attractive nature of the hydrogen bonds causes that hydrogen-bonded complexes to be the vacuum free-energy minima of both pairs.5,6 However, the addition of water molecules can change the situation, as described in the Introduction, and stacked structures become more favorable in water solution. The formation of stacked complexes in solution has been also demonstrated experimentally. In this arrangement, their contact does not interfere with the favorable solvation of the polar edges of the molecules. Here we show that the addition of a single layer of graphene can change the preference form stacking back to hydrogen bonding, despite full solvation. The preference for the hydrogen-bonded complexes cannot be simply explained by hydrophobic effect because both types of assembly contain the same number of desolvated hydrophobic patches. In hydrogenbonded complexes, one patch on the nucleobase and one patch on the graphene layer per each nucleobase is desolvated. In the π
π
π stacking geometry it is one patch on graphene, two on the middle nucleobase and one on the top nucleobase. Both geometries are therefore comparable in terms of size of the desolvated hydrophobic surface (i.e., four patches for both geometries). The fact that purely hydrophobic explanation fails does not mean that the entropy does not play any important role. In the presence of the graphene surface, the interaction between bases competes with the graphene
base interactions. From the simulations, it is difficult to decide whether this effect is caused by the stabilization of the nucleobases in positions suitable for the formation of hydrogen-bonded structures, owing to the hydrophobic nature of the face of graphene, or the reduction of the configurational space of interacting molecules from 3D to 2D. The QM calculations allow a partial decomposition of the binding free energy. The largest term to determine the binding mode is enthalpy. The bases interact strongly with the graphene sheet18 (interaction energy in the range of 80
110 kJ/mol), and the stacking between the two bases2 is weaker (mA
mT, 55 kJ/mol; mC
mG, 75 kJ/mol). 19460

dx.doi.org/10.1021/jp202491j |J. Phys. Chem. C 2011, 115, 19455–19462

The Journal of Physical Chemistry C The graphene also partially shields the hydrogen bonds from solvent, which allows conservation of the hydrogen-bonded configurations even in solution. This suggests that stacking in DNA supports the formation of hydrogen bonds by the same mechanism; the power of such a shielding is apparent in the stability of the termini of double-helical oligomers. It is difficult to interpret the entropic contribution from the QM calculation because it does not include the effects of the solvent that restrain movement of the bases by friction or the formation of bridges between the bases. However, the magnitude of this term is not likely to make a difference in the overall ranking of the configurations. Our calculations show that the adsorption of individual bases is preferred as compared with the adsorption of a stacked base pairs by several tens of kilojoules per mole. This is in agreement with the experiments showing the formation of monomolecular layers of nuclobases on the surface. Our results represent one of many steps necessary for understanding the formation and application of assemblies containing (nano)graphene and nucleic acid components. Nucleobases or nucleic acids show a great potential for the solubilization, sorting, and delivery of carbonaceous nanomaterials. Different types of carbon nanotubes have been successfully sorted by specially tailored DNA molecules that wind around the nanotube.42
44 The bases of these single-stranded DNAs interact with the nanotube by stacking, and the bases of the two chains interact by hydrogen bonding. Another potential application of graphene
DNA interactions is the idea of single-molecule DNA sequencing by a controlled pulling of single-stranded DNA through a pore in graphene with electrochemical monitoring.45
47 Despite the fact that nucleic acids and carbonaceous nanomaterials can differ from the situation in our calculations (bases in nucleic acids are covalently and noncovalently linked, carbon nanotube surfaces are curved, etc.), they show that thermodynamics of these interactions must be properly addressed in the modeling of these assemblies. Our results also suggest a possible role of the surfaces in the chemical origins of life. They indicate that a surface could have formed a platform for a base-specific (and thus replicable) WC pairing, contrary to the base-nonspecific π
π stacking in a water environment and that inorganic surfaces have acted as “pairing catalysts” in prebiotic evolution.

’ CONCLUSIONS Metadynamics based on the molecular mechanics forcefield allowed us to obtain the free-energy surfaces describing the pairing of mA with mT and mG with mC on the surface of a graphene in water. In both cases, a hydrogen-bonded pair adsorbed on a graphene surface is a significantly more stable arrangement as compared with the adsorption of a stacked base pair onto the surface. These results were confirmed by a quantum-mechanical calculation of the most important contributions to the binding free energy. A comparison with the results on free, micro-, or fully solvated bases shows that these interactions are highly influenced by the local environment and that the addition of a solvent, surface, or both can change between hydrogen bonding and π
π stacking. These changes in the interaction type manifested at the level of free energy; therefore, free energy modeling (or, in principle, very long unbiased simulations) must be applied to study this phenomenon.

ARTICLE

’ ASSOCIATED CONTENT

bS

Supporting Information. Illustrations and coordinates of complexes studied by quantum chemical methods (after SCCDFTB-D geometry-optimization) are included. This material is available free of charge via the Internet at http://pubs.acs.org.

’ AUTHOR INFORMATION Corresponding Author

*E-mail: [email protected].

’ ACKNOWLEDGMENT This work was a part of research project no. Z40550506 of the Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic, and was supported by grant nos. LC512 (P.H.), MSM6198959216 (P.H.), and MSM6046137305 (V.S.) from the Ministry of Education, Youth and Sports of the Czech Republic. The work was also supported by the operational program Research and Development for Innovations of the European Social Fund (project CZ.1.05/1.1.00/03.0058). The support of Praemium Academiae, Academy of Sciences of the Czech Republic, awarded to P.H. in 2007 is also acknowledged. ’ REFERENCES (1) Watson, J. D.; Crick, F. H. C. Nature 1953, 171, 737. (2) Jure
cka, P.; Hobza, P. J. Am. Chem. Soc. 2003, 125, 15608.
erny , J.; Kabela
c, M.; Hobza, P. J. Am. Chem. Soc. 2008, (3) C 130, 16055. (4) Kabela
c, M.; Hobza, P. J. Phys. Chem. B 2001, 105, 5804. (5) Kabela
c, M.; Hobza, P. Chem.—Eur. J. 2001, 7, 2067.
eha, D.; Hobza, P. J. Phys. Chem. B (6) Kabela
c, M.; Zendlova, L.; R 2005, 109, 12206. (7) Nir, E.; Kleinermanns, K.; de Vries, M. S. Nature 2000, 408, 949. (8) Yanson, I. K.; Teplitsky, A. B.; Sukhodub, L. F. Biopolymers 1979, 18, 1149. (9) Geim, A. K. Science 2009, 324, 1530. (10) Heckl, W. M.; Smith, D. P. E.; Binnig, G.; Klagges, H.; H€ansch, T. W.; Maddocks, J. Proc. Natl. Acad. Sci. U.S.A. 1991, 88, 8003. (11) Srinivasan, R.; Gopalan, P. J. Phys. Chem. 1993, 97, 8770. (12) Tao, N. J.; Shi, Z. J. Phys. Chem. 1994, 98, 1464. (13) Freund, J.; Edelwirth, M.; Kr€ obel, P.; Heckl, W. M. Phys. Rev. B 1997, 55, 5394. (14) Sowerby, S. J.; Petersen, G. B. J. Electroanal. Chem. 1997, 433, 85. (15) Uchihashi, T.; Okada, T.; Sugawara, Y.; Yokoyama, K.; Morita, S. Phys. Rev. B 1999, 60, 8309. (16) Sowerby, S. J.; Stockwell, P. A.; Heckl, W. M.; Petersen, G. B. Origins Life Evol. Biospheres 2000, 30, 81. (17) Husale, S.; Sahoo, S.; Radenovic, A.; Traversi, F.; Annibale, P.; Kis, A. Langmuir 2010, 26, 18078. (18) Antony, J.; Grimme, G. Phys. Chem. Chem. Phys. 2008, 10, 2722. (19) Laio, A.; Parrinello, M. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562. (20) Iannuzzi, M.; Laio, A.; Parrinello, M. Phys. Rev. Lett. 2003, 90, 238302. (21) Laio, A.; Gervasio, F. L. Rep. Prog. Phys. 2008, 71, 126601. (22) van der Spoel, D.; Lindahl, E.; Hess, B.; Groenhof, G.; Mark, A. E.; Berendsen, H. J. J. Comput. Chem. 2005, 26, 1701. (23) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. Comput. Phys. Commun. 2009, 180, 1961. (24) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee., T. J. Comput. Chem. 2003, 24, 1999. 19461

dx.doi.org/10.1021/jp202491j |J. Phys. Chem. C 2011, 115, 19455–19462

The Journal of Physical Chemistry C (25) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollamn, P. A.; Case, D. A. J. Comput. Chem. 2004, 25, 1157. (26) Wang, J.; Wang, W.; Kollman, P. A.; Case, D. A. J. Mol. Graphics Modell. 2006, 25, 247. (27) Gasteiger, J.; Marsili, M. Tetrahedron 1980, 36, 3219. (28) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. J. Chem. Phys. 1983, 79, 926. (29) Porezag, D.; Frauenheim, T.; K€ohler, T.; Seifert, G.; Kaschner, R. Phys. Rev. B 1995, 51, 12947. (30) Seifert, G.; Porezag, D.; Frauenheim, T. Int. J. Quantum Chem. 1996, 58, 185. (31) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. J. Chem. Phys. 2001, 114, 5149. (32) Riley, K. E.; Pito
nak, M.; Jure
cka, P.; Hobza, P. Chem. Rev. 2010, 110, 5023–5063.
erny, J.; Jure
cka, P.; Hobza, P.; Valdes, H. J. Phys. Chem. A (33) C 2007, 111, 1146. (34) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. J. Phys. Chem. B 2009, 113, 6378. (35) Aradi, B.; Hourahine, B.; Frauenheim, Th. J. Phys. Chem. A 2007, 111, 5678. (36) Ahlrichs, R.; B€ar, M.; Baron, H.-P.; Bauernschmitt, R.; B€ocker, S.; Crawford, N.; Deglmann, P.; Ehrig, M.; Eichkorn, K.; Elliott, S.; Furche, F.; Haase, F.; H€aser, M.; H€attig, C.; Hellweg, A.; Horn, H.; Huber, C.; Huniar, U.; Kattannek, M.; K€ohn, A.; K€olmel, C.; Kollwitz, € M.; May, K.; Nava, P.; Ochsenfeld, C.; Ohm, H.; Patzelt, H.; Rappoport, D.; Rubner, O.; Sch€afer, A.; Schneider, U.; Sierka, M.; Treutler, O.; Unterreiner, B.; von Arnim, M.; Weigend, F.; Weis, P.; Weiss, H. Turbomole, version 5.9; Universit€at Karlsruhe: Karlsruhe, Germany, 2006. http://www.turbomole.com. (37) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Montgomery, J. A., Jr.; Vreven, T.; Kudin, K. N.; Burant, J. C.; Millam, J. M.; Iyengar, S. S.; Tomasi, J.; Barone, V.; Mennucci, B.; Cossi, M.; Scalmani, G.; Rega, N.; Petersson, G. A.; Nakatsuji, H.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Klene, M.; Li, X.; Knox, J. E.; Hratchian, H. P.; Cross, J. B.; Bakken, V.; Adamo, C.; Jaramillo, J.; Gomperts, R.; Stratmann, R. E.; Yazyev, O.; Austin, A. J.; Cammi, R.; Pomelli, C.; Ochterski, J.; Ayala, P. Y.; Morokuma, K.; Voth, G. A.; Salvador, P.; Dannenberg, J.; Zakrzewski, V. G.; Dapprich, S.; Daniels, A. D.; Strain, C.; Farkas, O.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J, B.; Ortiz, J. V.; Cui, Q.; Baboul, A. G.; Clifford, S.; Cioslowski, J.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Gonzalez, C.; Pople, J. A. Gaussian 03, revision C.02; Gaussian, Inc.: Wallingford, CT, 2004. (38) Hobza, P.; Kabela
c, M.;
Sponer, J.; Mejzlík, P.; Vondra
sek, J. J. Comput. Chem. 1997, 18, 1136. (39) Grimme, S.; Muck-Lichtenfeld, C.; Antony, J. J. Phys. Chem. C 2007, 111, 11199. (40) Kong, L. M.; Bjelkevig, C.; Gaddam, S.; Zhou, M.; Lee, Y. H.; Han, G. H.; Jeong, H. K.; Wu, N.; Zhang, Z. Z.; Xiao, J.; Dowben, P. A.; Kelber, J. A. J. Phys. Chem. C 2010, 114, 21618. (41) de Jong, D. H.; Sch€afer, L. V.; de Vries, A. H.; Marrink, S. J.; Berendsen, H. J.; Grubm€uller, H. J. Comput. Chem. 2011, 32, 1919. (42) Zheng, M.; Jagota, A.; Strano, M. S.; Santos, A. P.; Barone, P.; Chou, S. G.; Diner, B. A.; Dresselhaus, M. S.; McLean, R. S.; Onoa, G. B.; Samsonidze, G. G.; Semke, E. D.; Usrey, M.; Walls, D. J. Science 2003, 302, 1545. (43) Zheng, M.; Jagota, A.; Semke, E. D.; Diner, B. A.; McLean, R. S.; Lustig, S. R.; Richardson, R. E.; Tassi, N. G. Nat. Mater. 2003, 2, 338. (44) Martin, W.; Zhu, W.; Krilov, G. J. Phys. Chem. B 2008, 112, 16076. (45) Garaj, S.; Hubbard, W.; Reina, A.; Kong, J.; Barton, D.; Golovchenko, J. A. Nature 2010, 467, 190. (46) Schneider, G. F.; Kowalczyk, S. W.; Calado, W. E.; Pandraud, G.; Zandbergen, H. W.; Vandersypen, L. M. K.; Dekker, C. Nano Lett. 2010, 10, 3163.

ARTICLE

(47) Merchant, C. A.; Healy, K.; Wanunu, M.; Ray, V.; Peterman, N.; Bartel, J.; Fischbein, M. D.; Venta, K.; Luo, Z.; Johnson, A. T. C.; Drndi, M. Nano Lett. 2010, 10, 2915.

19462

dx.doi.org/10.1021/jp202491j |J. Phys. Chem. C 2011, 115, 19455–19462