On the Association of the Base Pairs on the Silica Surface Based on

Article pubs.acs.org/JPCC

On the Association of the Base Pairs on the Silica Surface Based on Free Energy Biased Molecular Dynamics Simulation and Quantum Mechanical Calculations Susanta Haldar,†,‡ Vojtech Spiwok,§ and Pavel Hobza*,†,∥ †

Institute of Organic Chemistry and Biochemistry and Gilead Science Research Center, Academy of Sciences of the Czech Republic, Flemingovo nam. 2, 166 10 Prague 6, Czech Republic ‡ Department of Physical and Macromolecular Chemistry, Faculty of Science, Charles University in Prague, Albertov 6, 128 43 Prague 2, Czech Republic § Department of Biochemistry and Microbiology, Institute of Chemical Technology, Prague, Technicka 3, 166 28 Prague 6, Czech Republic ∥ Regional Center of Advanced Technologies and Materials, Department of Physical Chemistry, Palacký University, Olomouc, 771 46 Olomouc, Czech Republic ABSTRACT: The adsorption of the DNA bases and base pairs on the hydrophobic silica surface has been investigated by ab initio quantum mechanical (QM) methods (DFT−D) and molecular mechanics (MM) and also by biased molecular dynamics (MD) simulations (metadynamics). The structures of all the clusters (surface with single-bases and base pairs) predicted by means of the force field are compared with the results of direct QM calculations. The MM interaction energies for all clusters agreed well with the QM ones, which justifies the use of MM methods in the evaluation of accurate adsorption free energies. Rigid rotor−harmonic oscillator− ideal gas (RR−HO−IG) calculations based on QM and MM entropies as well as biased metadynamics (MTD) simulations based on MM demonstrated that mA−mT (Adenine−Thymine) and mG−mC (Guanine−Cytosine) base pairs are adsorbed on a fully solvated silica surface in different H-bonded forms. Both QM and MM techniques convincingly demonstrated that adsorptions of any H-bonded (Watson−Crick (WC), non−WC, and Hoogsteen) structures on the surface are stronger than that of any π−π stacked structures.

1. INTRODUCTION Hydrogen bond plays an important role in the field of biological chemistry. They are known to determine the structure (and thus also the function) of biomolecules. Watson and Crick1 explained the stability of the DNA duplex by the presence of hydrogen bonding (H-bonding) between purine and pyrimidine nucleobases, where the guanine (G) and cytosine (C) pair has three and the adenine (A) and thymine (T) pair has two hydrogen bonds. Later, it was shown that besides hydrogen bonding, stacking also (intra as well as interstrand) contributes to the stability of the DNA duplex.2−4 In our previous studies,5−7 we have shown that planar Hbonded structures of A−T and G−C pairs, as well as of their methylated analogues, are more stable than the T-shaped and stacked structures. Passing to the water environment with several water molecules (microsolvation) as well as bulk water, the scenario is changed and the stacked structures of mA−mT and mG−mC base pairs become the most stable structure. The question concerning the preferential structure of a base pair now remains when it is adsorbed on the solid surface. This question is of immense importance, because it can elucidate the © XXXX American Chemical Society

role of noncovalent interactions in the stabilization of nucleic acids and also in the stabilization of genetic codes. Structure and function of nucleic acids are determined by an interplay between hydrogen bonding and stacking interaction of their bases. The former type of interaction is responsible for the storage and transfer of genetic information, whereas the latter stabilizes nucleic acids and plays other roles. An important question is how abiotic factors in early evolution, such as mineral surfaces, enforced formation of one of these types of noncovalent forces.8,9 In our previous article,10 we investigated the adsorption of methylated DNA base pairs on a nonperiodic graphene surface. A classical Molecular Dynamics (MD) simulation was performed for a time scale of 1 ns in explicit water, where both base pairs showed the stacked assembly. When the system was exposed to an external bias potential, it was clearly shown that both base pairs prefer to stay in various H-bonded forms rather than in the stacked form. Received: January 7, 2013 Revised: April 17, 2013

A

dx.doi.org/10.1021/jp400198h | J. Phys. Chem. C XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry C

Article

and was shown to be a suitable tool for the calculations like of the base pair adsorption on the silica surface in the previous papers.15,19−21 The aim of the present article is to investigate the adsorption of various base pairs on the silica surface. Here, we introduced a nonperiodic 3D silica surface containing 194 atoms (Si60O106H28), on which the methylated purine and pyrimidine bases were adsorbed. Figure 1 shows a model of the base pairs with a silica surface.

In this article, another type of surface was considered, namely, a nonrigid hydrophobic silica surface. Silica and silicates are among the most widespread constituents on the earth, and this surface is likely to participate in the origin of life. Silica has also been one of the most frequently studied materials in surface science. Adsorption processes of several heterocyclic compounds, isolated bases, and base pairs on silica and aluminosilicates (mainly clay minerals) are therefore relevant from the environmental and evolutionary point of view.11−13 Silica-based materials have an extremely wide range of applications, ranging from electronics to catalysis, e.g., peptide-bond formation.14 The flexibility of the Si−O−Si bridge between the two layers is responsible for the structural varieties of several silica solid surfaces. Due to the flexibility of the base pairs during a biased metadydynamics (MTD) simulation,15 all kinds of translational (pyrimidine over purine or vice versa), rotational (base flipping), and vibrational degrees of freedom have been considered in this paper. Similarly to our previous study,10 here we consider purine and pyrimidine bases methylated at N9 and N1 position, respectively. A consequence of methylation is the fact that the average potential energy barrier between several H-bonded and π−π stacked (pyrimidine over purine or vice versa) base pairs becomes lower and it also restricts the formation of H-bonds associated with N9−H.5 Therefore, the number of H-bonds formed becomes lower. Due to the flexibility of base pairs, the transition between several H-bonded and also between Hbonded and surface−π−π stacking configurations or vice versa takes place more frequently. The pair association in terms of WC and non−WC Hbonding and π−π stacking is a complicated process. Problems arise from the fact that the thermodynamics of all the relevant interactions (base/base, base/surface, base/water, surface/ water, and water/water) are determined by classical MD. The prediction of the equilibria between H-bonded and stacked structures from unbiased molecular dynamics trajectories can be successful only if the studied process is fast enough to be sampled on reasonable time scales. Therefore, a fully solvated system containing any base pairs on an organic or inorganic surface is already too complex. The sampling hence requires the application of more advanced free energy modeling which is fast enough to construct the free energy profile of the base movement. A suitable method for such a study is the biased molecular dynamics simulation technique. Up to now, several methods have been developed to bias a classical molecular dynamics simulation to enhance the sampling of the various minima on the free energy surface (FES). The metadynamics method is defined by several collective coordinates and the introduction of a historydependent biasing potential. The biased potential mainly prevents the system from visiting the regions that it has already explored before. This helps the system escape from the traps of the free energy minima. There are many circumstances where the free energy surface (FES) has several local minima separated by larger barriers, such as protein folding,16,17 firstorder phase transitions, or chemical reactions.18 In such cases, a simulation biased by a harmonic potential with collective variables starts to find each minimum on the FES and moves spontaneously to the next minimum only under very favorable circumstances. The base pair adsorption on the surface is a very complicated process because of the existence of several minima and fast flipping processes sometimes connected with a large barrier. MTD has successfully been used to explore the FES

Figure 1. Example of one of the systems studied, a snapshot of the mG−mC base pairs with the silica surface taken from the MTD simulation (the water molecules are not shown here).

2. COMPUTATIONAL METHODS In all QM and MTD cluster calculations, a finite 3D hydrophobic silica surface with the methylated base pairs has been used. The hydrophobic surface was taken from the literature.22 The methylated base pairs were prepared manually replacing the H-9 with a methyl group. The Lennard-Jones parameters of the surface were taken from the literature.23 In order to compare the molecular mechanics (MM) potential energy surface with the QM one, we investigated single-base (purine and pyrimidine) adsorption on the surface. To find suitable Lennard-Jones parameters for silicon and oxygen, all surface-single-base geometry was optimized at the MM level and the resulting MM interaction energy was compared with QM interaction energy. The parameters for silicon (Bondi’s values for van der Waals radii) are shown in Table 1. The interaction energies for all the surface-single-base complexes were determined at two levels of theory: dispersion corrected density functional theory (DFT−D)24,25 with TPSS− D26/TZVP and B3LYP−D/TZVP functionals, and MM force field. The gas phase interaction energy between subsystems R and T is defined as follows: ΔEgas(RT ) = Egas(RT ) − Egas(R ) − Egas(T )

(1)

where RT stands for whole complex in the gas phase. The empirical potential calculations were based on the General Amber Force Field (GAFF).27,28 For GAFF calculations, the atomic partial charges for surface and also base pairs were calculated according to the restricted electrostatic potential (RESP) approach.29 The RESP fit was performed onto a grid of electrostatic potential points calculated at the HF/6-31G* level as recommended by the force-field designers.28 The RESP charges calculated at the HF/6-31G* level were considered in the MM calculation in explicit water and the B3LYP/cc−pVTZ RESP30 charges were considered in the gas phase interaction energy calculation . B3LYP/cc−pVTZ DFT B



Article

Table 1. Comparison of the Silica Parameters Obtained from the QM and MM Methodsa Surface

TPSS−D/TZVP

Si−O bond length Si−H bond length Si−O−Si angle O−Si−O angle a

0.164 0.146 165.5 110.1

± ± ± ±

0.0001 0.001 13.2 2.1

B3LYP−D/TZVP 0.163 0.146 160.3 109.1

± ± ± ±

0.0004 0.001 20.5 0.3

0.164 0.146 168.8 109.0

MM

Parameters in MM

± ± ± ±

εSi = 1.06b σSi = 0.382c εO = 0.48 σO = 0.311 εH = 0.066 σH = 0.106 rSi = 0.210d rO = 0.145 rH = 0.122

0.0002 0.001 11.4 0.4

Bond length is in nm; angles are in deg. bkJ/mol. cnm. dBondi’s values for van der Waals radii in nm.

accurate relative stabilization energy.41,42 The basis set superposition error (BSSE) has been systematically included.43 The flipping of the bases at the surface from WC to non-WC and from π−π stack to T-shape configuration was calculated at TPSS−D/TZVP//COMSO44 level. The binding free energies were determined within the applied force field by using the rigid rotor−harmonic oscillator−ideal gas (RR−HO−IG) approximation.45 The size of the complexes investigated did not make it possible to use the more sophisticated QM and MM methods. A generalized Born46 implicit solvent model was used to calculate the solvation free energy. The implicit solvation parameter (van der Waals radii) for silicon atoms was taken from the literature47 and is shown in Table 1. The overall binding free energy of complex RT in vacuum (ΔGgas(RT)) is constructed as follows:

method is more accurate than the HF/6-31G* because the latter lacks electron correlation. Consequently also the RESP charges calculated with DFT are expected to be more accurate than the HF charges. This was confirmed in a study of peptides in the gas phase.30 The HF/6-31G* RESP charges are larger (more negative) than the B3LYP/cc−pVTZ ones.31 The DFT−D calculations were performed in the Turbomole 6.3 program package;32 the RESP charge and MM calculation were performed with the Gaussian0933 and AMBER packages,34 respectively. For the MTD simulation, the surfaces including base pairs were placed into a rectangular box of 3.5 × 2.6 × 2.4 nm3 along with 661 TIP3P35 explicit water molecules. No molecules were restrained or constrained to their initial positions and the box was used solely to confine water molecules. The geometry of the whole system was optimized and equilibrated by a 50 ns unbiased classical molecular dynamics simulation. All simulations were performed with a 1.0 fs time step at a temperature of 300 K and overall a constant pressure and temperature NPT ensemble was utilized with a Nose−Hoover36,37 thermostat in each case to couple the system to a heat reservoir. The classical dynamics simulation was followed by a 130 ns biased metadynamics simulation, during which the hills of the biased potential were added at every picosecond (130 000 hills in total). The bias potential was acting on two collective variables, namely, on stacking and WC hydrogen bonding. For the stacking (Sstack) and hydrogen bonding (SWC) of the base pairs on the silica surface, the following coordination numbers were defined as Npurine Npyrimidine

Sstack =

∑ ∑ i

j

ΔGgas(RT ) = Ggas(RT ) − Ggas(R ) − Ggas(T )

The free energy is defined in the following manner Ggas = Egas + Hgas(0 → T ) + ZPVE − TSgas

ΔGsol = ΔGgas + ΔΔGsol

Swc =

∑ ∑ i

j

(2)

3. RESULTS 3.1. Interaction Energy. Table 2 shows the DFT−D and MM interaction energies of the silica surface with single nucleic acid bases. The MM stabilization energies are larger on average by 8.53 and 7.40 kJ/mol than both the DFT−D (TPSS−D and B3LYP−D) energies, which are very similar. The DFT−D

1 − (rij/r0)10 1 − (rij/r0)12

(6)

ΔΔGsol represents the total solvation free energy, which is the difference in the solvation free energy between complex RT and corresponding subsystems calculated with the Generalized Born implicit solvation model.

12

and Npurine Npyrimidine

(5)

where Egas represents the electronic energy, ZPVE is the zeropoint vibrational energy, H is the enthalpy when passing from 0 to 300 K, and the last term in eq 5 stands for the entropy at room temperature. The solvation free energy is simply added to ΔGgas to calculate the binding free energy in implicit solvent (ΔGsol), which is defined as

1 − (rij/r0)6 1 − (rij/r0)

(4)

(3)

Sstack was calculated from the intermolecular distances rij between the nine ring atoms of purine (i) and the six ring atoms of pyrimidine (j). SWC acts on the donors/acceptors of WC H-bonding (two pairs of mA−mT and three pairs of mG− mC). The reference distance r0 were equal to 0.4 nm for Sstack and 0.3 nm for SWC. The MTD biased simulation along with the classical MD simulation was performed in Gromacs38 4.5.4 with a Plumed 1.1.0 plug-in.39 To verify the performance of the empirical potential used, all the minima obtained in the MTD simulations were reoptimized at the TPSS−D/TZVP level (only base pairs were considered). The DFT−D method used has been shown to provide reliable stabilization energies for H-bonded and stacked DNA base pairs40 and also for the peptides where intramolecular dispersion correction is necessary for the investigation of

Table 2. Interaction Energies of the Silica-Single-Base Complexes (in kJ/mol) Systems

TPSS−D/TZVP

B3LYP−D/TZVP

MM

S-Aa S-Gb S-Cc S-Td S−Ue

−59.83 −64.77 −46.10 −45.48 −43.40

−60.63 −65.81 −47.45 −46.74 −44.69

−64.10 −70.37 −53.77 −60.54 −53.47

a

silica-adenine. bsilica-guanine. csilica-cytosine. dsilica-thymine. esilicauracil.

C



Article

Figure 2. Free energy surface of the association of the mA−mT and mG−mC base pairs on the silica surface calculated with the negative of the MTD biased potential. The axis depicts the collective variables describing base stacking and WC H-bonding. All minima from A−K are highlighted.

Figure 3. All the snapshots (A−K) obtained from metadynamics. The minima A−E correspond to the mA−mT base pairs and the minima F−K correspond to the mG−mC base pairs. The minima D and H represent WC; B is the reverse WC for the mA−mT pair; A and F are the dissociated pairs; E and K depict the stack pairs for the mA−mT and mG−mC pairs, respectively.

obtained after the simulation are shown in Figure 2, where the different colors at each minimum correspond to the depths of free energy accumulated by the populations of base pair association in different configurations. The free energy calculated at each minimum represents the free energy of a set of configurations having similar values of CVs rather than the free energy of a single unique configuration. The selection of the snapshots (the most populated configuration) from each minimum (Figure 3) obtained by MTD with similar values of CVs close to the predicted minimum is therefore the most straightforward way of interpreting the results of MTD. To better illustrate the mutual orientation of the two bases on the surface, the purine bases were superimposed. The assembly of the base pairs, the free energy convergence of the stacked pair with respect to the H-bonded pair over time, and their flipping are shown in Figures 4, 5, and 6. A total number of 45 and 40 transitions between the H-bonded pair and the three-layer surface−π−π stacked assembly were observed during the 130 ns MTD simulation for both mA− mT and mG−mC base pairs, respectively. A number of transitions was lower in the case of the base pairs adsorbed at the graphene surface.10 This provides evidence of higher flexibility of base pairs at the silica surface. Multiple base flipping has been observed during the simulation. The orientation of the bases changes from the cis-orientation, suitable for WC pairing, to the trans-orientation, suitable for

stabilization energies are clearly more reliable than the MM ones. 3.2. MTD Free Energy. Base pair association and dissociation have been investigated using a biased MTD simulation. The biased MTD simulations predict that in the fully solvated system, the base pairs prefer H-bonded form rather than π−π stacked form (Table 4). The MTD biased simulation efficiently describes the association states using the collective variables (CV), where CV corresponds to the coordination number defined in eqs 2 and 3. The addition of the biased potential during the MTD simulation promoted transitions between the states where both nucleobases were stuck on the silica surface as H-bonded to the silica−purine−pyrimidine surface−π−π stacked assembly. While the silica−purine−pyrimidine stacked arrangement requires the simulations in 3D, that for silica H-bonded pair is reduced to 2D. Sampling in 2D is evidently much simpler than in 3D. The coordination numbers are high in the associated pair and close to zero when the pair is dissociated, no matter what their distance and mutual orientation are. A visual inspection of the trajectories revealed that all the association processes took place in 2D space on the silica surface. After the 130 ns of MTD simulation, different WC and non-WC base pairs were observed along with a three-layer surface−π−π stacked configuration. Non-WC base pairings often contain an H-bond of the WC pairing or at least a close proximity of the WC H-bond donor/acceptors. The FESs D



Article

bonded to the π−π stacked form and back to the reversed WC configuration. Several kinds of base flipping have been observed, the most common among which is the WC to the non-WC configuration, and the silica−purine−pyrimidine to the silica−pyrimidine−purine π−π stacked complex or vice versa. The pyrimidine base is highly movable and jumps very often (sometimes purine base also) from ∼0° to ∼180° to form a stack assembly of one upon another or from a WC to a nonWC configuration. The silica−pyrimidine−purine assembly has been observed during the simulation, but this flipping occurs much less than the previous base flipping. There were 51 and 58 flipping events observed for the mA−mT and mG−mC pairs, respectively. There was much more base flipping on the silica surface in comparison with the previously investigated base flipping on nanographene.10 This might arise from the more hydrophilic nature of the silica surface, containing several oxygen atoms, which create H-bonds with the bases during the flipping. Two kinds of the most probable base flipping for mA−mT pair in the aqueous media are shown in Figures 7 and 8 with the corresponding energy profiles from the WC to non-WC, and pyrimidine−purine π−π stack to purine−pyrimidine Tshaped (mostly the transition-state) transformation. Another frequent base flipping occurs between the WC H-bonded and π−π stacked complex, where the pyrimidine base directly flips by ∼180° over the purine base. During the flipping, the purine base forms a strong H-bond with the surface (with an average distance of ∼0.21 nm between the surface oxygen atom and N6−H). Table 3 shows the interaction energies of WC and non-WC structures. The non-WC structure is slightly less stable than WC one (Figure 7) due to the loss of the dispersion energy (3.21 kJ/mol) which might be explained by the different orientation of methyl groups in both structures. The T-shaped configuration is mainly formed when the purine base flips from ∼0° to ∼90° (Figure 7) or both base pairs flip simultaneously from ∼0° to ∼90° on the surface (Figure 8). The flipping between H-bonded and π−π stacked mA−mT base pair in solution was calculated41 at the CCSD(T)/CBS level in the absence of silica surface and barrier of ∼50 kJ/mol was found. In the present paper (TPSS−D/TZVP level) the same kind of flipping barrier (however in the presence of silica surface) was found to be higher, about 84 kJ/mol. The second barrier (from π−π stacked to the T-shaped) amounts to about 130 kJ/mol. The barrier height in the former case is almost exclusively determined by the dispersion energy difference (72 kJ/mol). Contribution of the dispersion energy is much lower in the second case (5 kJ/mol). Reliable description of base flipping can thus only be obtained if the computational method properly covers the dispersion energy. Another aim of this study is to calculate the free energy difference between the H-bonded and stacked structures of base pairs on the solvated surface. To obtain an accurate prediction of association free energy, it is necessary to take into account the effective concentration of the bases (the number of the bases and the size of the surface). After a long 130 ns MTD simulation for the mA−mT pair, five minima were observed on the FES (see Figure 3). Besides the WC base pairing, all the other H-bonded forms were found to be more stable than surface−π−π stacking configurations. The dissociated pair (A) turned out to be one of the most populated minimum for the mA−mT base pair. The WC base pairing (D) is less populated than the reverse WC (B) pair. Another minimum, called C, characterized by an interaction via N7 of the purine base as an

Figure 4. Development of Sstack (A) and the free energy difference between the stacked and WC structures (WC to stack free energy (B)) during the 130 ns MTD simulation of the association of the mA−mT base pair on the silica surface.

Figure 5. Development of Sstack (A) and the free energy difference between the stacked and WC structures (WC to stack free energy (B)) during the 130 ns MTD simulation of the association of the mG−mC base pair on the silica surface.

Figure 6. Orientation of the bases of the mA−mT and mG−mC pair during 130 ns MTD simulation on the silica surface.

reverse WC H-bonding, or vice versa. The angles specifying the plane orientation are 0° and 180° for the cis- and transorientation, respectively. Due to the flexibility of the silica surface and base pairs, there is a spontaneous flip from HE



Article

Figure 7. Flipping between the WC base and the non-WC base of the mA−mT pair on the surface (the silica surface is not shown here).

Figure 8. Flipping between the pyrimidine−purine surface−π−π stacked and surface-T configuration (the silica surface is not shown here).

acceptor (i.e., Hoogsteen and Hoogsteen-like base pairing) also turned out to be a more effective base pair association than WC pairing. The minima A, B, and C are all found to be in the cisorientation, whereas several stacked pairs (minimum E) were found in both cis- and trans-orientation. In the case of the mG−mC pair (see Figure 3), we found six free energy minima during the 130 ns MTD simulation. The dissociated pair (F) also turned out to be one of the most populated minima along with other minima toward the surface. The other minima similar to these in the mA−mT pair are found to be the WC H-bonded configuration (H) and the surface−π−π stacked configuration (K). The FES of the mG− mC pair is more complicated than that of the mA−mT pair which is caused by the presence of various H-bonded motifs. For example, the minimum (G) is associated with two H-bonds

between N1−H of mG and O2 of mC and between N2−H of mG and N3 of mC. The minimum (I) is characterized by the presence of two H-bonds between N1−H and, N2−H of mG, and O2 of mC. The minimum (J) is associated with only one H-bond between N2−H of mG and O2 of mC, which means that the pyrimidine base may flip and adopt different positions. The number of minima found in the present case is smaller than that in the case of base pairs adsorbed at the graphene surface.10 The minima G and J are quite dense as seen in Figure 2, and both of them possess higher binding free energy (more negative) than the WC configuration. Therefore, these minima contribute significantly to the binding between the base pairs and the surface and cannot be neglected. We have shown previously10 that in the case of the mG−mC pair, structures similar to G and J were less important than the WC H-bonded F



Article

Table 3. Interaction Energies Calculated with QM and MM for All Snapshots (in kJ/mol) TPSS−D/TZVPc(graphene)

TPSS−D/TZVP(silica) Eint

Systems

Edis

Aa B C D(WC) E

−101.67 −103.34 −97.86 −105.73 −61.17

(−1.67)b (3.81) (−4.06) (40.50)

−125.89 −127.58 −120.10 −130.79 −71.18

F G H(WC) I J K

−99.50 −99.03 −105.01 −100.96 −99.37 −71.59

(0.43) (−5.51) (−1.46) (0.13) (27.91)

−123.42 −122.38 −128.77 −121.82 −127.12 −76.60

a

MM (silica) mA−mT −137.23 −139.95 −137.82 −140.37 −81.88 mG−mC −130.37 −134.60 −135.64 −132.30 −134.30 −84.10

b

(−2.72) (−0.59) (−3.14) (55.35)

(−4.23) (−5.27) (−1.93) (−3.93) (46.27)

Eint

Edis

−112.48 −107.39 −107.15 −114.64 −59.87

−192.96 −198.23 −202.11 −200.92 −108.82

−114.92 −117.88 −119.88

−197.67 −200.73 −201.69

−75.80

−123.25 c

cf. Figure 3. Values in parentheses represent the interaction energy differences relative to the dissociated states A and F. Interaction energies were calculated at TPSS−D/TZVP level on the SCC−DFTB−D optimized geometry.

4. DISCUSSION The present calculations demonstrate that methylated bases under fully solvated conditions prefer to bind to the silica surface as H-bonded pairs rather than as stacked pairs. The empirical potential produces reliable results in comparison with the much more expensive QM results, which means that MM can be successfully used for the description of the base−base interactions on the silica surface. We have shown previously that the standard empirical potentials provide a relatively accurate modeling of the H-bonding and stacking interactions between the bases48,49 in the gas phase as well as on the standard graphene surface.10 An interesting result concerns the coexistence of multiple configurations at each minimum. The collective variables used in the MTD simulation can distinguish not only between separated, H-bonded, and stacked bases but also between different WC and non-WC pairings. The existence of multiple minima on the FES for the mA−mT and mG−mC pairs showed the necessity to introduce collective variables. Such minima as B, C, G, and J were found to be denser than WC pairing. Therefore, the binding free energy of different base pairs on the surface can be influenced by existence of several non-WC base pair associations. Experimentally, it has already been proven that the stacked mG−mC pair does not exist in the gas phase, whereas an addition of water molecules changes the configuration of their association as described in the Introduction.6,7 Accurate CCSD(T)/CBS2 calculations show that the WC pairing of the mA−mT and mG−mC complex is stronger than the respective stacked pairing by 13.4 and 44.0 kJ/mol, respectively. This finding suggests that the gas phase binding free energies of H-bonded complexes could be higher than these stacked pairs. Table 3 shows the TPSS−D/TZVP and MM interaction energies calculated for all the snapshots for each minimum found during the MTD simulation. The numbers in parentheses in columns 2 and 4 depict the difference in interaction energy relative to the dissociated state in both the mA−mT and mG−mC base pairs. The average difference in interaction energies (shown in parentheses in Table 3) calculated for the mA−mT and mG−mC pairs with the TPSS−D and MM methods amount to 9.65 and 12.23, and 4.30 and 7.73 kJ/mol, respectively. The average differences of

structure, which clearly demonstrates that the association of the bases strongly depends on the type of the surface. The QM TPSS−D/TZVP optimization was applied to all the minima obtained from the MTD simulation (see Figure 2). The main geometrical motifs determined at the MTD level remain unchanged upon passing to the QM level. The QM optimized geometries were further reoptimized at the MM level and interaction energies obtained were again compared with the QM one. Table 3 shows the interaction energies of all clusters along with the dispersion contribution calculated at the TPSS− D/TZVP and MM levels. The dissociated modes for both base pairs are supported by the water bridging if the explicit water model in the MTD simulation is considered, but they are missing if the gas phase QM optimization is performed. Table 4 Table 4. Binding Free Energies Calculated with MM(RR− HO−IG)sol in Solvent and the MTD Simulation for All Snapshots (in kJ/mol) Systems

MM(RR−HO−IG)sol

MTD

mA−mT A B C D(WC) E

−53.89 −61.50 −56.02 −57.49 −17.28

F G H(WC) I J K

−44.06 −50.29 −46.90 −46.86 −47.07 −12.76

(−4.01)a [−7.61]b [−2.13] [−3.60] (40.21) [36.61] mG−mC (−3.39) [−6.23] [−2.84] [−2.80] (−0.17) [−3.01] (34.14) [31.30]

−53.97 −46.32 −47.53 −44.47 −32.76

(−1.85) [7.65] [6.44] [9.50] (11.71) [21.21]

−42.76 −36.82 −35.10 −36.53 −37.66 −24.23

(−1.72) [5.94] [7.66] [6.23] (−0.84) [5.10] (10.87) [18.53]

a Values in parentheses depict the binding free energy differences relative to the WC pairs D and H. bValues in square bracket depict the binding free energy differences relative to the dissociated pair.

shows the binding free energy calculated with the MTD in explicit water and MM force field with implicit solvent. Harmonic frequencies were evaluated at the MM level (MM optimized geometries). Table 4 shows that the adsorption of the H-bonded structures of both base pairs is strongly preferred over the adsorption of π−π stacked structures. G



Article

interaction energies between these two methods for the mA− mT and mG−mC pairs are 2.58 and 3.43 kJ/mol, respectively, which suggests good agreement between the MM and QM methods. For the TPSS−D/TZVP method, the difference in the interaction energy of H-bonded pairs with respect to the dissociated pair varies from −4.06 to 3.81 kJ/mol and from −5.51 to 0.13 kJ/mol for the mA−mT and mG−mC pairs, respectively. For the MM method, it varies from −3.14 to −0.59 kJ/mol and from −5.27 to −1.93 kJ/mol for the mA− mT and mG−mC pairs, respectively. DFT−D as well as MM data show that in the case of mA−mT and mG−mC pairs, the WC and non-WC structures, and WC and other H-bonded structures, respectively, are more stable than the dissociated state. This means that both base pairs prefer to stay in the different H-bonded form rather than in the dissociated (stabilized by water-bridging) or stacked forms. Table 3 also shows interaction energies and their dispersion parts for the analogues structures of mA−mT and mG−mC pairs adsorbed at the graphene surface.10 It is evident that in most cases the stabilization energies at the graphene surface are larger and average difference between silica and graphene surfaces is 6.4 and 13.4 kJ/mol for mA−mT and mG−mC pairs, respectively. It is further evident that at both surfaces, the dispersion energy forms a dominant part of the total stabilization energy; the role of dispersion energy is even more pronounced at the graphene surface. Our calculation clearly predicts that the adsorption of individual bases on the surface is penalized on average by an energy of ∼40 kJ/mol with respect to the H-bonded base pair adsorption. The single base optimization with TPSS−D/TZVP and B3LYP−D/TZVP methods pointed out to the stack nature of binding. In this case, the bases are not fully parallel to the surface and contain one H-bond to the surface. The sum of the individual contributions of adenine and thymine, and guanine and cytosine, gives interaction energy which is almost equal to the interaction energy of the H-bonded pair. Table 4 depicts the free energy of the adsorption/binding of base pairs in different configurations calculated with MM using RR−HO−IG approximation. The results were thoroughly compared with the binding free energies calculated with the MTD simulation. Both RR−HO−IG calculations and MTD simulations show that, for the mA−mT pair, the minimum B (non-WC) is more attractive than the WC (D) pair by 4.01 and 1.85 kJ/mol (shown in parentheses in Table 4), whereas surface−π−π stacking (E) is less attractive by 40.21 and 11.71 kJ/mol, respectively than the WC pair. In the case of the mG− mC pair, the minima G and J are found to be more attractive by 3.39 and 1.72, and 0.17 and 0.84 kJ/mol, respectively, whereas surface−π−π stacking (K) was less attractive by 34.14 and 10.87 kJ/mol, respectively, than the WC pair (H) for both methods. Table 4 further shows that, in the case of mA−mT pair, both calculations prefer non-WC structure over the WC one and in the case of mG−mC pair, minima G and J are preferred over the WC one. The numbers in the square bracket in Table 4 refer to the binding free energy difference relative to the dissociated pair. Analyzing the MM values, we found that H-bonded structures of both base pairs are more stable than the dissociated pair, while the opposite is true for stacked structures. Let us add here that the same was found previously at the graphene paper.10 Table 4 (MM as well as MTD) tells us that non-WC structures of both base pairs are more stable than the respective WC structures. The opposite was found in our previous paper10

where instead of silica surface the graphene surface was considered. On the basis of SCC−DFTB−D50,51 (RR−HO− IG) calculations, it was shown there that the WC structures of both base pairs were more stable than the other planar Hbonded structures. The respective differences (between WC and other H-bonded structures) are, however, rather small, and they are similar (4 and 7 kJ/mol for silica and graphene base pairs, respectively). The binding free energies calculated in MM show that all WC and non-WC (Hoogsteen) base pairs have higher binding affinity toward the surface than the stacked pairs. The dissociated pairs A and F show lower binding affinity in comparison with other H-bonded configurations. MM calculations also show that the major contribution to the total binding free energy comes from the enthalpy part whereas the entropic contribution is almost negligible. It was difficult to consider the entropic contribution of the whole system (surface, base pairs, and also the solvent molecules) together in both QM and MM calculations in the aqueous medium. For example, in explicit water solvation, the dissociated pair is restrained by water molecules (water-bridging). The implicit solvent model does not include this effect of water bridging between base pairs and therefore it is not easy to compute the entropic contribution. Because the total entropy term is rather small, it will not make an important contribution to the total binding free energy.

5. CONCLUSIONS QM and MM calculations as well as biased metadynamics simulations demonstrate that the mG−mC and mA−mT base pairs are adsorbed on the fully solvated silica surface in an Hbonded form. Biased MTD simulations as well as RR−HO−IG calculations convincingly demonstrate that the adsorption of all the H-bonded structures of both base pairs is stronger than that of any π−π stacked structures. It was also shown that Hbonded structures of both base pairs are in the presence of fully solvated silica and graphene surfaces more stable than the stacked structures. Further, the WC, non-WC, and Hoogsteen H-bonded structures are adsorbed at solvated graphene surface more strongly than at solvated silica surface.

■

AUTHOR INFORMATION

Corresponding Author

*Tel.: +420 220 410311. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This work was supported by the Institute of Organic Chemistry and Biochemistry, Academy of Sciences of the Czech Republic [RVO: 61388963], and the Czech Science Foundation [P208/ 12/G016]. This work was also supported by the Operational Program Research and Development for Innovations−European Science Fund (CZ.1.05/2.1.00/03.0058). The support of Premium Academia of the Academy of Sciences of the Czech Republic awarded to P.H. in 2007 is also acknowledged. We further thank Prof. Piero Ugliengo for providing the structures of the silica surface.

■

REFERENCES

(1) Watson, J. D.; Crick, F. H. C. A Structure for Deoxyribose Nucleic Acid. Nature 1953, 171, 737−738.

H



Article

(2) Jurecka, P.; Hobza, P. True Stabilization Energies for the Optimal Planar Hydrogen-Bonded and Stacked Structures of Guanine···Cytosine, Adenine···Thymine, and Their 9- and 1-Methyl Derivatives: Complete Basis Set Calculations at the MP2 and CCSD(T) Levels and Comparison with Experiment. J. Am. Chem. Soc. 2003, 125, 15608− 15613. (3) Cerny, J.; Kabelac, M.; Hobza, P. Double-Helical→Ladder Structural Transition in the B-DNA is Induced by a Loss of Dispersion Energy. J. Am. Chem. Soc. 2008, 130, 16055−16059. (4) Kolar, M.; Kubar, T.; Hobza, P. On the Role of London Dispersion Forces in Biomolecular Structure Determination. J. Phys. Chem. B 2011, 115, 8038−8046. (5) Kratochvil, M.; Sponer, J.; Hobza, P. Global Minimum of the Adenine···Thymine Base Pair Corresponds Neither to Watson−Crick Nor to Hoogsteen Structures. Molecular Dynamic/Quenching/ AMBER and ab Initio beyond Hartree−Fock Studies. J. Am. Chem. Soc. 2000, 122, 3495−3499. (6) Kabelac, M.; Hobza, P. At Nonzero Temperatures, Stacked Structures of Methylated Nucleic Acid Base Pairs and Microhydrated Nonmethylated Nucleic Acid Base Pairs are Favored over Planar Hydrogen-Bonded Structures: A Molecular Dynamics Simulations Study. Chem.Eur. J. 2001, 7, 2067−2074. (7) Kabelac, M.; Zendlova, L.; Reha, D.; Hobza, P. Potential Energy Surfaces of an Adenine−Thymine Base Pair and Its Methylated Analogue in the Presence of One and Two Water Molecules: Molecular Mechanics and Correlated ab Initio Study. J. Phys. Chem. B 2005, 109, 12206−12213. (8) Ferris, J. P. Montmorillonite Catalysis of RNA Oligomer Formation in Aqueous Solution. A Model for the Prebiotic Formation of RNA. J. Am. Chem. Soc. 1993, 115, 12270−12275. (9) Ferris, J. P. Mineral Catalysis and Prebiotic Synthesis: Montmorillonite-Catalyzed Formation of RNA. Elements 2005, 1, 145−149. (10) Spiwok, V.; Hobza, P.; Rezac, J. Free-Energy Simulations of Hydrogen Bonding versus Stacking of Nucleobases on a Graphene Surface. J. Phys. Chem. C 2011, 115, 19455−19462. (11) Pelmenschikov, A.; Leszczynski, J. Adsorption of 1,3,5Trinitrobenzene on the Siloxane Sites of Clay Minerals: Ab Initio Calculations of Molecular Models. J. Phys. Chem. B 1999, 103, 6886− 6890. (12) Rimola, A.; Civalleri, B.; Ugliengo, P. Physisorption of Aromatic Organic Contaminants at the Surface of Hydrophobic/Hydrophilic Silica Geosorbents: a B3LYP−D Modeling Study. Phys. Chem. Chem. Phys. 2010, 12, 6357−6366. (13) Mignon, P.; Ugliengo, P.; Sodupe, P. Theoretical Study of the Adsorption of RNA/DNA Bases on the External Surfaces of Na+Montmorillonite. J. Phys. Chem. C 2009, 113, 13741−13749. (14) Aquino, A. J. A.; Tunega, D.; Gerzabek, M. H.; Lischka, H. Modeling Catalytic Effects of Clay Mineral Surfaces on Peptide Bond Formation. J. Phys. Chem. B 2004, 108, 10120−10130. (15) Laio, A.; Parrinello, M. Escaping Free-Energy Minima. Proc. Natl. Acad. Sci. U.S.A. 2002, 99, 12562−15566. (16) Bussi, G.; Gervasio, F. L.; Laio, A.; Parrinello, M. Free Energy Landscape for ß Hairpin Folding from Combined Parallel Tempering and Metadynamics. J. Am. Chem. Soc. 2006, 128, 13435−13441. (17) Piana, S.; Laio, A. A Bias-Exchange Approach to Protein Folding. J. Phys. Chem. B 2007, 111, 4553−4559. (18) Ensing, B.; Vivo, M, D.; Liu, Z.; Moore, P.; Klein, M. L. Metadynamics as a Tool for Exploring Free Energy Landscape of Chemical Reactions. Acc. Chem. Res. 2006, 39, 73−81. (19) Barducci, A.; Bonomi, M.; Parrinello, M. Metadynamics. WIREs. Comput. Mol. Sci 2011, 1, 826−843. (20) Laio, A.; Gervasio, F. L. Metadynamics: A Method to Simulate Rare Events and Reconstruct the Free Energy in Biophysics, Chemistry and Material Science. Rep. Prog. Phys. 2008, 71, 126601. (21) Sutto, L.; Marsili, S.; Gervasio, F. L. New Advances in Metadynamics. WIREs. Comput. Mol. Sci 2012, 1, 771−779.

(22) Tosoni, S.; Civalleri, B.; Ugliengo, P. Hydrophobic Behavior of Dehydroxylated Silica Surfaces: A B3LYP Periodic Study. J. Phys. Chem. C 2010, 114, 19984−19992. (23) Cruz−Chu, E. R.; Aksimentiev, A.; Schulten, K. Water−Silica Force Field for Simulating Nanodevices. J. Phys. Chem. B 2006, 110, 21497−21508. (24) Hohenberg, P.; Kohn, W. Inhomogeneous Electron Gas. Phys. Rev. 1964, 136, 864−871. (25) Kohn, W.; Sham, L. J. Self-Consistent Equations Including Exchange and Correlation Effects. Phys. Rev. 1965, 140, 1133−1138. (26) Cerny, J.; Jurecka, P.; Hobza, P.; Valdes, H. Resolution of Identity Density Functional Theory Augmented with an Empirical Dispersion Term (RI−DFT−D): A Promising Tool for Studying Isolated Small Peptides. J. Phys. Chem. A 2007, 111, 1146−1154. (27) Wang, J. M.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (28) Wang, J.; Wang, W.; Kolmann, P. A.; Case, D. A. Automatic Atom Type and Bond Type Perception in Molecular Mechanical Calculations. J. Mol. Graphics Modell. 2006, 25, 247−260. (29) Bayly, Ch. I.; Cieplak, P.; Cornell, W.; Kollman, P. A. A WellBehaved Electrostatic Potential Based Method Using Charge Restraints for Deriving Atomic Charges: the RESP Model. J. Phys. Chem. 1993, 97, 10269−10280. (30) Valdes, H.; Pluhackova, K.; Pitonak, M.; Rezac, J.; Hobza, P. Benchmark Database on Isolated Small Peptides Containing an Aromatic Side Chain: Comparison between Wave Function and Density Functional Theory Methods and Empirical Force Field. Phys. Chem. Chem. Phys. 2008, 10, 2747−2757. (31) Duan, Y.; Wu, C.; Chowdhury, S.; Lee, M. C.; Xiong, G.; Zhang, W.; Yang, R.; Cieplak, P.; Luo, R.; Lee, T.; et al. A Point-Charge Force Field for Molecular Mechanics Simulations of Proteins Based on Condensed-Phase Quantum Mechanical Calculations. J. Comput. Chem. 2003, 24, 1999−2012. (32) Ahlrichs, A.; Bar, M.; Haser, M.; Horn, H.; Kolmel, C. Electronic Structure Calculations on Workstation Computers: The Program System Turbomole. Chem. Phys. Lett. 1989, 162, 165−189. (33) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Mennucci, B.; Petersson, G. A. Gaussian 09, revision A.1; Gaussian, Inc.; Wallingford, CT, 2009. (34) Case, D. A.; Darden, T. A.; Cheatham, III, T. E.; Simmerling, C. L.; Wang, J.; Duke, R. E.; Luo, R.; Crowley, M.; Walker, R. C.; Zhang, W. et al. AMBER 10; University of California:: San Francisco, 2008. (35) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (36) Nosé, S. A Molecular Dynamics Method for Simulations in the Canonical Ensemble. Mol. Phys. 1984, 52, 255−268. (37) Hoover, W. G. Canonical Dynamics: Equilibrium Phase-space Distributions. Phys. Rev. A 1985, 31, 1695−1697. (38) Hess, B.; Kutzner, C.; Van der Spoel, D.; Lindahl, E. GROMACS 4: Algorithms for Highly Efficient, Load-Balanced, and Scalable Molecular Simulation. J. Chem. Theory Comput. 2008, 4, 435− 447. (39) Bonomi, M.; Branduardi, D.; Bussi, G.; Camilloni, C.; Provasi, D.; Raiteri, P.; Donadio, D.; Marinelli, F.; Pietrucci, F.; Broglia, R. A.; Parrinello, M. PLUMED: A Portable Plugin for Free-Energy Calculations with Molecular Dynamics. Comput. Phys. Commun. 2009, 180, 1961−1972. (40) Cerny, J.; Hobza, P. Energy Barriers between H-bonded and Stacked Structures of 9-Methyladenine...1-Methylthymine and 9Methylguanine...1-Methylcytosine Complexes. Chem. Commun 2010, 46, 383−385. (41) Valdes, H.; Pluhackova, K.; Hobza, P. Phenylalanyl-GlycylPhenylalanine Tripeptide: A Model System for Aromatic-Aromatic Side Chain Interactions in Proteins. J. Chem. Theory Comput. 2009, 9, 2248−2256. I



Article

(42) Gloaguen, E.; Valdes, H.; Pagliarulo, F.; Pollet, R.; Tardivel, B.; Hobza, P.; Piuzzi, F.; Mons, M. Experimental and Theoretical Investigation of the Aromatic−Aromatic Interaction in Isolated Capped Dipeptides. J. Phys. Chem. A 2010, 114, 2973−2982. (43) Boys, S. F.; Bernardi, F. The calculation of Small Molecular Interactions by the Differences of Separate Total Energies. Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−566. (44) Klamt, A.; Schuurmaan, G. COSMO: A New Approach to Dielectric Screening in Solvents with Explicit Expressions for the Screening Energy and its Gradient. J. Chem. Soc., Perkin Trans. 2 1993, 5, 799−805. (45) Podolsky, B. Quantum-Mechanically Correct Form of Hamiltonian Function for Conservative Systems. Phys. Rev. 1928, 32, 812−816. (46) Hawkins, G. D.; Cramer, C. J.; Truhlar, D. G. Parametrized Models of Aqueous Free Energies of Solvation Based on Pairwise Descreening of Solute Atomic Charges from a Dielectric Medium. J. Phys. Chem. 1996, 100, 19824−19839. (47) Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on the Generalized Born Approximation with Asymmetric Describing. J. Chem. Theory Comput. 2009, 5, 2447−2464. (48) Riley, K. E.; Pitonak, M.; Jurecka, P.; Hobza, P. Stabilization and Structure Calculations for Noncovalent Interactions in Extended Molecular Systems Based on Wave Function and Density Functional Theories. Chem. Rev. 2010, 110, 5023−5063. (49) Hobza, P.; Kabelac, M.; Sponer, J.; Mejzlik, P.; Vondrasek, J. Performance of Empirical Potentials (AMBER, CFF95, CVFF, CHARMM, OPLS, POLTEV), Semiempirical Quantum Chemical Methods (AM1, MNDO/M, PM3), and ab Initio Hartree−Fock Method for Interaction of DNA Bases: Comparison with Nonempirical beyond Hartree−Fock Results. J. Comput. Chem. 1997, 18, 1136−1150. (50) Elstner, M.; Hobza, P.; Frauenheim, T.; Suhai, S.; Kaxiras, E. Hydrogen Bonding and Stacking Interactions of Nucleic Acid Base Pairs: A Density-Functional-Theory Based Treatment. J. Chem. Phys. 2001, 114, 5149−5155. (51) Elstner, M.; Poreza, G. D.; Jungnickel, G.; Elstner, J.; Haugk, M.; Frauenheim, T.; Suhai, S.; Seifert, G. Self-Consistent-Charge DensityFunctional Tight-Binding Method for Simulations of Complex Materials Properties. Phys. Rev. B 1998, 58, 7260−7268.

J


On the Association of the Base Pairs on the Silica Surface Based on

Recommend Documents