Charge transport in a nonlinear, three-dimensional DNA model with disorder

CHARGE TRANSPORT IN A NONLINEAR, THREE–DIMENSIONAL DNA MODEL WITH DISORDER JFR ARCHILLA Nonlinear Physics Group of the University of Sevilla, Departamento de F´ısica Aplicada I, ETSI Inform´atica, Avda Reina Mercedes s/n, 41012– Sevilla, Spain Email: archil[email protected] D HENNIG AND J AGARWAL Freie Universit¨at Berlin, Fachbereich Physik, Institut f¨ur Theoretische Physik, Arnimallee 14, 14195–Berlin, Germany We study the transport of charge due to polarons in a model of DNA which takes in account its 3D structure and the coupling of the electron wave function with the H–bond distortions and the twist motions of the base pairs. Perturbations of the ground states lead to moving polarons which travel long distances. The influence of parametric and structural disorder, due to the impact of the ambient, is considered, showing that the moving polarons survive to a certain degree of disorder. Comparison of the linear and tail analysis and the numerical results makes possible to obtain further information on the moving polaron properties. 1. Introduction Charge transport along DNA is a subject of particular interest for two main reasons: on the one hand, it plays a fundamental role in biological functions as repair and biosynthesis; on the other hand, because of possible applications in molecular electronics and as molecular wires 1. Results on experimental DNA conductivity are controversial. It has been reported that it is good conductor 2, insulator 3and semiconductor 4. A possible explanation, among others, for these striking differences can be the different long-range correlations of the DNA sequences 5. In this work we propose a nonlinear mechanism for charge transport along DNA in the framework of the base pair picture 6,7, taking into account its spatial structure and the coupling of the spatial and electron variables. It turns out to be an efficient mechanism for charge transport which survives to a certain degree of diagonal and structural disorder. 1 2 2. Model We consider a semi–classical three–dimensional, tight–binding model for the DNA molecule, the sketch of the system being shown in Fig. 1. The lattice rn l0+dn n−1 θ0+θn n−1 θ0 R0 l0 Figure 1. Sketch of the helical structure of the DNA model, the bases being represented by bullets. Geometrical parameters R0,θ0,l0and the radial and angular variables rn and θn,n−1, respectively are indicated. oscillators are treated classically while the charge is described by a quantum system. The justification is that the nucleotides are large molecules with mass numbers of about 300. This also implies that molecular motions are small and slow compared to the one of the charge particle. The Hamiltonian of the whole system is given by ˆ H=ˆ Hel +ˆ Hrad +ˆ Htwist , where ˆ Hrad and ˆ Htwist are, de facto, classical Hamiltonians and we can omit the hat on them. The Hamiltonians corresponding to the distances between nucleotides in each base pair Hr, and the twist motion Htwist are given by Hrad =X n (pr n)2 2M+MΩ2 rr2 n 2;Htwist =X n¡pθ n,n−1¢2 2J+JΩ2 θθ2 n,n−1 2(1) where {rn}represent the stretchings from the equilibrium distance, Mthe reduced mass or each base pair, Ωrthe linear radial frequency, {θn,n−1} are the angles between two consecutive base pairs with respect to their equilibrium value θ0= 36◦,Jis the inertia moment, and Ωθ, the twist linear frequency. 3 The electronic part is given by a tight–binding system of the form ˆ Hel =X n En|nihn| − Vn−1,n|n−1ihn| − Vn+1,n|n+ 1ihn|.(2) In the state |nithe charge carrier is localized at the nth base pair. The quantities Vn−1,n represent the nearest–neighbour transfer integrals along the base pairs and Enare the on–site matrix–elements. We write the electronic state as |Ψi=Pncn(t)|ni, where cn(t) is the probability amplitude of finding the electron at the state |ni. The interaction between the electronic variables and the structure variables rnand θn,n−1arises from the dependence of the electronic parameters Enand Vn,n−1on the spatial coordinates, given by En=E0 n+k rnand Vn,n−1=V0(1 −α dn,n−1), where dn,n−1is the first order Taylor expansion of the distances between consecutive nucleotides, with respect to their equilibrium value. Is is given by dn,n−1≃R0 l0 [(1 −cos(θ0)) (rn+rn−1) + sin(θ0)R0θn,n−1].(3) Realistic parameters for DNA molecules are given in Refs. 8,9,10. We scale the time according to t→Ωrt, which allow us to write the Hamiltonian in terms of dymensionless quantities 10. The classical Hamiltonian is defined as Hclass =hΨ|ˆ H|Ψi, the dynamical equations for rnand θn,n−1are ˙pr n=M¨rn=−∂Hclass/∂rnand ˙pθ n,n−1=J¨ θn,n−1=−∂Hclass/∂θn,n−1, while the evolution equations for the electron variables cnare obtained from the Schr¨odinger equation i~(∂Ψ/∂t) = ˆ Hel|Ψi, which is equivalent to ˙cn=−(i/~)(∂Hclass/∂c∗ n). In this way we can obtain the scaled dynamical equations, which are iτ˙cn= (E0 n+k rn)cn −(1 −α dn+1,n)cn+1 −(1 −α dn,n−1)cn−1(4) ¨rn=−rn−k|cn|2−αR0 l0 (1 −cos θ0) ש[c∗ n+1cn+cn+1c∗ n] + [c∗ ncn−1+cnc∗ n−1]ª(5) ¨ θn,n−1=−Ω2θn,n−1−α V R2 0 l0 sin θ0[c∗ ncn−1+cnc∗ n−1],(6) where τ=~Ωr/V0represents the time scale separation between the fast electron motion and the slow bond vibrations. The values of the scaled parameters are 10 τ= 0.2589, Ω2= [0.709 − 1.417]×10−2,V= 0.0823, R0= 34.862 and l0= 24.590, the time unit being ∼1.6 ps. Note that the time scales for the different variables differ in an order of magnitude: the fastest variables are the {cn}, with a characteristic 4 frequency of order 1/τ ∼4, followed by {rn}with the unity, and {θn,n−1} with Ω ∼0.08 There are no reliable values for the electron–radial and electron–twist coupling parameters, kand α. The criterion we have taken for their values is the consistency of the numerical simulations with the hypothesis, i.e., the deformations of the helix are small, and the linear approximations of the distance dn,n−1and the trigonometric functions remain valid. Typical values are k∼1 and α∼0.002. 3. Stationary states We suppose initially that rnand θn,n−1are constant, i.e., the BornOppenheimer approximation. This allows the obtention of expressions for them from Eqs. (5)–(6) which are inserted in Eq. (4), leading to a nonlinear Schr¨odinger equation for the electronic amplitudes: iτ˙cn=·E0 n−k2|cn|2−k α R0 l0 (1 −cos θ0) ש[c∗ n+1cn+cn+1c∗ n]+[c∗ ncn−1+cnc∗ n−1]ª¤cn −( 1 −α dn+1 n)cn+1 −( 1 −α dn,n−1)cn−1,(7) whith {dn,n−1}given by Eq. (3), depending algebraically on the {cn} through {rn}and {θn,n−1}. To obtain stationary localized solutions we use a numerical method 11,12. We substitute cn= Φnexp(−iE t/τ) in Eq. (7), with Φnconstants, and obtain a nonlinear difference system EΦ = ˆ A φ, with Φ = (Φ1, . . . , ΦN). Its solutions are attractors of the map: Φ→Φ0=ˆ AΦ/kˆ AΦk,(8) k·kbeing the quadratic norm. We start with a completely localized state at a site n0, i.e., Φn=δn,n0and apply the map above until Φ0= Φ. In this way we obtain stationary localized solutions and their energy E. In the ordered case, with E0 n=E0,∀n,E0can be made zero with the ansatz cn→cnexp(−iE0t/τ). We obtain symmetric polarons with width of about 20 sites, combined with a local compression of order ∼0.15 and a local unwinding of ∼1◦and energies E.−2∼ −0.2 eV. Considering static diagonal disorder with random E0 nthe localization is enhanced, due to Anderson localization, the polaron being asymmetric, with a specific shape that depends on each particular disorder implementation. Similar results are obtained with structural disorder, i.e., random equilibrium values R0 nand θ0 n. The Floquet analysis shows that these polarons are linearly stable 10. 5 4. Charge transport To activate the polaron motion we need to perturb the zero velocities of the ground state with a localized, spatially antisymmetric mode. This is obtained either with the time–consuming pinning–mode method 13, either with the simpler discrete gradient method, i.e., perturbing the variables rnwith velocities parallel to {rn+1 −rn−1}14. The main difference is that the first guarantees the mobility, while the second does not, but in practice works well very often. Values of the modulus of the kick velocity λr= (P˙r2 n(0))1/2∼0.02, equivalent to a kinetic energy of ∼200 meV are appropriate to obtain good mobility with low radiation. We obtain propagation of the localized electron amplitudes with constant velocity. A complementary compression and unwinding of the helix travels with it, while a localized oscillation of the angular variables remains at the initial positions. If we increase the diagonal disorder with En∈[−∆E, ∆E], the charge transport still takes place, until values of ∆E∼0.05. Thereafter mobility becomes impossible and the polaron is pinned to the lattice. With disorder, however, the velocity of the first momentum of the electronic occupation amplitude nc(t) = Pnn|cn(t)|2is not uniform, as shown in Fig. 2. We also consider structural disorder with random distributions of the base pair spacings and angles around their ordered case values, with different mean standard deviations ∆. The conclusion is that, up to values of ∆ = 10%, it does not significantly affect polaron mobility, although there is a slight reduction of its velocity. It is interesting to remark that for the values of αcoherent with our small amplitude hypotheses it is not possible to move the polaron by kicking the angular variables. These results are coherent with the fact that the Floquet analysis shows that the pinning mode appears only in the radial variables and its eigenvalue separates from the optimal value around 1 when the disorder is increased 10. 5. Linear and tail analysis Although our system is studied in the nonlinear regime, and the results obtained throughout this work are essentially nonlinear, there is a number of linear techniques that can be applied as a reference and as a tool for obtaining useful information. They are the study of the linear system itself and the tail analysis. 6 0 100 200 300 400 500 600 700 800 900 1000 −80 −70 −60 −50 −40 −30 −20 −10 0 10 20 30 t nc Figure 2. Diagonal disorder. The position of the center of the electron breather as a function of time and for different amounts of disorder. Full line (Ordered case), short dashed line (∆E= 0.025), dashed-dotted line (∆E= 0.050), dotted line (∆E= 0.100) and long dashed line (∆E= 0.500). 5.1. The linear system If in Eqs. (4)–(6) we cancel out the nonlinear terms, we obtain: iτ˙cn=−cn+1 −cn−1,¨rn=−rn,¨ θn n−1=−Ω2θn n−1.(9) In the first equation the terms En, which are all equal to some value E0 in an homogeneous system, have been eliminated through the ansatz cn→ cnexp(−iE0t/τ). The system becomes decoupled and the variables rn and θn n−1correspond to independent linear oscillators with frequencies wr= 1 and Ω = 0.0842, respectively, in the scaled units. To obtain the ground state we substitute in Eq. (9) cn(t) = φnexp(−iE t/τ), φnbeing time independent. We obtain the stationary discrete Schr¨odinger equation: E φn=−φn+1 −φn−1. Substitution of the linear modes φn= exp(i q n) leads to: E=−2 cos q. Therefore, the linear energy spectrum runs from −2 to 2. The minimum energy, E=−2, corresponds to the wave vector q= 0, i.e., all the oscillators vibrating in phase. The nonlinear ground states described above have E.−2 and derive from this mode. 7 5.2. Tail analysis For sites far enough from the polaron center, which can be only a few sites, we can still apply the linear analysis. Now, the substitution of the tail mode cn(t) = φnexp(−ξ n −iq n), with n > 0 and ξ > 0, leads to: E=−2 cos(q) cosh(ξ) and 2 sin(q) sinh(ξ) = 0. The second equation, implies that only two wave vectors are possible, q= 0, with negative energy, and q=π, with positive one. The energy of the first, which is the one that appears throughout this paper because it is movable, is given by: E=−2 cos(ξ)<−2. Therefore, the nonlinearity, produces localized states with larger values of the frequency wp=|E|/τ outside the linear spectrum. The distance of the energy calculated in the full system from Eq=0 =−2, is a measure of the degree of nonlinearity. It is interesting to relate the values of the energies with the polaron breadth. Here, we define it loosely as three times the number of sites necessary for |φn|2to decrease to a 5%. The factor 3 is there to allow for both and for the nonlinear center of the polaron. This leads to ∆n∼4.5/ξ. An approximate table is: ξ2 1 0.5 0.15 0.075 ∆n2.25 4.5 9 30 60 E -7 -3 -2.25 -2.02 -2.006 where the dimensionless Eunits are equivalent to 0.1 eV. This makeshift method gives a good estimate of the breadth of the excitation and produces a good fitting a few units far from its center. The numerically calculated polarons found in this work correspond to ∆n∼20. The tail analysis can also be applied to the moving polaron, far enough from its center. Let us propose a localized traveling wave of the form: cn= exp{−ξ(n+v t)−i (q n+wpt)}, which represents a traveling localized wave moving to the left behind the polaron, with positive ξ,q, and velocity v. Substitution in Eq. (9) leads to the following equations: E≡wpτ= −2 cos(q) cosh(ξ) and τ ξ v = 2 sin(q) sinh(ξ). They are undetermined, as there are four unknowns: E,q,ξand v, and do not allow us to determine the velocity. This is little surprise, as the polaron velocity depends on the energy given with the perturbation to the static one, and this is done to the spatial variables. However, these equations can be used to obtain an estimate of the moving polaron energy and wave number, using the decay length ξcorresponding to the static one, as there is no appreciable change of shape, and the velocity observed within the simulations. The conclusion is that there is only a negligible increase of the energy, i.e., the energy given by the kick is stored in the spatial variables. The fact that there is 8 a imaginary part of the energy −i ξ v τ is, as it is known, consequence of the fact that we are dealing with only a part of the system 15, in which the energy is decreasing as the polaron moves to the left. 6. Conclusions We have considered a model for charge transport along DNA, taking into account the 3D structure. We obtain linearly stable, localized stationary states. We are able to move them using the pinning mode and the discrete gradient methods. For the ordered case, the translational velocity is constant, while for the disordered case the electron motion is irregular and if the disorder is high enough it is impossible. We have performed a linear and tail analysis of the system, which makes possible to obtain the phonon spectrum, measure the degree of nonlinearity, obtain estimates of the polaron breath and energy. It shows that the charge transports very little energy, being most of it associate with the spatial coordinates. The main conclusion of the whole work can be that our proposed mechanism is an efficient means for charge transport along ordered and disordered DNA. Acknowledgments The authors are grateful to partial support under the LOCNET EU network HPRN-CT-1999-00163. JFRA acknowledges DH and the Institut f¨ur Theoretische Physik for their warm hospitality References 1. M Ratner. Nature, 397:480, 1999. 2. B Alavi P Tran and G Gruner. Phys. Rev. Lett., 85:1564, 2000. 3. E Braun, Y Eichen, U Sivan, and G Ben-Yoseph. Nature, 391:775, 1998. 4. D Porath, A Bezryadin, S de Vries, and C Dekker. Nature, 403:635, 2000. 5. P Carpena, P Bernaola-Galvan, PC Ivanov, and HE Stanley. Nature, 418:955, 2002. 6. S Cocco and R Monasson. J. Chem. Phys., 112:10017, 2000. 7. M Barbi, S Cocco, and M Peyrard. Phys. Lett. A, 253:358, 1999. 8. M Peyrard and MD Kruskal. Physica D, 14:88, 1984. 9. L Stryer. Biochemistry. Freeman, New York, 1995. 10. D Hennig, JFR Archilla, and J Agarwal. Physica D, in press, 2003. 11. G Kalosakas, S Aubry, and GP Tsironis. Phys. Rev. B, 58:3094, 1998. 12. NK Voulgarakis and GP Tsironis. Phys. Rev. B, 63:14302, 2001. 13. D Chen, S Aubry, and GP Tsironis. Phys. Rev. Lett., 112:139, 1996. 14. M Iba˜nes, JM Sancho, and GP Tsironis. 65:041902–041914, 2002. 15. C Cohen-Tannoudji, B Du, and F Lalo´e. Quantum Mechanics, volume 1. John Wiley & Sons Inc, 1977.