Charge transport in poly(dG)-poly(dC) and poly(dA)-poly(dT) DNA polymers

arXiv:nlin/0308003v1 [nlin.PS] 1 Aug 2003 Charge transport in poly(dG)-poly(dC) and poly(dA)-poly(dT) DNA polymers D. Hennig1, E.B. Starikov2, J.F.R. Archilla3and F. Palmero3 1Freie Universit¨at Berlin, Fachbereich Physik, Institut f¨ur Theoretische Physik Arnimallee 14, 14195 Berlin, Germany 2Karolinska Institute, Center for Structural Biochemistry NOVUM S - 14157 Huddinge, Sweden 3Group of Nonlinear Physics, Departamento de F´ isica Aplicada I ETSI Inform´atica, Avda Reina Mercedes, s/n. 41012 - Sevilla, Spain August 1, 2003 Abstract We investigate the charge transport in synthetic DNA polymers built up from single types of base pairs. In the context of a polaron-like model, for which an electronic tight-binding system and bond vibrations of the double helix are coupled, we present estimates for the electron-vibration coupling strengths utilizing a quantum-chemical procedure. Subsequent studies concerning the mobility of polaron solutions, representing the state of a localized charge in unison with its associated helix deformation, show that the system for poly(dG)-poly(dC) and poly(dA)-poly(dT) DNA polymers, respectively possess quantitatively distinct transport properties. While the former supports unidirectionally moving electron breathers attributed to highly efficient long-range conductivity the breather mobility in the latter case is comparatively restrained inhibiting charge transport. Our results are in agreement with recent experimental results demonstrating that poly(dG)-poly(dC) DNA molecules acts as a semiconducting nanowire and exhibits better conductance than poly(dA)-poly(dT) ones. PACS numbers: 87.-15.v, 63.20.Kr, 63.20.Ry 1 Introduction Particularly with view to possible applications in molecular electronics based on biomaterials electronic transport (ET) through DNA has recently become of intensified interest [1],[2]. Although the debate whether DNA constitutes a conductor is still ongoing, there exist already strong experimental evidence that DNA forms an effectively one-dimensional molecular wire [3]. Among several theoretical attempts to describe the charge transport mechanism in DNA the polaron approach has turned out lately to be a promising candidate for modeling constructive interplay between the charge carrying system and vibrational degrees of freedom of the DNA conspiring to establish coherent ET [4]-[8]. Recent 1 experiments are in support of the polaron mechanism for ET in DNA polymers [9]. The present study deals with the theoretical description of ET in synthetically produced DNA polymers consisting of a single type of base pairs, i.e. either poly(dG)-poly(dC) or poly(dA)-poly(dT) DNA polymers [10]. Utilizing a nonlinear approach based on the concept of breather and polaron solutions we explore whether conductivity depends on the type of the DNA polymer which might be of interest for the design of synthetic molecular wires. Moreover, we ameliorate preceding studies of ET in DNA [11],[12] in the sense that, instead of adjusting the coupling parameters, we use now credible estimates for them derived with the help of quantum-chemical methods. 2 Model for polaron-like charge transport in DNA Our model of charge transport in DNA is based on the finding that the charge migration process is dominantly influenced by the transverse vibrations of the bases relative to each other in radial direction within a base pair plane [13]. In fact, the impact of other vibrational degrees of freedom (e.g. twist motions, helical pitch changes and longitudinal acoustic phonons along the strands which are significantly restrained by the backbone rigidity) can be expected to be negligible with respect to ET in DNA. Hence, the motion can be viewed as confined to the base pair planes [14]. The Hamiltonian for the ET along a strand in DNA comprises two parts H=Hel +Hvib ,(1) with Hel is the part which describes the ET over the base pairs and Hvib represent the dynamics of radial vibrations of the base pairs . The electronic part is given by a tight-binding system Hel =X n En|cn|2−Vn n−1c∗ ncn−1+cnc∗ n−1.(2) The index ndenotes the site of the n−th base on a strand and |cn|2determines the probability to find the electron (charge) residing at this site. Enis the local electronic energy and Vnn−1is the transfer matrix element mediating the transport of the electron along the stacked base pairs. We make the usual assumption that ET takes place only along the base pair sequence on a strand excluding inter-strand ET. The vibronic part of the Hamiltonian Hvib models dynamical changes of the radial equilibrium positions of the bases. Supposing that these radial vibrations can be treated classically and harmonically we represent Hvib as Hvib =1 2X n1 M(pr n)2+MΩ2 rr2 n.(3) The radial coordinates rnquantify the radial displacements of the base units from their equilibrium positions along the line bridging two bases of a base pair within the base pair plane. Mdenotes the reduced mass and Ωris the harmonic frequency. Due to the constraint enforced by the sugar-phosphate backbone radial displacements lead not only to changes of the equilibrium lengths of the 2 rn l0+dn n−1 θ0 θ0 R0 l0 Figure 1: Sketch of the structure of the DNA model. The bases are represented by bullits and the geometrical parameters R0,l0,θ0,rnand dn n−1are indicated. hydrogen bonds (bridging two bases of a base pair) but are also connected with a change of the three-dimensional distance between two consecutive bases on a strand entailing deformations of the corresponding covalent bonds. In an expansion up to first order around the equilibrium positions [15] this stackingdistance is determined by dn n−1=R0 l0 ( 1 −cos θ0) (rn+rn−1).(4) The geometrical parameters R0,θ0and l0determine the equilibrium value of the radius, the twist angle between two adjacent base pairs, and the equilibrium length of the covalent bond linking two consecutive bases along a strand, respectively. The latter is given by l0=qa2+ 4R2 0sin2(θ0/2) ,(5) with abeing the distance between neighboring base pairs measured along the orientation of the helix axis. A sketch of the structure of the DNA model and the designation of the geometrical parameters R0,l0,θ0,rnand dn n−1is presented in Figure 1. With respect to the interaction between the electronic and the vibrational degrees of freedom variable, cnand rn, respectively, it is assumed that the electronic parameters Enand Vn n−1are modified by displacements of the bases within the base pairs. As a quantum-chemical computation of the geometry dependence of the electronic parameters Enand Vn n−1reveals their most significant modulation originates from radial distortions of the helix which are related with hydrogen and covalent bond deformations, respectively. However, the influence of small angle deformations on the values of the electronic parameters can be discarded (see also further below). 3 The modulation of the on-site electronic energy E0by the radial vibrations of the base pairs is expressed as En=E0+k rn.(6) On the other hand the actual charge occupation has its impact on the local radial distortion of the helix. Furthermore, the transfer matrix elements Vn n−1 are assumed to depend on the three-dimensional stacking-distance between two consecutive bases on a strand as follows Vn n−1=V0(1 −α dn n−1).(7) The quantity αregulates how strong Vn n−1is influenced by the distance. As typical parameters for DNA molecules one finds [13],[15]: a= 3.4˚ A, R0≈10˚ A,θ0= 36◦,E0= 0.1eV , Ωr= 6.252 ×1012 s−1,V0≃0.1eV and M= 4.982 ×10−25kg. After scaling the time as t→Ωrtone passes to the dimensionless quantities: ˜rn=sMΩ2 r V0 rn,˜ k=k pMΩ2 rV0 ,˜ E0=E0 V0 ,(8) ˜α=sV0 MΩ2 r α , ˜ R0=sMΩ2 r V0 R0,(9) and for the sake of convenience the tildes are omitted in the following. The values of the scaled parameters are obtained as R0= 34.862 and l0= 24.590. The set of coupled equations of motion read as i τ ˙cn= (E0+k rn)cn −(1 −α dn+1,n)cn+1 −(1 −α dn n−1)cn−1(10) ¨rn=−rn−k|cn|2−R0 l0 (1 −cos θ0) ×α[c∗ n+1cn+cn+1c∗ n] + [c∗ ncn−1+cnc∗ n−1](11) and the ratio τ=~Ωr/V0= 0.0411 determines the time scale separation between the slow electron motion and the fast bond vibrations. (Notice that any E0cnterm on the r.h.s. of Eq. (10) can be eliminated by a gauge transformation cn→exp(−iE0t/τ)cn.) Contrary to previous studies we use for our computations credible values for the electron-mode coupling strengths kand αas a result of our a quantumchemical computational procedure. In this work, we perform quantum-chemical calculations on symmetrical homodimers consisting of two nucleoside Watson-Crick base pairs (adenosinethymidine (AT) and guanosine-cytidine (GC) base pair steps (BPS)) stacked over each other to mimic the conventional A-DNA and B-DNA conformations (for a schematic view see [16]). Taking into account DNA backbone at least in form of intact sugar moieties, instead of substituting it by protons or methyl groups, is necessary for the consistency of the calculations [17] and to correctly describe charge transfer through DNA duplexes [18]. We use semiempirical all-valence-electron PM3 Hamiltonian [19] within the MOPAC7 version of CI 4 (configuration interaction) approximation, as described in detail in [20]. We chose here PM3-CI method, but not an ab initio one, since a) The molecular fragments involved are very large; b) Similar ab initio calculations even for smaller segments experience difficulties with SCF convergence [21]; c) The work [17] used similar semiempirical Hamiltonian (AM1) for analogous molecular fragments; d) PM3-CI approximation is good at describing excited states of nucleoside base pair steps [20]; e) CI approximation is indispensable when charged states of nucleic acid bases are considered using semiempirical quantum chemistry (cf., e.g., [22] and references therein). Here the effects of small, but non-negligible (up to 0.1 ), radial stretching and compressing of Watson-Crick hydrogen bonds (WC H-bonds) on the MO energies of the BPS under study have been estimated. Specifically, we perturbed in the above sense the equilibrium lengths of the WC H-bonds in only one of the base pairs in all the BPS involved and monitored the resulting changes in the energies of the highest occupied molecular orbital and of the occupied molecular orbital next to the highest one (HOMO and HOMO-1, respectively), as compared with the equilibrium values of these energies. According to the Koopmans theorem well known in quantum chemistry (cf., e.g., [21] and references therein), the HOMO energy is approximately equal to molecular ionisation potential and in the tight-binding approximation could be viewed as the site energy, whereas the difference between the HOMO and HOMO-1 energies is approximately twice the hopping integral. With this in mind, we were able to estimate kand αparameters of our model by calculating linear regression of the corresponding site energy and hopping integral changes, respectively, onto WC H-bond distance perturbations. We were able to arrive at a very tight linear correlation between the former and the latter ones. Interestingly, we failed to reveal any correlation between the tight-binding Hamiltonian parameter changes and BPS twist angle/helical pitch perturbations (’helical pitch’ is the distance ain Eq. (5)). The perturbations in the two latter parameters were also relatively small (up to 10 degrees and 0.1 , respectively). For the coupling parameters of the poly(dA)- poly(dT) DNA polymer we obtain: k= 0.0778917 eV/ ˚ Aand α= 0.053835 ˚ A−1. The corresponding values for the poly(dG)-poly(dC) DNA polymer are determined as k=−0.090325 eV/ ˚ Aand α= 0.383333 ˚ A−1. That we are equipped with the quantum-chemical estimates of the coupling parameters αand kis a definite step forward and distinguishes the present study from previous ones [11],[12]. We emphasize that the quantum-chemical estimates for the coupling parameters differ significantly from the ones used for the model study in [11],[12] where these parameters, in lack of reliable values for them, were treated as adjustable. As a consequence there is a difference between the ET scenario described in [11],[12] and the one we are going to illustrate in the following. 3 Stationary localized electron-vibron states Caused by the nonlinear interplay between the electronic and the vibrational degrees of freedom of the helix the formation of polaronic electron-vibration 5 Figure 2: The spatial pattern of the polaronic electron-vibration compound. (a) The electronic part. Full (dashed) line: poly(dA)-poly(dT) (poly(dG)- poly(dC)) DNA polymer. (b) Unitless radial deformation pattern. Assignment of line types as in (a). compounds is possible [4]-[8],[11],[12]. We construct such localized stationary solutions of the coupled system (10),(11) with the help of the nonlinear map approach explained in detail in [23]. In Figure 2 we depict the profiles of the (standing) polaron states for the poly(dA)-poly(dT) and the poly(dG)-poly(dC) DNA polymer, respectively. In both cases the polarons are of fairly large extension (width). Regardless of the DNA polymer type the electronic wave function is localized at lattice site and the envelope of the amplitudes decays monotonically and exponentially with growing distance from this central site (base pair). However, the electronic wave function of the poly(dG)-poly(dC) DNA polymer is stronger localized than the one of its poly(dA)-poly(dT) counterpart. Accordingly, the attributed radial displacement patterns are exponentially localized at the central lattice site. Concerning the resulting static radial helix deformations we find that there is a drastic difference between the poly(dA)- poly(dT) and the poly(dG)-poly(dC) DNA polymers. In the former case the overall non-positive radial amplitudes imply that the H-bridges experience contractions. In contrast, in the latter case the H-bridges get stretched. Nevertheless, these deformations are rather weak, i.e. on the order of 1.5 10−3˚ A. 4 Charge transport We study now the ET achieved by moving polarons. In fact, since the constructed polaron solutions are of fairly large extension (the half-width involves .100 lattice sites) we can expect them to be mobile. In order to activate polaron motion we used the discrete gradient method [24] to obtain suitable initial perturbations of the momentum coordinates pr nwhich initiate coherent motion of the polaron compound. We consider first the case of the poly(dG)-poly(dC) DNA polymer. The propagation features are illustrated in Figure 3 where the spatio-temporal evolution of the electronic and the vibrational radial breather are shown. For a DNA lattice consisting of 300 sites (base pairs) the set of coupled equations (10),(11) was integrated using a fourth-order Runge-Kutta method and periodic boundary conditions were imposed. Maintenance of the norm conservation Pn|cn(t)|2= 1 served to assure accurate computations. 6 Figure 3: Breather motion along the DNA for the poly(dG)-poly(dC) DNA polymer.(a) The electronic breather. (b) The vibrational breather. The radial deformations rn(t) are given in dimensionless units. 7 Figure 4: Time evolution of the first momentum of the electronic occupation probability. Full (dashed) line: poly(dA)-poly(dT) (poly(dG)-poly(dC)) DNA polymer. The electronic component sets off to move directionally along the lattice in unison with the vibrational amplitude pattern. However, there remains a small amplitude vibrational breather at the initial position. In both cases the localized structures are practically preserved apart from the arising amplitude breathing indicating a periodic energy exchange between the electronic and vibrational degrees of freedom. Hence, long-range ET is achievable. With the initial injection of kinetic energy the (vibrational) system possesses now increased energy content so that we observe vibrational breathers with amplitudes being larger than those of their static equivalents (cf. Figures 2 (b) and 3 (b)). As the poly(dA)-poly(dT) DNA polymer is concerned we found that it does not exhibit such good conductivity as its poly(dG)-poly(dC) counterpart. The different propagation scenarios are properly illustrated using the time evolution of the first momentum of the electronic occupation probability defined as ¯n(t) = Pnn|cn(t)|2. One can clearly observe that in the(dG)-(dC) case electron propagation proceeds unrestrictedly and unidirectionally with uniform velocity (see Figure 4). Distinctly, the (dA)-(dT) electron moves itinerantly in a confined region, comprising not more than five bases, around the starting site restraining conductivity. Regarding the energy storing capacity we monitored the normalized participation number defined as p(t) = P(t) P(0) (12) with P(t) = 1 Pn|cn(t)|4.(13) Since the electronic wave function is normalized the electron breather is completely confined at a single site if p= 1 and is uniformly extended over the lattice if pis of the order N, viz. the number of lattice sites. Hence, pmeasures how many sites are excited to contribute to the electronic breather pattern. 8 Figure 5: The normalized participation number p(t). Full (dashed) line: poly(dA)-poly(dT) (poly(dG)-poly(dC)) DNA polymer. From Figure 5 we infer that the (dG)-(dC) electron breather extension performs oscillations the maximal amplitudes of which correspond to slight growth of the spatial width up to merely ≃1.2 times its starting value. In comparison, the width of the (dA)-(dT) electronic amplitude pattern experiences stronger broadening leading to a more extended electron state. Considering the energy exchange between the electronic and vibrational subsystems we observe that this process evolves with frequency Ωr. Hence, the dynamical distribution between the electronic and vibrational energy is driven by the radial harmonic vibrations. In particular, the amount of kinetic energy flowing from the vibrational system into the electronic subsystem (initially the electronic kinetic energy is zero) is governed by the strength of the coupling α. Apparently, in the (dG)-(dC) case the energy sharing is by far more pronounced than in the (dA)-(dT) case due to the fact that the coupling strengths differ by almost an order of magnitude, that is α(dG)−(dC)= 0.110 and α(dA)−(dT )= 0.0154. This explains the large content of kinetic energy injected from the vibrational degrees of freedom into the electronic ones inducing high electron mobility. Furthermore, the k’s have opposite sign, i.e. k(dA)−(dT )=−0.2591 and k(dG)−(dC)= 0.2234, which leads to the sign difference in the radial distortion patterns r(dA)−(dT ) nand r(dG)−(dC) n. An open question which is subject of further research is whether the polarons survive at ambient temperature. In spite of the values of the radial variables being small the polaron is a compound object and the fact that the transfer integral elements Vn n−1≈0.1 eV are larger than kBT≈0.025 eV at T= 300 K, suggests that they would survive. Equally the effect of an electric field during the whole movement of the polaron is currently investigated. In conclusion, we have found that conductivity in synthetically produced DNA molecules depends on the type of the single base pair of which the polymer is built of. While a polaron-like mechanism, relying on the nonlinear coupling between the electron amplitude and radial vibrations of the base pairs, is responsible for long-range and stable ET in (dG)-(dC) DNA polymers, the 9