J. Phys. Chem. B 2007, 111, 11083-11089
11083
A Molecular Dynamics Calculations of Hole Transfer Rates in DNA Strands Maksym Volobuyev,*,† Humberto Saint-Martin,‡ and Ludwik Adamowicz† Nano-Biomolecular Engineering Science and Technology (n-BEST) Program and Department of Chemistry, UniVersity of Arizona, Tucson, Arizona 85721, and Centro de Ciencias Fsicas, UniVersidad Nacional Auto´ noma de Me´ xico, Apartado Postal 48-3, CuernaVaca, Morelos 62251, Me´ xico ReceiVed: NoVember 29, 2006; In Final Form: January 29, 2007
A computational model, which includes both tunneling and thermal hopping mechanisms, has been applied to study the charge transfer in DNA (GC)n and (AT)n strands. The calculations revealed the crucial role palyed by the A or G NH2-group vibrations in the hole transfer in both types of strands. Charge-transfer rates in the two strands have been determined based on the molecular dynamics calculations. They are in good agreement with the available experimental data. The modeling appoach used here may be employed in the theoretical study of the charge transfer in natural and artificial DNA strands containing AT and GC pairs.
1. Introduction The charge transfer through DNA strands has been intensively studied in the past decade. This is because some important applications of this process play a crucial role in biology and nanotechnology. The experimental studies of this phenomenon have provided a large body of data concerning various aspects of the transfer process. At this point a refined model of the DNA charge-transfer process is needed to describe the phenomenon and explain the accumulated experimental results. The work on such a model has been carried out in our research group and the results presented here describe the progress we have made in this development. Due to high interest in charge transfer in aperiodic organic polymers like DNA, a model needs to be developed that not only explains the mechanism for the chargetransfer process but can also be used as a predictive tool to design molecular systems based on the DNA template with a specific rate for the charge transfer. Modeling and subsequent synthesis of such systems may provide biomolecular materials for various applications, particularly in the area of nanoelectronics. The charge in a molecular system can be transferred through the movement of negative or positive charges or both. In the former case electrons are the charge carriers and in the latter the holes. If the molecular system contains atoms or groups of atoms that strongly attract electrons, they can affect or even stop the electron transfer by electron trapping. Similar trapping can also occur for holes. In DNA, the nucleic acid bases such as cytosine, thymine, and uracil have higher electron affinities than purine bases and the electron transfer occurs mainly through those pyrimidine bases. In contrast, the transfer of holes occurs primarily through guanine (or adenine) molecules because that molecule has the lowest ionization potential of all DNA nucleobases. One of the differences between the two processes of electron transfer and hole transfer is the channels the two charge carries propagate in a molecule. In the former case the channels consist of HOMOs and in the latter case of LUMOs. Due to this difference, the electron and hole transport problems have been treated separately in theoretical calculations,1,2 even † ‡
University of Arizona. Universidad Nacional Auto´noma de Me´xico.
though the two mechanisms in some situations may operate simultaneously. The charge transfer at the molecular level should be described as a quantum mechanical (QM) event since, at that level, the available energy states are usually discrete and this important feature has to be taken into consideration. There are different levels of theory that can be used in quantum mechanical calculations concerning the charge transfer. Obviously, since rather large molecular systems are being considered, that level cannot be very high. An important element of a quantum calculation of the charge-transfer rate in DNA is to determine the electronic coupling between adjacent nucleobases, since the rate constant for charge-transfer process depends on this coupling. Within the quantum mechanical approach there are different models that one can apply. One of the most commonly used is the model involving polarons.3 As mentioned, the most serious limitation of the quantum approach is the size of the molecular system that can be considered. Taking into account only one Watson-Crick nucleobase pair with a sugar part gives a system containing more than 60 atoms,1 which is rather large for a QM calculation. In the experiment, however, the charge transfer is usually investigated in DNA fragments with at least 10 Watson-Crick pairs. Thus, in order for a theoretical calculation to be relevant, the models have to describe systems of similar lengths. There is also an additional difficulty that appears in a QM modeling of the charge transfer and is related to the fact that in the QM approach the hole or electron movement between molecular fragments of the studied system occurs through tunneling. The charge tunneling operates over a relatively short distance since its probability decays exponentially. The charge transfer over a long distance can be described by taking into account the thermal hopping mechanism.4 This mechanism operates through structural fluctuations of the DNA system that intensify where the temperature is raised. These fluctuations affect the internal vibrational motion for each base and, to an even larger extent, the interbase motions. The interbase and intrabase fluctuations may create conditions when an electron or a hole can tunnel between two adjacent bases. To describe both the structural dynamics and the electron or hole transfer dynamics quantum mechanically in a single model is a complicated task that is difficult to solve even for very small systems.
10.1021/jp067908h CCC: $37.00 © 2007 American Chemical Society Published on Web 08/28/2007
11084 J. Phys. Chem. B, Vol. 111, No. 37, 2007 Due to the complexity of the charge-transfer problem in a structurally flexible system like DNA, we have focused on an approach that is a hybrid of molecular dynamics and quantum mechanical methods. The main advantage of such an approach is the inclusion of the time dependency in the description of the state of the system, yet treating the important aspects of the problem at the quantum mechanical level. This feature is important to describe the dynamics of the charge-transfer problem. Also due to its reduced complexity the molecular dynamics (MD)/QM approach can be applied to relatively large systems using modest computational resources. The low cost of the method allows us to describe the interaction of the charge transferring system with its molecular surroundings, which is an important effect that needs to be accounted for in a realistic molecular-conduction model. Thus, the MD/QM approach can be used to describe a long-distance charge-transfer that occurs via a thermal hopping mechanism. This type of mechanism dominates in structurally flexible, aperiodic organic polymers such as DNA. Our first work on the electron transfer in DNA5 concerned (AT)n + water systems. The transfer rate constants obtained in that work were in good agreement with data from other studies. The following work6 focused on the hole transport through (GC)n strands submerged in water. It was a more computationally demanding task because guanine as a purine nucleobase is larger than thymine, whose structure is based on the pyrimidine ring. In the present work we consider mixed strands that contain both adenine and guanine nucleobases. The presence of GG or GGG sites in the strand plays a crucial role in the hole transfer because such sites trap positive charges due to the relatively low ionization potential of guanine clusters.7,8 2. The Model Our charge-transfer calculations are based on MD simulations. We consider systems that are single DNA strands surrounded by water molecules. The presence of water stabilizes the strand and represents the natural DNA environment in biological systems and systems where the DNA charge transport is usually studied in experiments. The size of the MD box in the calculations has been selected using the following rules: the distance between any atom of the DNA molecule and the surface of box cannot be less than 8 Å, and the number of water molecules should be sufficient to fill up the whole space within the box not occupied by the DNA strand. The initial position of the DNA strand in the box was always near the center. The solvation effects have been explicitly included in the MD simulation. For generating the MD trajectory we have used the AMBER package.9 The Cornell et al.10 force field has been used for the DNA strand and the TIP3 force field for the water molecules. Before the simulation of the charge-transfer event was performed, the system was first fully equilibrated. The time step in the simulation was 1 fs, and the chosen thermal bath parameter was 2 ps. As shown before,2 a cation state (hole) in a DNA strand localizes on a single purine base. Since the interplay between cation and neutral states of the purine base (guanine and adenine) in the strand is crucial in the dynamic model of the hole transport in DNA, these states have to be more accurately represented in the simulation than in a standard MD implementation. Our approach to this was described in our previous papers5,6 where we developed an approximate method to deal with this problem. Below we briefly describe the approach in relation to the calculations we have performed in this work. Important properties that affect the process of the transport of a hole along the DNA strand are the electron affinity of the
Volobuyev et al.
Figure 1. N9-Methyladenine (left) and N9-methylguanine (right).
purine base, where the hole is located at the particular moment in time in the MD simulation, and the ionization potentials of the two neighboring bases. The relation between these three quantities determines whether the hole will or will not move. Both the electron affinity of the cationic base and the ionization potential of the neutral base depend on the geometry of the bases in the particular moment in time in the MD simulation (the geometry of the base undergoes thermal fluctuations that intensify at higher temperatures; those fluctuations occur on the time scale of the MD simulation). In order to determine the ionization potential of an adenine base or a guanine base embedded in the DNA strand, we extracted the base from the strand and we compared its structure with the equilibrium structure of an isolated base. In this comparison we used an isolated adenine or a guanine molecule methylated at the N9 position, because such a system more closely resembles the base in the DNA strand attached at the N9 position to the DNA backbone. The equilibrium structures of the N9-methylated adenine and guanine molecules are shown in Figure 1. Let us suppose that the geometry of the nucleobase molecule extracted from the DNA strand at particular moment of the MD simulation is different from the vacuo-equilibrium geometry of that base. We can say that the nucleobase was deformed by the thermal motion. Such a deformation can be always described as a result of a structural vibration. The vibration is either a distortion along a single normal mode or, as it is in most cases, a simultaneous distortion along several (or all) normal modes. Thus, the thermally distorted geometry of the base molecule, g(t), at a particular moment in time, t, can be represented as a sum of the equilibrium geometry of the N9-methylated base, geq, and contributions from distortions along all (M) normal modes of the system, gidisp, multipled by the distortion parameters, Ri(t) M
g(t) ) geq +
Ri(t)gidisp ∑ i)1
(1)
If such a decomposition is known, it can be used to approximate the total energy of the distorted system provided that the energy functions, Ei(Ri) ) f(Rigdisp), are also known for each normal mode within the necessary interval of the distortion parameter Ri. In this we assume that the contribution from terms that couple two or more normal modes can be neglected. The multipliers, Ri, can be determined by projecting the mass weighted vector describing the displacement of the geometry from the equilibrium structure
∆g ) xm(g(t) - geq)
(2)
onto the normal coordinate, Qi
Ri ) Qi∆g
(3)
The functions Ei(Ri) ) f(Rigdisp) were generated for each normal mode of N9-methyladenine. The following sets of
Hole Transfer Rates in DNA Strands
J. Phys. Chem. B, Vol. 111, No. 37, 2007 11085
TABLE 1: N9-Methylguanine (G) and N9-Methyladenine (A) Energies in Vacuo at Different Equilibrium Geometriesa
G0
G+ ∆E A0 A+ ∆E
at neutral equilibrium geometry
at cation equilibrium geometry
-581.879453 -581.602296 7.54 -506.646030 -506.357153 7.86
-581.868719 -581.615639 6.89 -506.637388 -581.366298 7.38
neutral molecule
a ∆E values are IPs for molecules at neutral geometry and EA for cation geometry. (Total energies in au; ∆E in eV.)
functions were first determined: (1) the energy of the neutral N9-methyladenine as a function of the displacement along each of its normal modes from the N9-methyladenine equilibrium geometry; (2) the energy of the neutral N9-methyladenine as a function of the displacement along each normal mode of the N9-methyladenine cation from the equilibrium geometry of the cation; (3) the energy of the N9-methyladenine cation as a function of the displacement along each normal mode of the neutral N9-methyladenine from the equilibrium geometry of the neutral; (4) the energy of the N9-methyladenine cation as a function of the displacement along each of its normal mode form the cation equilibrium geometry. To determine the above functions, the energy of the system was calculated at at least 31 single points for each normal mode in QM calculations. The range for Ri for each normal mode, i, was obtained from a preliminary MD run that lasted just over 300 ps. Then geometries for the single-point QM calculations were generated by adding the coordinates of the particular normal mode scaled by incrementally varied value of R to the coordinates of the respective equilibrium geometry in the “standard orientation” of the molecule
gji ) geq + RjQi,
Rj ∈ [Rmin ‚‚‚ Rmax]
TABLE 2: Atomic Charges (in au) of Adenine Used in the Standard AMBER Parametrization and the Charges Used in Our Model
(4)
The single-point calculations were done for each mode. Due to the deviations from the equilibrium geometry being generally small during thermal fluctuations at room temperature, a polynomial fitting for each of the E(R) functions was used. Since there is usually some degree of anharmonicity of the energy for each mode, a polynomial fit beyond a simple quadratic function had to be used. For including terms with powers higher than 2, we used the Chebyshev polynomial method. As is well-known, the Chebyshev orthogonal polynomials provide the most optimal polynomial representation of an analytical function.11 In the energy expansions in terms of R, we included powers determined based on the respective Chebyshev representations. The polynomial fits determined this way had the highest standard deviations for all modes around 10-6 au with the highest degree of the polynomial used being 10. When the power of the polynomial excited 6, the number of single point energy calculations for generating the fit was increased to 81. The quantum mechanical calculations for generating the energy points for both the neutral N9-methyladenine molecule and its cation were performed in the same way as in the calculations for N9-methylguanine in our previous work.6 The B3LYP/6-31G** level of theory implemented in the GAUSSIAN 98 program package12 was used here both for the energy single-point calculations and for the vibrational analysis. The original atomic charge distribution in the AMBER parametrization of the adenine molecule was corrected by the standard RESP procedure implemented in the R.E.D-II program13
N9 C8 H8 N7 C5 C6 N6 H61 H62 N1 C2 H2 N3 C4
AMBER
our model
cation our model
-0.026800 0.160700 0.187700 -0.617500 0.072500 0.689700 -0.912300 0.416700 0.416700 -0.762400 0.571600 0.059800 -0.741700 0.380000
-0.046650 0.051072 0.152157 -0.495222 0.122811 0.445556 -0.707049 0.349653 0.364776 -0.601295 0.422427 0.054035 -0.591329 0.373758
-0.065996 0.221505 0.194469 -0.434116 0.154101 0.471059 -0.481567 0.369508 0.402066 -0.569951 0.552541 0.102562 -0.534037 0.512556
(Table 2). All other force field parameters for guanine and adenine remained unchanged. The results of the QM calculations for N9-methylguanine were taken from our previous work. To describe the charge migration along the DNA chain we used a simple physical model based on two elemental chargetransfer steps. The first step involved the charge tunneling, which is the main mechanism for charge transfer over a short distance. The probability to cross the barrier V for a particle with energy E (see Figure 1 in Supporting Information) can be estimated by the factor14
χ)
16 exp(-ka) (1 + Z2)(1 + Z-2)
(5)
where k2 ) (2m(V - E))/p2, Z ) ((V - E)/V)1/2 and where a is the barrier width. As mentioned above, an excess unit charge (positive or negative) in a DNA strand is localized on a single nucleobase so in our model a hole will migrate to the next nucleobase if χ is greater than 0.5. If it is not greater than 0.5, no charge migration will occur. The distance a between the two adjacent nucleobases transferring a hole between each other in the relation (5) was chosen as the distance of two closest atoms of the two bases minus the atomic radii of the atoms. Our calculations showed that the tunneling permits only carrying a hole from guanine to adenine. Also, due to an exponential decay of the tunneling probability we allow the charge transfer to only occur between two adjacent bases. No tunneling was allowed between bases separated by one or more other bases. The second mechanism for the charge-transfer considered here was the thermal hopping. In this process a hole transfer occurs only if the following transition condition is satisfied in the simulation
|VEA(A)| - |VIP(B)| > td
(6)
where td is a threshold, VEA(A) is the vertical electron affinity of the base with the hole, and VIP(B) is the vertical ionization potential for either the neutral nucleobase placed to the left of the base with the hole or the one placed to the right. If condition 6 was satisfied for both neighboring nucleobases, the hole was moved to the one with the lower VIP. Using different thresholds for charge transfer in the forward and backward directions, we can reproduce the experimental fact that, while the hole may move in both directions,15 it will move more likely toward the negative field direction than in the opposite direction, if an external electric field is present.
11086 J. Phys. Chem. B, Vol. 111, No. 37, 2007
Volobuyev et al. TABLE 3: Sequences of DNA Used in the Simulations
Figure 2. Energy difference (au) as a function of the nucleobase position in strand AAAGGGAAAA.
3. Results and Discussion 3.1. Static Properties. Two important characteristics of a nucleobase in a DNA strand are its ionization potential (the energy associated with removing an electron from the base) and its electron affinity (the energy associated with putting an excess electron on the base). Thymine has the highest electron affinity of all four nucleobases, A, C, G, and T, forming DNA. Thus, if one considers transfer of an excess electron through a DNA strand, the interaction of the electron with thymines in the strand is the key property. Inversely, since the positive charge (hole) has preference to migrate over nucleobases with low ionization potential like guanine and adenine, the interaction of the hole with those two bases has to be considered. The preference of excess electrons to move through the pyrimidine bases and of excess holes to move through the purine bases in DNA has been experimentally studied by the Giese group7,16 among others. Our calculated values for ionization potentials (IPs) for adenine and guanine (Table 1) are in good agreement with both the experimental and the previously calculated data.17-19 It should be noted that guanine has the lowest IP of all the DNA nucleobases. Also one can notice that adenine has a lower difference between VIP and VEA than guanine. This makes it easier for a hole to transfer between two adenine molecules by the thermal hopping mechanism. The difference between VIP and VEA for adenine, |VEA(A+)| - |VIP(A0)| ) 0.16 eV, is three times smaller than that for guanine, |VEA(G+)| - |VIP(G0)| ) 0.48 eV. If such a VEA/VIP relation between adenine and guanine occurs not only at the equilibrium geometry but also at geometries distorted from the equilibrium, it makes the hole adenine f guanine hopping more likely than the adenine f adenine hopping. For example, if a guanine molecule is located in the DNA strand before a positively charged adenine molecule and an neutral adenine molecule is located after the one with the positive charge (GA+A), it is highly probable that the positive charge will migrate in the backward direction to the guanine, even if this direction of migration is against the direction of the external electric field applied to the DNA molecule. The VEA/VIP differences for a hole transfer from a guanine cation to another guanine and to an adenine are |VEA(G+)| - |VIP(G0)| ) 0.65 eV and |VEA(G+)| - |VIP(A0)| ) 0.97 eV, respectively. At arbitrary geometries of the cation molecule (guanine) and the neighboring neutral adenine or guanine molecules, these differences may vary to some extent. But even with those variations the thermal hopping of a hole between G+ and A0 sites is almost impossible. Hence, only tunneling may move a hole from G+ to A0. Figure 2 shows energy differences, ∆E (VIP for the neutral molecule and VEA for the cation), that occur in the first
5′3′-
A T
A T
A T
A T
compound 1 A A A T T T
A T
A T
A T
-3′ -5′
5′3′-
G C
G C
G C
G C
compound 2 G G G C C C
G C
G C
G C
-3′ -5′
5′3′-
A T
A T
A T
G C
compound 3 G G A C C T
A T
A T
A T
-3′ -5′
nucleobase molecule in the stand (compound 3; see Table 3) averaged over 100 ps. In the calculation that produced those results the position of the hole was fixed so even when a favorable condition for tunneling (5) or for hopping (6) occurred the charge was not moved and remained at the same nucleobase. Since there are ten nucleobases in the strand we studied, halting the hole at each of those positions would produce ten different ∆E profiles. However, as we determined, within a statistical error of ≈0.001 au the profiles are the same. This means that hole presence or absence on a particular nucleobase does not change the average energy difference between the cation and the neutral state of the base. In our model the cation state of adenine differs from the neutral state only by the charge distribution (Table 2). This led to a difference only in the Coulombic part of the total energy while the other parameters involved in the potential such as the bonds, the bond angles, the dihedral angels, and the Lennard-Jones constances remained the same. One could consider varying also those fixed parameters to improve the model, but this would lead to much more costly calculations. The ∆E plots shown in Figure 2 clearly reveal the difference between the adenine and guanine molecules in the strand. The average ∆E for adenine is always higher than that for guanine, and a hole almost always prefers to be located on a guanine molecule in correspondence to the minimum-energy principle. Hole transfer from adenine to guanine is a likely event, but at the same time hole transfer in the opposite direction is impossible by thermal hopping. An averaged gap of 0.02 au is too big for such a mechanism to operate, and this makes guanine molecules likely traps for holes in the strand.8 The somewhat unexpected behavior of ∆E at the first adenine molecule in the strand can be explained by edge effects. Both 3′ and 5′ ends of the strand are more exposed to the influence of the surrounding medium (in present calculations we used water as the medium), so vibrations of atoms in nucleobases near the strand ends can be somewhat different from the vibrations of the nucleobases inside the strand. We found such differences only for the 5′ nucleobases located at the edge of the strand in all the systems we studied. The edge effect can be reduced or even eliminated by placing the hole not at the first nucleobase in the strand but at the second one in the first step of the charge-transfer modeling calculation. The ∆E function provides an approximate profile of the potential for a hole transport in the strand. Such information can be also provided, for example, by the minimal and maximum ∆E values. In Figure 3, the dash line represents the well-known shape structure obtained by shifting the maximal value of ∆E. We changed the positions of the maxima by shifting them by 0.5 to the left to obtain a more readable plot. Without shifting the values of the maxima and the minima, the curve could belong to the same strand position and, in such a case, the figure would contain only separated vertical lines instead of spline curves. At the same plot we also included the averaged ∆E from the previous figure. It should be noted that the precision of the
Hole Transfer Rates in DNA Strands
J. Phys. Chem. B, Vol. 111, No. 37, 2007 11087 time moments are not the same. Thus, the average ∆Ei for each normal mode was calculated as
∆Ei ) 1/t
Figure 3. Energy difference profile (dash line) based on ∆E limit values at 300 K averaged over 100 ps. Solid line is the same as shown in Figure 2.
values of the maximum and the minimum is lower than that for the averaged values. But even at that level of precision the values of the maximum and the minimum can be used to draw the following qualitative conclusions. 1. ∆E for both A and G nucleobases can vary in a relatively wide (≈0.1 au) range. This shows the importance of the molecular vibrations for the charge transfer. 2. In the previous work6 we showed that in some circumstances there is a possibility for a backward charge transfer against the direction of the external field. We described the conditions for such an event to occur. For some positions of the hole (for instance, 3, 5, and 8) the hole transfer to the previous position lowers the total energy due to an effect that we called “thermodynamic”. This effect was related to lowering the ∆E value for each nucleobase in the strand. At the same time charge transfer to the next position from, for example, the position 5 needs to overcome a much lower barrier than a charge transfer in the backward direction due to the factor which we called “kinetic”. The combination of the thermodynamic and kinetic factors that occur at a particular moment in time of the MD simulation may result in either forward or backward hole transfer that can occur due to either thermal hopping or tunneling. An increase of the temperature of up to 325 K or even 350 K does not change the average ∆E values and their upper and lower points for either nucleobase in the strand. We attribute this lack of the sensitivity in the temperature to the fact that ∆E is dependent on the strong electrostatic and covalent forces of the electron binding. This last conclusion is consistent with the previous observation that the hole presence or absence (i.e., which leads to different atomic charge distributions in the nucleobase; see Table 2) does not noticeably change ∆E. It should also be noted that the adenine and guanine molecules are aromatic systems that can strongly interact in the stacked orientation due to the overlapping of their π electrons. Also, the NH2 group which is not a part of the aromatic ring can create a strong hydrogen bond with an oxygen atom of the conjugated nucleobase (cytosine for guanine and thymine for adenine). Temperature change affects the vibration intensity for each normal mode. Some intensities increase and some decrease. Also ∆E can decrease or increase due to the different ways the temperature affects different normal modes behavior. Using the geometry decomposition (1), we can calculate ∆Ei for each normal mode. After time averaging of ∆Ei it is possible to estimate the contribution of each normal mode to the total ∆E value and, hence, to the charge-transfer process. However a direct averaging is not possible because ∆E values at different
∆Ei(t)
∑t ∆E(t)
(7)
where ∆E(t) ) (∑i∆E2i (t))1/2. After such normalization, ∆Ei provides the contribution of the mode i in the total ∆E as shown in Figure 2. In the figure the modes are ordered in the increasing energy values. In the values of ∆E shown in the figure, we excluded the constant part that corresponds to the energy differences at the respective equilibrium geometries (see Table 1). Hence only deviations from the equilibrium ∆E values due to the vibration deformation of the molecule along normal modes were taken into account. Due to using the same force-field non-Coulombic parameters for the cation and the neutral molecules for both adenine and guanine, the ∆E curves for the two systems have almost the same shape. There is one exception: the end molecule from the 5′ side of the chain where the curves for the two system differ. We mentioned this “edge” effect before. Since such an effect plays no role in the charge transfer involving the whole chain, the curve corresponding to this end position in the chain was not included in the figure. Without the curves corresponding to the end chain positions, the ∆E curves have very similar shapes for adenine and guanine for all positions of the charge in the strand. For adenine, as well as for guanine, some modes give positive contributions and some give negative ones. The negative contributions are typically much bigger (for example, mode 31 for adenine and mode 13 for guanine) than the positive contributions. This results in ∆E for an isolated molecule being higher than that for a molecule in the strand. Since a hole jump between two nucleobases depends on the relation between the cation VDE and the VIPs of the two neighboring neutral molecules (relation 6) a negative contribution to VIP should enlarge the probability of the charge migration. At the same time a negative contribution to VDE should lower the probability. Usually only intense vibrations (resulting in large negative contributions to VIP) are needed to cross the energy gap between VIP and VDE and enable a charge hopping to occur in the thermal hopping model of the charge transfer. The normal-mode analysis showed that some of the modes cause deformations to the aromatic ring (these modes are outof-plane molecular motions) while other modes involve peripheral side groups (CH3, NH2) or side atoms (H8 and H3 for adenine and O for guanine) (Figure 1). Among modes that have large contributions to ∆E for both adenine and guanine, the modes involving vibrations of the NH2 group play a crucial role. In the case of N9-methyladenine, the modes 39-41 and 4748 are almost pure NH2-group vibrations while in mode 38 vibrations of the NH2 group and vibrations of the aromatic ring mix. It should be mentioned that mode 41 is one of the most intense modes of the isolated molecule. For N9-methylguanine modes 42, 43, and 49 also involve the NH2 group vibrations and play an important role. But, while for methylguanine the most contributing mode 13 corresponds to the symmetric NH2 group vibrations and mode 8 includes rotation of that group, the normal mode 31 for methyladenine does not include atoms of the amino group. This mode contains vibrations of the H3 and H8 atoms instead. We can conclude that the NH2 contribution is important for both the adenine and guanine molecules. The same result, although in a different form, was obtained before by other authors.20,21 In the present implementation of the model, we only
11088 J. Phys. Chem. B, Vol. 111, No. 37, 2007
Volobuyev et al.
Figure 4. ∆E distribution over the normal modes for adenine (top) and guanine (bottom) molecules.
consider individual nucleobases for their ability to attract excess holes or electrons and we do not take into account how the base interactions affect the hole and electron affinities of the bases. However, it is clear that those effects are important. This, for example, can be seen through the analysis of the influence of the interbase interactions on the NH2 group vibrations and on the ∆E values. 3.2. Dynamics Properties. Dynamics properties of the charge transfer through DNA strands have been investigated in this work using some simple model systems. Calculations of the charge-transfer time have been performed for the following strands: (AT)n and (GC)n (n ∈ [4, ‚‚‚, 10], see Table 3). For such strands the dependency ln(τ) vs ln(N) should be linear due to the relation
Figure 5. Total transfer time t as a function of strand length n in logarithmic scale for (AT)n (above) and (GC)n (below) strands.
TABLE 4: Transfer Rates (ps-1) for First, Second, and Last Hole Jump in (AT)n and (GC)n Strandsa n strand
jump
4
5
6
7
8
9
10
ATn
1f2 2f3 (n - 1) f n average 1f2 2f3 (n - 1) f n average
0.03 0.05 0.10 0.06 0.46 1.16 1.38 1.00
0.06 0.11 0.16 0.11 0.79 1.75 1.52 1.54
0.11 0.13 0.32 0.18 0.90 1.89 1.70 1.60
0.11 0.20 0.44 0.24 0.93 1.62 2.59 2.00
0.13 0.30 0.23 0.26 1.00 1.59 0.90 1.68
0.19 0.43 0.44 0.34 1.05 1.84 1.91 1.63
0.16 0.20 0.37 0.28 0.96 2.00 1.60 1.61
GCn
a
-1 η
k≈τ N
Averaged rates were calculated for transfer over the whole strands.
(8)
For mixed GC/AT strands, there is no universal relation such as (8), so we have not considered such systems at the initial stage of the calculations. For the homogeneous strands, (AT)n and (GC)n, we have first determined the dependency of the rate (or the charge-transfer time) on the strand length. In these calculations we used the following values of the thresholds for all stands: tf ) 0.017 au and tb ) ∞. For such thresholds only charge transfer in the forward direction can occur. By setting tb ) ∞, we assume that the backward charge jumps do not provide important contributions to the total charge-transfer rate. As the calculations showed, this works better for the (GC)n strands than for the (AT)n strands. For each system we performed 100 calculations at the temperature 300 K. Each calculation started with the excess charge being placed at the first nucleobase in the strand and stopped when the charge reached the last nucleobase. In the first calculations the initial atomic positions and velocities were randomly chosen, and in each of the following calculations the initial positions and velocities were taken from the last time step of the previous calculation. For each calculation we determined the time of the charge transfer between each two adjacent nucleobases and the total transfer time being the sum of all pair transfer times. The results averaged over all the calculations performed for each strand are shown in Figure 5. In both cases ((AT)n and (GC)n) the dependency is close to linear. But for strands with the number of bases equal to four (n ) 4), a deviation from the linear behavior appears. In the theoretical models considered so far, n ) 4 has been a border case between the pure tunneling charge-transfer mechanism and the thermal-hopping mechanism.
The deviation from linearity is less significant for the (AT)4 system than for the (GC)4 system where the tunneling plays a larger role. Due to the stochastic character of the thermal-hopping charge transfer in the molecular dynamics simulations, deviations from linearity may increase with increasing n. However, as the calculations show, this does not alter the overall linear behavior of the rate. For example, when one nucleobase pair was added to the strand (AT)10 the slope of the regression line remained virtually unchanged. Estimations for the parameter η in (8) derived from the results of the calculations are 2.8 for the (AT)n strands and 1.6 for the (GC)n strands. These values seem somewhat higher than they should be. However, if we exclude the point corresponding to the smallest n value (n ) 4), the values decrease to 2.5 and 1.2, respectively. So, the most likely values of η for the (AT)n and (GC)n strands should be around 2.0 and 1.0. These values are consistent with what one can expect for the unbiased diffusive hopping ((AT)n) and the acceptor-direction-biased random walk ((GC)n). The hole transfer rates for all the strands were calculated from the transfer times and are shown in Table 4. The values of the averaged rates are in agreement with the available data.22 However, a direct comparison cannot be made because at this stage we performed the calculations only for homogeneous strands ((AT)n and (GC)n) while the experimental data correspond to mixed strands. Large discrepancies in the values of the constants are likely due to the stochastic character of the charge transfer. The discrepancy decreases slowly with the increasing number of the base pair in the strand. For both (AT)n and (GC)n the values of
Hole Transfer Rates in DNA Strands all constants are already converged starting from n ) 8. This is particularly the case for the average rates. These are the rates that can be determined from the experiment. As this is the case in the experiment,7 we need to know four different constants, KAA, KGG, KAG, and KGA, to predict the charge-transfer rate in a strand consisting of A’s and G’s. The first two constants, KAA and KGG, can be estimated using the results of the simulations performed in this work. The other two constants, KAG and KGA, can be calculated (less precisely) from the modeling of mixed adenine/guanine strands and this will be the subject of the next work. 4. Conclusions Some new features have been implemented in our model of the charge transfer through DNA strands. Since a charge can migrate through a DNA strand by either a thermal-hopping or a tunneling mechanism, both mechanisms are now implemented in our model. A hole can move by an adenine/adenine, guanine/ guanine, or adenine/guanine hopping or tunneling, but the preferred location of a hole is a guanine molecule, i.e., a hole at guanine is more stable than that at adenine. An excess charge (both negative and positive) can migrate in a DNA strand in either forward or backward direction even if a potential bias is applied to the strand. This is due to “kinetic” and “thermodynamic” factors that control the charge migration. The following conclusions can be drawn based on the results of the simulations performed in this work: The higher hole affinity of guanine compared to adenine makes a sequence of three guanines in the strand very difficult for a hole to cross. Consequently, any GGG fragment in the strand always plays the role of a hole trap and effectively stops the hole transfer. This observation is consistent with both experimental and previous theoretical reports. Among the vibrational normal modes of adenine and guanine, those involving the NH2 group vibrations play the most important role in thermal variability of the VIP and the VDE and, hence, they affect the thermally activated charge-transfer process the most. The charge migration in DNA strands can be characterized as either an unbiased diffusive hopping (mostly operating in(AT)n) or the acceptor-direction-biased random walk (mostly operating in(GC)n). The calculated rate constants are in good agreement with available experimental data. We will use those rates in future work to predict charge-transfer rates in mixed AT/GC strands. Supporting Information Available: A figure showing tunneling penetration through a potential barrier between two
J. Phys. Chem. B, Vol. 111, No. 37, 2007 11089 adenine molecules. This material is available free of charge via the Internet at http://pubs.acs.org. References and Notes (1) Gu, J.; Xie, Y.; Schaefer, H. J. Phys. Chem. B 2005, 109, 13067. (2) Voityuk, A. J. Phys. Chem. B 2005, 109, 10793. (3) Yoo, K.-H.; Ha, D. H.; Lee, J.-O.; Park, J.; Kim, J.; Kim, J. J.; Lee, H.-Y.; Kawai, T.; Choi, H. Phys. ReV. Lett. 2001, 87, 198102. (4) Jortner, J.; Bixon, M.; Langenbacher, T.; Michel-Beyerle, M. E. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 12759. (5) Smith, D. M. A.; Adamowicz, L. J. Phys. Chem. B 2001, 105, 9345. (6) Volobuyev, M.; Adamowicz, L. J. Phys. Chem. B 2005, 109, 1048. (7) Giese, B.; Biland, A. Chem. Commun. 2002, 667. (8) Giese, B. Acc. Chem. Res. 2000, 33, 631. (9) Case, D. A.; Pearlman, D. A.; Caldwell, J. W.; Cheatham, T. E., III; Ross, W. S.; Simmerling, C. L.; Darden, T. A.; Merz, K. M.; Stanton, R. V.; Cheng, A. L.; Vincent, J. J.; Crowley, M.; Tsui, V.; Radmer, R. J.; Duan, Y.; Pitera, J.; Massova, I.; Seibel, G. L.; Singh, U. C.; Weiner P. K.; Kollman, P. A. AMBER 6, University of California, San Francisco, 1999. (10) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. J. Am. Chem. Soc. 1995, 117, 5179. (11) Weisstein, Eric W. Chebyshev Polynomial of the First Kind. From MathWorld-A Wolfram Web Resource. http://mathworld.wolfram.com/ ChebyshevPolynomialoftheFirstKind.html. (12) Frish, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Zakrzewski, V. G.; Montgomery, J. A.; Stratmann, R. E.; Burant, J. C.; Dapprich, S.; Millam, J. M.; Daniels, A. D.; Kudin, K. N.; Strain, M. C.; Farkas, O.; Tomasi, J.; Barone, V.; Cossi, M.; Cammi, R.; Mennucci, B.; Pomelli, C.; Adamo, C.; Clifford, S.; Ochterski, J.; Petersson, G. A.; Ayala, P. Y.; Cui, Q.; Morokuma, K.; Malick, D. K.; Rabuck, A. D.; Raghavachari, K.; Foresman, J. B.; Cioslowski, J.; Ortiz, J. V.; Stefanov, B. B.; Liu, G.; Liashenko, A.; Piskorz, P.; Komaromi, I.; Gomperts, R.; Martin, R. L.; Fox, D. J.; Keith, T.; Al-Laham, M. A.; Peng, C. Y.; Nanayakkara, A.; Gonzalez, C.; Challacombe, M.; Gill, P. M. W.; Johnson, B. G.; Chen, W.; Wong, M. W.; Andres, J. L.; Head-Gordon, M.; Replogle, E. S.; Pople, J. A. Gaussian 98; Gaussian Inc.: Pittsburgh, PA, 1998. (13) Dupradeau, F.-Y.; Pigache, A.; Zaffran, T.; Cieplak, P. Automatic, highly reproducible and effective RESP and ESP charge Derivation: Implementation of multiconformation and/or multiorientation RESP fit in the R.E.D. program. Submitted for publication. (14) Pilar, F. Elementary quantum chemistry; Dover Publications Inc.: Mineola, NY, 1990. (15) Furrer, E.; Giese, B. HelV. Chim. Acta 2003, 86, 3623. (16) Giese, B. Curr. Opin. Chem. Biol. 2002, 6, 612. (17) Hush, N.; Cheung, A. Chem. Phys. Lett. 1975, 34, 11. (18) Dolgunicheva, O.; Zakrzewski, V. G.; Ortiz, J. V. J. Am. Chem. Soc. 2000, 122, 12304. (19) Wetmore, S. D.; Boyd, R. J.; Eriksson, L. A. Chem. Phys. Lett. 2000, 322, 129. (20) Alexandre, S.; Artacho, E.; Soler, J.; Chacham, H. Phys. ReV. Lett. 2003, 91, 108105. (21) Danilov, V.; Anisimov, V. Biomol. Struct. Dyn. 2005, 22, 471. (22) Lewis, F.; Letsinger, R.; Wasielewsky, M. Acc. Chem. Res. 2001, 34, 159.
