Multiscale model of electronic behavior and localization in stretched dry DNA

TL;DR

This paper develops a multiscale model combining quantum calculations and tight-binding approaches to study how stretching DNA affects its electronic properties and electron localization, revealing significant localization changes.

Contribution

It introduces a multiscale modeling framework that links atomic-scale electronic states with mesoscopic structural distortions in stretched DNA.

Findings

Localization length decreases significantly in stretched DNA

Electronic states are sensitive to mechanical deformation

Model captures environmental effects on electron transport

Abstract

When the DNA double helix is subjected to external forces it can stretch elastically to elongations reaching 100% of its natural length. These distortions, imposed at the mesoscopic or macroscopic scales, have a dramatic effect on electronic properties at the atomic scale and on electrical transport along DNA. Accordingly, a multiscale approach is necessary to capture the electronic behavior of the stretched DNA helix. To construct such a model, we begin with accurate density-functional-theory calculations for electronic states in DNA bases and base pairs in various relative configurations encountered in the equilibrium and stretched forms. These results are complemented by semi-empirical quantum mechanical calculations for the states of a small size [18 base pair poly(CG)-poly(CG)] dry, neutral DNA sequence, using previously published models for stretched DNA. The calculated electronic states are then used to parametrize an effective tight-binding model that can describe electron hopping in the presence of environmental effects, such as the presence of stray water molecules on the backbone or structural features of the substrate. These effects introduce disorder in the model hamiltonian which leads to electron localization. The localization length is smaller by several orders of magnitude in stretched DNA relative to that in the unstretched structure.

Figures

Tables

Table 1. Table 1: Eigenvalues (in eV) of the frontier states for the DNA base pairs and the base-pair combinations, at the equilibrium configurations for the backbone distance, the axial distance (at zero relative angle of rotation) and the angle of rotation (at the equilibrium axial distance). The column labeled “min” gives the values of the distances and the angle at the equilibrium configurations. Due to symmetry the values for the minima at rotation angles larger than 180 o are similar to those given here and are not shown.

	min	HOMO	LUMO	gap
Backbone distance
AT	$8.67$ Å	$-1.63$	1.60	3.23
CG	$8.73$ Å	$-0.80$	1.31	2.11
Axial distance
AT-AT	$3.67$ Å	$-1.33$	1.37	2.70
CG-CG	$3.52$ Å	$-0.46$	0.95	1.41
AT-CG	$3.36$ Å	$-0.71$	1.00	1.71
Rotation angle
	${36}^o$	$-1.48$	1.58	3.06
AT-AT	${108}^o$	$-1.45$	1.68	3.13
	${180}^o$	$-1.55$	1.63	3.18
	${36}^o$	$-0.52$	1.22	1.74
CG-CG	${108}^o$	$-0.64$	1.54	2.18
	${180}^o$	$-0.94$	1.60	2.54
	${36}^o$	$-0.86$	1.51	2.37
CG-GC	${108}^o$	$-0.66$	1.43	2.09
	${180}^o$	$-0.60$	1.12	1.72
	${36}^o$	$-0.73$	1.38	2.11
AT-CG	${108}^o$	$-0.59$	1.27	1.86
	${180}^o$	$-0.81$	1.25	2.06

Table 2. Table 2: Parameters for the on-site ( ε 𝜀 \varepsilon ) and hopping matrix elements ( t i , i = 1 , 2 , 3 formulae-sequence subscript 𝑡 𝑖 𝑖 1 2 3 t_{i},i=1,2,3 ), for the HOMO and LUMO states of unstretched poly(CG)-poly(CG) DNA.

	HOMO	LUMO
$\epsilon$(eV)	3.12	$-0.09$
$t_1$(meV)	14.0	$-0.29$
$t_2$(meV)	2.60	0.04
$t_3$(meV)	0.09	0.26

Equations

$$|{\psi }^{(n)}⟩=\sum _{\nu }{c}_{\nu }^{(n)}|{\phi }_{\nu }⟩$$

$$|{\tilde{\psi }}^{(i)}⟩=\sum _n⟨{\psi }^{(n)}|{\tilde{\psi }}^{(i)}⟩|{\psi }^{(n)}⟩$$

$$\zeta =\sum _i\left(⟨{\tilde{\psi }}^{(i)}|{\hat{z}}^2|{\tilde{\psi }}^{(i)}⟩-⟨{\tilde{\psi }}^{(i)}|\hat{z}|{\tilde{\psi }}^{(i)}{⟩}^2\right)$$

$$⟨{\tilde{\psi }}^{(n)}|\hat{z}|{\tilde{\psi }}^{(m)}⟩(z_n-z_m)=0$$

$$z_n=⟨{\tilde{\psi }}^{(n)}|\hat{z}|{\tilde{\psi }}^{(n)}⟩.$$

$$M_{nm}=⟨{\psi }^{(n)}|\hat{z}|{\psi }^{(m)}⟩$$

$$⟨{\psi }^{(n)}|\hat{z}|{\psi }^{(m)}⟩$$

$$t_{ij}=⟨{\tilde{\psi }}^{(i)}|\mathcal{H}|{\tilde{\psi }}^{(j)}⟩=\sum _n{\epsilon }^{(n)}⟨{\tilde{\psi }}^{(i)}|{\psi }^{(n)}⟩⟨{\psi }^{(n)}|{\tilde{\psi }}^{(j)}⟩$$

$$\mathcal{H}=\epsilon \sum _n{c}_n^{†}{c}_n$$

$$c_k=\frac{1}{\sqrt{N}}\sum _n{e}^{-ikn}{c}_n$$

$$E_k^{\pm }=\epsilon +2t_3\cos (2k)\pm$$

$$g(\omega )=\frac{1}{N}\sum _{k,n}\delta (\omega -E_k^{(n)})$$

$${\mathcal{H}}_{\mathrm{dis}}=\sum _n{U}_n{c}_n^{†}{c}_n$$

$$P(U)=\frac{1}{\gamma \sqrt{2\pi }}\exp \left(-\frac{U^2}{2{\gamma }^2}\right).$$

$$L_i={\left[⟨{\Psi }^{(i)}|{\hat{n}}^2|{\Psi }^{(i)}⟩-⟨{\Psi }^{(i)}|\hat{n}|{\Psi }^{(i)}{⟩}^2\right]}^{1/2}$$

$$\hat{n}=\sum _{n}n{c}_n^{†}{c}_n.$$

$$t=\max \left(\min (|t_1|,|t_2|),|t_3|\right).$$

Full text

I Introduction

Soon after Watson and Crick’s discovery of the DNA double-helix structure WATSON1953 , Eley and Spivey ELEY1961 introduced the notion of efficient charge transport along the stacked $\pi$ orbitals of the bases. The mechanism of charge transport has been the subject of numerous studies in the intervening years, with renewed interest fuelled recently by both biological and technological considerations. Over a decade ago, Barton and co-workers observed distance-independent charge transfer between DNA-intercalated transition-metal complexes MURPHY1993 and argued that it would be relevant for biology and biotechnology. More recent electron transport experiments on DNA have yielded widely varying results, showing alternatively insulating behavior BRAUN1998 ; DEPABLO2000 ; STORM2001 ; ZHANG2002 ; BABA2006 , semiconducting behavior PORATH2000 , Ohmic conductivity FINK1999 ; CAI2000 ; TRAN2000 ; COHEN2005 , and proximity induced superconductivity KASUMOV2001 . The large number of relevant variables endemic to such experiments, like the DNA-electrode contact, and the rich variety of structures that DNA can assume, are the causes of variability in the experimental measurements (for a recent review of transport theory and experiments see Ref. ENDRES2004 ).

Specifically, there is a large diversity of the DNA forms in terms of its composition, length, and structure. Experiments done long ago, suggested that DNA substantially longer than its natural length (also referred to as “overstretched DNA”) can undergo a transition to an elongated structure up to twice the length of relaxed DNA WILKINS51 . This was also confirmed by recent single molecule stretching experiments SMITH1996 ; CLUZEL1996 ; STRICK2000 , which showed that the molecule can be reversibly stretched up to 90% of its natural length. Such important deformations of the double helix may occur in biological environments. Stretching of DNA is also related to cellular processes, such as transcription and replication. For example, proteins often induce important local distortions in the double helix while they diffuse along the molecule in search of their target sequences. The electronic and transport properties of DNA are directly influenced by its different conformations as well as by environmental factors, such as counterions, impurities or temperature. A full account of these effects based on a realistic, atomic scale description of the structure and the electronic properties challenges the capabilities of theoretical models.

Theoretical efforts to understand the electronic behavior and transport in DNA can be divided into two general categories:

(i) Model calculations that use effective hamiltonians and master equations to describe the dynamics of electrons and holes in DNA (see, for instance, Refs. YAMADA2004 ; IGUCHI2004 ; Yi2003 ; Chang2004 ). Recent results CAETANO2005 have led to considerable insights concerning the sequence-independent delocalization of electronic states in DNA. The main limitation of such approaches lies in the difficulty of determining accurate values for the parameters in the effective hamiltonians.

(ii) Ab initio calculations that can provide an accurate and detailed description of the electronic features RNBARNETT2001 ; ARTACHO2003 ; ALEXANDRE2003 . These approaches are typically limited to a small number of atoms due to computational costs, and cannot readily handle the full complexity of DNA molecules in various conformations. In particular, stretching of DNA can induce a very significant deviation from the B form which is stable under normal conditions in aqueous solution. Such structural distortions are bound to have a profound effect on the electronic behavior. A realistic description of these effects makes it necessary to handle both the atomic scale features and the overall state of the macromolecule.

In the present work, we address the problem of DNA stretching effects on the electronic states and the electron localization by providing a bridge between the two extremes of the length scale; a similar methodology was recently used to study hole transfer in DNA SENTHI05 . Theoretically, there are different ways of pulling the opposite ends of the DNA strands, leading to different stretched DNA forms, which are determined largely by base pair reorientations. Here, we use the poly(CG)-poly(CG) structures obtained in the pioneering study of Lebrun and Lavery LEBRUN1996 as the representative structure for stretching effects. This study modeled the adiabatic elongation of selected DNA molecules in two modes of stretching, corresponding to pulling on opposite 3’-3’ ends or 5’-5’ ends of the molecule: In the 3’-3’ stretching mode, the DNA helix is unwound leading to a ribbon-like structure, while in the 5’-5’ stretching mode the DNA helix contracts.

We begin with a set of detailed calculations for the electronic structure of DNA bases (A,T,C,G) and representative base pairs (AT-AT, CG-CG, AT-CG, CG-GC) in various relative configurations, as they are likely to appear in the stretched forms, These calculations are based on density-functional theory DFT ; LDA and serve to set the stage for more extensive calculations which employ successive levels of approximations necessary to handle the computational demands. Specifically, we extract the salient features of electronic structure of the individual DNA bases and base pairs from the ab initio calculations; these are compared to an efficient and realistic semi-empirical model ELSTNER1998a , in order to establish the validity of the latter approach. At this intermediate scale, we consider an 18 base pair poly(CG)-poly(CG) DNA sequence which has been stretched by 30%, 60% and 90% relative to the natural length of the unstretched B form. The atomic structure of these forms has been established by Lebrun and Lavery LEBRUN1996 , using empirical interatomic potentials. We next use the information from this approximate description to build an effective hamiltonian for the electronic behavior at much larger scales. This allows us to describe electron localization, due to the combined effects of stretching and environmental factors, over mesoscopic to macroscopic length scales. The essence of the approach and the different scales involved are shown schematically in Fig. 1. We emphasize that we address here issues related only to dry and neutral DNA structures, where the negatively charged groups on the backbone are passivated by protons, conditions that are relevant to the experiments we consider for comparison to our theoretical results; water molecules or counterions (such as Na+) are not considered in our calculations.

II Theoretical Methods

II.1 Ab initio calculations

As our first step toward establishing the electronic behavior of dry, neutral DNA, we study the nature of electronic states in individual bases and in base pairs. For these calculations we used three different implementations of density-functional theory DFT : a method that uses atomic-like orbitals as the basis SIESTA , one that uses plane waves VASP and a third that uses a real-space grid HARES . In all three approaches, we used the same exchange-correlation functional in the local-density approximation LDA , for consistency and simplicity. More elaborate approximations to exchange-correlation effects, such as the generalized gradient approximation GGA , do not provide any improvement in describing the physics of these weakly interacting units. In each method we used pseudopotentials to represent the atomic cores, of the Trouiller-Martins type Trouiller in SIESTA, the Vanderbilt ultrasoft type Vanderbilt in VASP and the Hammann-Schluter-Chiang type HSC in HARES, with computational parameters (number of orbitals in basis, plane-wave kinetic energy cutoff and grid spacing) that ensure a high level of convergence. These calculations provide a thorough check on the consistency of various computational schemes to reproduce the electronic features of interest. The results are in excellent agreement across the three approaches. Since in these calculations there are no adjustable parameters, we refer to them in the following as ab initio results.

II.2 Construction of semi-empirical model

The stretched forms contain a large number of atoms, typically beyond what can be efficiently treated with the ab initio methods used for the DNA bases and base pairs. Accordingly, for the electronic structure calculations of these structures we use an efficient semi-empirical quantum-mechanical approach which employs a minimal basis set ELSTNER1998a . The consistency of this approach is then verified against the ab-initio calculations. Within the semi-empirical scheme, the electronic eigenfunctions are expressed as

[TABLE]

where the basis set $|{\varphi_{n}}\rangle$ includes the $s$ and $p$ atomic orbitals for each atom in the system. The coefficients $c^{(n)}_{\nu}$ are numerical constants, with $|c^{(n)}_{\nu}|^{2}$ giving the weight of orbital $|{\varphi_{n}}\rangle$ to the electronic wavefunction. This method uses a second order expansion in the electronic density to obtain the total energy and takes into account self-consistently charge transfer effects which are important for biological systems. The method gives results for the band gaps that are in excellent agreement with those of the ab initio approaches described above (see Refs. MARAGAKIS2002 ; DEPABLO2000 ).

The highest occupied and lowest unoccupied molecular orbitals (HOMO and LUMO, respectively, also referred to collectively as “frontier states” in the following) are extended over the entire structure in Bloch-like wave functions. In order to describe electron hopping and localization, we need to express these in terms of a basis of Wannier-like states that are localized on the individual bases. To this end, we construct maximally localized states on single base pairs by taking linear combinations of the HOMO and LUMO states from the wavefunctions of Eq. (1). The maximally localized states will then be used to calculate the hopping parameters in the effective 1D hamiltonian. Using the extended electronic states $|{\psi^{(n)}}\rangle$ of the frontier states, with corresponding energies $\varepsilon_{(n)}$, we define the maximally localized states $|{\tilde{\psi}^{(i)}}\rangle$ through the unitary transformation

[TABLE]

which minimizes the sum of the variances

[TABLE]

under the constraint $\langle{\tilde{\psi}^{(i)}}|\tilde{\psi}^{(j)}\rangle=\delta_{ij}$ where $z$ is the position along the helical axis. Similar and more general methodologies have been developed in the past for obtaining maximally localized states from extended ones MARZARI97 ; SGIAVORELLO01 . Due to the invariance of the trace, the first term in Eq. (3) is independent of the unitary transformation and the problem is simplified to one of maximizing the second term on the right-hand side with the same orthonormality constraint. Carrying out the minimization, we arrive at the equation

[TABLE]

where

[TABLE]

By inspection, we see that $\zeta$ is maximized when $z_{n}=z_{m}$ for all $m$ and $n$, corresponding to maximally delocalized states. On the other hand, $\zeta$ is minimized when the states $|{\tilde{\psi}^{(n)}}\rangle$ are the eigenfunctions of the position operator $\hat{z}$ within the HOMO or LUMO subspace. Therefore, the problem is further reduced to constructing and diagonalizing the matrix

[TABLE]

which has the eigenvectors $\langle{\psi^{(n)}}|{\tilde{\psi}^{(i)}}\rangle$ that provide the desired transformation given in Eq. (2). The eigenvalues $z_{n}$ are the positions of the localized states. To evaluate the matrix elements we use the approximation

[TABLE]

where $S_{\mu\nu}=\langle{\varphi_{\mu}}|\varphi_{\nu}\rangle$ is the overlap matrix between the two atomic orbitals and $z_{\mu\nu}=\frac{z_{\mu}+z_{\nu}}{2}$ is the average $z$-value for the atoms located at sites given by the labels $\mu$ and $\nu$. Once the localized states are constructed, the hopping parameters can be computed as

[TABLE]

recalling that the quantities $\langle{\psi^{(n)}}|{\tilde{\psi}^{(i)}}\rangle$ are determined from the transformation described above.

Having defined the maximally localized states in terms of the electronic wavefunctions from the all-atom calculations, we next produce an effective tight-binding hamiltonian, which allows us to study electron hopping along the DNA double helix. This approach has also been used in a recent study on functionalized carbon nanotubes LEE05 . In our effective hamiltonian, we consider hopping between first and second neighbors along the helix, and denote the hopping matrix elements according to the scheme shown in Fig. 2 for the HOMO state of the poly(CG)-poly(CG) structure (all other frontier states involve exactly the same type of hopping matrix elements):

[TABLE]

where $n$ represents the $n^{th}$ base pair along the helical axis and we have neglected spin indices because they are unimportant for our analysis. Note that there is a difference between hopping elements connecting even and odd sites to their neighbors ($t_{1}$ and $t_{2}$ terms in the effective hamiltonian of Eq. (9)), due to the asymmetry in the structure illustrated in Fig. 2. Performing a Fourier transform on the electron creation and annihilation operators

[TABLE]

gives a hamiltonian which has coupling between momenta $k$ and $k+\pi/a$. By doubling the unit cell (and reducing the Brillouin Zone by a factor of two), this can finally be diagonalized to obtain the eigenvalues

[TABLE]

with the momentum sum carried out over the reduced Brillouin Zone. With these expressions for the band structure energies, the density of states (DOS)

[TABLE]

can be readily obtained. These quantities are essential in describing electron localization along the DNA double helix under different conditions.

II.3 Disorder and Localization length

In order to quantify the amount of localization that is expected in stretched DNA forms, we add a term to the hamiltonian in Eq. (9) of the form

[TABLE]

which is meant to emulate disorder arising from a variety of sources such as interaction of the DNA bases with stray water molecules and ions, or interaction with the substrate. $U_{n}$ are uncorrelated random energy variations chosen according to a Gaussian distribution of zero mean and width $\gamma$

[TABLE]

Once the disorder hamiltonian is constructed with a specific set of random on-site energies, by direct diagonalization we find the eigenstates $|{\Psi^{(i)}}\rangle$ of ${\cal H}+{\cal H}_{\rm dis}$ (we use capital symbols to denote the new wavefunctions from the hamiltonian that includes the disorder term) and then calculate the localization length defined as

[TABLE]

where

[TABLE]

For a single-hopping model with weak disorder, the localization length scales as $L\sim\left(t/\gamma\right)^{2}$ for electrons near the middle of the band THOULESS1972 , with $t$ the hopping matrix element which determines the band width. The more complicated effective hamiltonian considered here is not amenable to simple analytic treatment.

III Results and Discussion

We begin our discussion with an overview of electronic states in single bases and isolated base pairs. The structure of the base pairs is shown in Fig. 3 with the atoms in each base labeled for future reference. These calculations will set the stage for a proper interpretation of the behavior in the stretched and unstretched dry, neutral DNA helix.

III.1 Frontier states

The frontier states in the base pairs are related to only one component of the pair for both AT and CG. This is shown in Fig. 4. Specifically, the HOMO state of the AT pair is exactly the same as that of the HOMO state of the isolated A, and the LUMO state of AT the same as that of the isolated T. Similarly, the HOMO state of CG is identified with that of the isolated G and the LUMO state with that of the isolated C. Thus, the purines (A or G) give rise to the HOMO state, while the pyrimidines (T or C) are responsible for the LUMO states of each pair. It is clear from the same figure, that essentially all atomic $p_{z}$ orbitals which belong to a purine or pyrimidine contribute to the respective HOMO or LUMO $\pi$ state of the base pair. This is in agreement with calculations on the optical absorption spectra of DNA bases and base pairs VARSANO2006 . A closer inspection of Fig. 4 shows that the molecular frontier states of both AT and CG can be identified as similar contributions (up to sign changes) from specific groups of carbon and nitrogen atoms. Specifically, in the purines (A and G) three distinct groups of atoms are mainly involved in forming the HOMO orbital and include atoms (C8-N7), (C2-N3) and (N1-C6-C5-C4-N9), respectively. In the pyrimidines (T and C) the main groups involved in forming the LUMO orbital are two, (C4-C5-N1) and (N3-N7-C6). In both base pairs the atoms that are less involved in the frontier molecular states are the carbon atoms that form a double bond with an oxygen atom, such as C2 of A and C and the four-fold bonded C7 atom of A.

The frontier states are very little affected when the two components of the base-pair are separated along the direction in which they are hydrogen-bonded. To demonstrate this, we show in Fig. 5 the change in the eigenvalues of the frontier states in AT and CG as a function of the distance between the two atoms that are bonded to the two backbones (we call this the backbone distance). For both base pairs the nitrogen atoms labeled N1 and N9, are the ones attached to the backbone (see Fig. 3). In order to obtain realistic structures, for each value of the backbone distance we hold the atoms of each base that are bonded to the backbone fixed and allow all other atoms to relax fully. These calculations were performed with the SIESTA code SIESTA and the relaxed configurations were used as input to calculate the electronic structure with the other two methodologies VASP ; HARES . In Fig. 5 we show complete results from the SIESTA calculations and selected results from one of the other two approaches.

The results of Fig. 5 show clearly that only in the region where the backbone distance becomes significantly smaller than the equilibrium value, interaction between the two bases shifts the eigenvalues of the electronic states appreciably, but even then the shifts are relatively small for the frontier states. It is also noteworthy that the band gap of the AT pair is significantly larger ($\sim 3$ eV) than that of the CG pair ($\sim 2$ eV) and that the frontier states of CG lie within the band gap of the AT pair. This observation is important because it indicates that in an arbitrary sequence of base pairs, the frontier states will be associated with those of the CG pairs. This statement is verified by calculations of electronic states in the AT-AT, CG-CG and AT-CG base pair combinations, to which we turn next.

For more detailed comparisons, we collect in Table 1 the eigenvalues of the frontier states for the DNA pairs and the pair combinations, at different equilibrium configurations in the three relevant variables, the backbone distance, the axial distance and the rotation angle. Some results on the CG-GC base pair combination are also shown, to allow for comparison to the poly(C)-poly(G) sequence.

When two base pairs are stacked on top of each other, there are two degrees of freedom for motion of one relative to the other: a separation along the helical axis, which we will call axial distance, and a relative rotation around the helical axis. We take the helical axis to be that which corresponds to stacking of successive base pairs in the B form of the DNA double helix. According to the notation of Fig. 3, the helical axis for both base-pairs is normal to the line connecting atoms C4 and C6 and is closer (about one third of their distance) to the purine atom C6. For each configuration we fix the atoms that are bonded to the backbone at a given relative position and allow all other atoms to relax, as was done in the calculations involving the backbone distance discussed above. In Fig. 6 we show the behavior of electronic eigenvalues as a function of the axial distance and the rotation angle. As above, the eigenvalues show little dependence on these two variables, except for rather small values of the axial distance which correspond to unphysically small separation between the two base pairs.

What is also remarkable in the above results, is that in the AT-CG combination, the frontier states are clearly identified with those corresponding to the CG pair exclusively, which has the smaller band gap (see Fig. 6). Moreover, we note that the band gap of the poly(C)-poly(G) sequence, as calculated by the semi-empirical method based on a minimal atomic orbital basis ELSTNER1998a is in excellent agreement with the value obtained from the SIESTA calculation (2.0 eV and 2.1 eV, respectively). The band gap is expected to be significantly smaller in the case of wet DNA and in the presence of counterions, as shown in Ref. GERVASIO2002 , for a Z-DNA helix. The band gaps between all three ab initio methods are identical within the accuracy of these methods. The nature of electronic wavefunctions obtained by the different methods is also in good qualitative agreement. Accordingly, in the rest of this paper we focus our attention to electron localization in the dry, neutral poly(CG)-poly(CG) sequence, and employ the results of the semi-empirical electronic structure method.

III.2 Hopping electrons

In Fig. 7, we show the unstretched and the three stretched forms of the poly(CG)-poly(CG) sequences at 30%, 60%, 90% elongation, along with the features of the frontier states. For visualization purposes, we represent the calculated wavefunction magnitude of the frontier states by blue (HOMO) and red (LUMO) spheres, centered at the sites where the atomic orbitals are located. The radius of the sphere centered on a particular atom is proportional to the magnitude of the dominant coefficient $|c^{(n)}_{\nu}|^{2}$ at this site (see Eq. (1)), which is essentially proportional to the local electronic density. It is evident from this figure that the nature of the orbitals themselves, represented by the radii of the colored spheres, does not change much in the different stretched DNA forms, but the overlap between orbitals at neighboring bases is affected greatly by the amount of stretching. For the poly(CG)-poly(CG) sequence shown, the HOMO orbitals are always associated with the G sites for all the stretching modes, while the LUMO orbitals are related to the C sites. However, as the DNA becomes more elongated, the orbitals overlap even less and become localized for high stretching modes. The elongation to the overstretched form is achieved by changing the dihedral angle configuration of the DNA backbone, which leaves the local part of the orbitals essentially intact. Note how the orbitals rotate and spread out as the structure is being ovestretched, following the rotation of bases.

We now turn to a discussion of the results for the hopping matrix elements of Eq. (9). Our discussion here is relevant to what happens when the occupation of a frontier state is changed from complete filling (for the HOMO) or complete depletion (for the LUMO), that is, the physics of small amounts of hole or electron doping. In Table 2 we give the values for $\varepsilon,t_{1},t_{2},t_{3}$ (see Fig. 2) for the two frontier states of the unstretched poly(CG)-poly(CG) DNA form. The hopping matrix elements for the HOMO state involve only the G sites; those for the LUMO state involve only the C sites. As a consistency check, we have also calculated matrix elements for farther neighbors and found those to be much smaller in magnitude. We have calculated the values of $t_{1},t_{2},t_{3}$ by repeating the same procedure as above for the stretched forms of the poly(CG)-poly(CG) DNA sequence. We note that if $t_{2}=t_{3}=0$ electrons will not be able to migrate along the DNA molecule even if $t_{1}$ is quite large, because at least one of the other two hops is necessary for migration (see Fig. 2). From this simple picture, it is evident that the conductivity will be determined by which matrix element dominates. Quantitatively, the “bottleneck” hopping matrix element is given by

[TABLE]

In Fig. 8 we show the value of the “bottleneck” hopping matrix element calculated as a function of stretching. This indicates that hopping conductivity will dramatically decrease by several orders of magnitude upon stretching the molecule and that the hopping will decrease more from stretching in the 3’-3’ mode than in the 5’-5’ mode. This is due to the conformational changes induced by the different stretching modes, described earlier.

III.3 Localization length

The significant dropping of the hopping matrix elements upon stretching as described in the previous section is indicative of electron localization with a weak amount of disorder. To investigate this possibility in detail, we focus on effects of stretching in the 3’-3’ mode. The evolution of the density of HOMO states upon stretching is shown in Fig. 9; similar behavior is observed for the LUMO states. The dramatic narrowing of the DOS width (equivalent to reduced dispersion in a band-structure picture) is strongly suggestive of electron localization ANDERSON1958 , in this case induced by stretching. This localization length is controlled by the hopping elements $t$, since $\varepsilon$ is the same at each site.

For a more quantitative description, we show in Fig. 9 the localization length $L_{i}$ for each eigenstate for a $1500$ base-pair DNA strand under different amounts of stretching. The value of $L^{(i)}$ for each state is obtained from Eq.(15), with disorder strength $\gamma=0.3$ meV, which determines the width of the gaussian given in Eq. (14). This disorder strength is much smaller than the band width of the unstretched DNA, but becomes comparable to the band width as the molecule is stretched. The magnitude of such variations in on-site energies is consistent with those produced by the dipole potential terms, for instance, due to the presence of a stray water molecule situated on the substrate roughly 15 Å away from the DNA bases. We find that changing the value of $\gamma$ by an order of magnitude (either smaller or larger) does not affect the qualitative picture presented here. Note that the localization length is not a strict function of the energy, as it depends on the disorder near where a given state happens to be localized. As the molecule is stretched, the localization length dramatically decreases until, for 60% stretching, the eigenstates are completely localized on single base pairs.

The charge localization length as a function of DNA stretching has been recently studied in the experiment of Heim et al. Heim2004 . This study focuses on $\lambda$-DNA which has an irregular sequence of base pairs, and can be compared to our theoretical results for poly(CG)-poly(CG) recalling that the frontier states even for a random sequence are associated with those of the CG base-pairs. In the experiment, ropes of $\lambda$-DNA on a substrate are overstretched by a receding meniscus technique. The DNA ropes in this experimental setup are slightly positively charged, corresponding to a depletion of a few electrons per 1000 base pairs. We suggest that this situation is approximated by the structures of dry and neutral DNA that we considered above. Electrons were injected into the DNA and the resulting localization length was measured by an electron force microscope. For the unstretched DNA, the charge was found to delocalize across the entire molecule, extending over a length of several microns. On the other hand, the charge injected into the overstretched DNA is localized, extending over a few hundred nanometers only. This is qualitatively consistent with the picture that emerges from our theoretical analysis, and is even in reasonable quantitative agreement: the degree of localization in experiment, measured by the ratio of length scales going from unstretched to stretched DNA structures, is approximately two orders of magnitude, while the same quantity in our calculations, going from unstretched to 60% stretched DNA is $\sim 10^{3}$.

IV Summary

We have described and implemented a multiscale method to derive effective hamiltonian models that are able to capture the dynamics of conduction and valence electrons in stretched DNA, starting from ab initio, all-atom quantum mechanical calculations. The ab initio simulations revealed that the frontier states in the base pairs are related to only one component of the pair. The purines were found to be associated with the HOMO states while the pyrimidines with the LUMO states. In the AT-CG combination the frontier states are identified with those of the CG pair. For all combinations of bases and base pairs studied here, the nature of these states was not affected by separation of the bases or base pairs along different directions or rotation along the helical axis.

Turning to the next length scale and the semi-empirical calculations, we have calculated the “bottleneck” matrix elements for electron hopping along the DNA molecule, as a function of stretching. These show a significant decrease with elongation of DNA, which is stronger for stretching in the 3’-3’ mode than in the 5’-5’ mode. We were able to show quantitatively that stretching of DNA dramatically narrows the DOS width of frontier states. A small amount of disorder produced by environmental factors will naturally lead to localization of the electrons along the DNA. Our estimate for the degree of localization, based on a reasonable (and quite small) amount of disorder in the on-site energies for the electron states, is in very good agreement with recent experimental observations. This provides direct validation for the consistency and completeness of the multiscale method presented here.

Acknowledgements: The authors are grateful to Richard Lavery for providing the overstretched structures. MF acknowledges support by Harvard’s Nanoscale Science and Engineering Center, funded by the National Science Foundation, Award Number PHY-0117795.