Base-Pairs’ Correlated Oscillation Effects on the Charge Transfer in Double-Helix B-DNA Molecules

By introducing a suitable renormalization process, the charge carrier and phonon dynamics of a double-stranded helical DNA molecule are expressed in terms of an effective Hamiltonian describing a linear chain, where the renormalized transfer integrals explicitly depend on the relative orientations of the Watson–Crick base pairs, and the renormalized on-site energies are related to the electronic parameters of consecutive base pairs along the helix axis, as well as to the low-frequency phonons’ dispersion relation. The existence of synchronized collective oscillations enhancing the π-π orbital overlapping among different base pairs is disclosed from the study of the obtained analytical dynamical equations. The role of these phonon-correlated, long-range oscillation effects on the charge transfer properties of double-stranded DNA homopolymers is discussed in terms of the resulting band structure.


Introduction
In physiological conditions, the DNA double helix exhibits a full-fledged three-dimensional (3D) geometry, where every two consecutive Watson-Crick base pairs (bps) stand nearly parallel to each other, and they are twisted by a certain angle (θ 0 36 • in equilibrium conditions) around the helix axis. In their pioneering work, Eley and Spivey pointed out that a double-stranded DNA (dsDNA) molecule might behave as a one-dimensional (1D) aromatic crystal displaying a π-π based electrical conductivity along the helical axis [1]. The reasoning behind this proposal was that dsDNA's nucleobases adenine (A), guanine (G), cytosine (C), and thymine (T) are aromatic compounds whose atomic p z orbitals perpendicular to the plane of the base can form rather delocalized π bonding and π * antibonding molecular orbitals. If the orbital overlap between the bps is strong enough, this could lead to extended electronic states along the helical axis, thereby promoting charge transfer (CT) between consecutive bps in an efficient way over long distances through the aromatic bp stack within the DNA helix. Accordingly, CT depends on the intimate coupling among stacked bases, as determined by their relative separation and twist angle, and hence, any perturbation in that stacking, altering the optimal overlapping face-to-face configuration (θ 0 = 0), will significantly affect DNA charge migration. Consequently, one expects structural fluctuations to be an important factor, influencing charge carriers' transport through dsDNA molecules [2][3][4]. In fact, at physiological temperatures, the relative orientation of neighboring bases becomes a function of time, thereby modifying their mutual overlapping in an oscillatory way. The motion of bases can either occur in a synchronized manner (normal modes' propagation) or incoherently. The role of thermal fluctuations on the CT efficiency has been discussed in a number of previous works, where the structural fluctuations of the DNA double helix are described by sampling the initial angular velocities and twist angles from a Boltzmann distribution at a given temperature [5][6][7][8][9][10][11][12][13]. Not surprisingly, it was found that the uncorrelated motion of bps, randomly twisting back and forth around the helix axis due to thermal fluctuations, generally reduces π-π stacking overlap, hence degrading the CT efficiency.
On the contrary, the presence of synchronized, collective twist motions of the Watson-Crick bps in DNA duplexes can efficiently enhance the π − π orbital overlapping between non-consecutive bps via a long-range, phonon-correlated tunneling effect [14,15]. In this work, we will focus on coherent charge transport promoted by the coupling between both twist and radial vibration modes and charge motion through duplex DNA, thereby extending previously obtained results [16]. In order to analyze the interplay between the dsDNA low frequency bps' dynamics and CT efficiency, we will study the coupling between the oscillations of complementary bases along the transversal direction and the twisting motion of each bp as a whole through the helical sugar-phosphate backbone structure, explicitly taking into account its characteristic helical geometry, which has been shown to be very important in biological processes, such as denaturation and transcription [17]. To this end, we will explicitly take into account the stacking interaction, mediated by the orbital overlapping between adjacent bps along the helix, as well as hydrogen bond stretch motions, as described in the Peyrard-Dauxois-Bishop (PDB) model [18][19][20], in the phonon dynamical equations. In our approach, the charge carrier dynamics through a helical dsDNA molecule is expressed in terms of an effective renormalized Hamiltonian describing a diatomic linear chain, where the renormalized transfer integrals explicitly depend on the relative orientations of the Watson-Crick bps and the renormalized on-site energies are related to both the electronic parameters of consecutive codon units along the helix axis, as well as the low-frequency phonon dispersion relation. The corresponding effective hopping terms include both helical and dynamical effects in an intertwined fashion, allowing for a unified treatment of charge-lattice coupled dynamics in a fully analytical way. Thus, we disclose a number of remarkable symmetries of the motion equations themselves, which may be implemented with accurate charge transfer parameters derived from quantum chemistry and molecular dynamics approaches in a straightforward way. Our main conclusion is that a significant improvement of CT can occur in dsDNA via charge-phonon coupling mediated by synchronized helical waves stemming from collective, long-range correlated bps' oscillation modes.
The paper is organized as follows. In Section 2, we introduce the model Hamiltonian describing the lattice and electronic dynamics. The lattice contribution is expressed in cylindrical coordinates, in order to explicitly take into account the 3D geometry of the double helix DNA molecule. Then, we derive the π − π electronic coupling term describing the molecular orbitals' overlap through the helix axis in terms of these cylindrical coordinates. This term, describing the transfer integral between successive bps, allows one to relate the dynamical behavior to the CT process by means of a suitable tight-binding electronic model. A convenient feature of our adopted approach is that the resulting 3D fishbone model can be properly mapped into a mathematically simpler effective 1D chain model, still retaining much physico-chemical information in the corresponding renormalization parameters. In Section 3, we obtain the linearized canonical equations of motion for twist and radial lattice variables. In doing so, we introduce a number of characteristic frequencies along with their related time scales. For the sake of simplicity, in Section 3.2, we focus on the dynamics of homopolymer dsDNA molecules, showing the presence of collective oscillations in the form of helical waves. The related dispersion relations for the acoustic and optical branches are analytically derived, and the obtained results are compared to some available experimental results. In this section, we also disclose a very interesting relationship between the dynamics of the dsDNA molecule as a whole and that corresponding to its codon building blocks. In Section 4, we solve the Schrödinger equation for the effective 1D Hamiltonian previously introduced, making explicit use of the helical wave solutions in the transfer integral term. In this way, the charge-phonon coupling effect is fully incorporated in the resulting 3D energy spectrum. Finally, the main conclusions of this work are summarized in Section 5.

DNA Model Hamiltonian
Two kinds of order coexist in biological DNA, each one related to two separate subsystems in the DNA helix, namely the nucleobase and backbone systems [21]. The informative chemical order determined by the sequence of Watson-Crick bps can be suitably characterized by ab-initio quantum chemistry calculations [22,23], which properly highlight the emergence of molecular orbitals, as is shown in Figure 1b. In order to describe most basic properties of dsDNA molecules, we must consider a model Hamiltonian accounting for different scales of time and space by means of an adequate choice of generalized coordinates including both electronic and dynamic degrees of freedom. According to the Born-Oppenheimer approximation, the lattice and charge dynamics of a dsDNA molecule can be split in terms of the general Hamiltonian H = H l + H e , where H l describes the double strand dynamics, and H e describes the CT across nucleobases, as illustrated in the structural lattice model and the electronic tight-binding model depicted in Figure 1a,c, respectively. Version November 9, 2020 submitted to Materials 3 of 18

82
Two kinds of order coexist in biological DNA, each one related to two separate subsystems in 83 the DNA helix, namely the nucleobase and backbone systems [21]. The informative chemical order 84 determined by the sequence of Watson-Crick bps can be suitably characterized by ab-initio quantum 85 chemistry calculations [22,23], which properly highlight the emergence of molecular orbitals, as it is 86 shown in Figure 1b. In order to describe most basic properties of dsDNA molecules, we must consider  Diagram of the dsDNA electronic model showing the π-π channel transfer integrals t XY n,n±1 (cylindrical rods), the interstand transfer integrals between complementary bases t XY ⊥ (perpendicular bars), the glucosidic transfer integrals t P between nucleobases (spheres) and sugar-phosphate groups (cubes), and the on-site energies of these groups γ j .

94
In our lattice model we treat each nucleotide (base + sugar + phosphate) as a point mass, helically arranged and mutually connected by means of elastic rods describing: (1) the sugar-phosphate backbone along a given strand, and (2) the interstrand H-bonding between complementary bases (see Figure 1a) [24][25][26]. We explicitly take into consideration the mass difference among the four nucleobases, namely, m G = 347.05, m C = 307.05, m A = 331.06, m T = 322.05 amu, and so we realize that the mass of each bp as a whole is essentially the same, i. e., M ≡ m G + m C ∼ = m A + m T ∼ = 653.5 ± 0.5 amu. Adopting the reference frame indicated in Figure 1a the position of the nth nucleobase can be expressed as x n = r n cos ϕ n , y n = r n sin ϕ n , and z n = cϕ n , where n labels the considered bp along the dsDNA, r n and ϕ n are usual cylindrical coordinates, and c = h 0 /θ 0 , where h 0 0.34 nm is the equilibrium distance between two successive bp planes along the Z axis in the B-DNA form, and θ 0 is the equilibrium relative angular separation between neighboring bps. We note that in this model Full-atom dsDNA model with surfaces of constant charge density for the states corresponding to the π-like lowest unoccupied molecular orbital (LUMO, in red) of the C bases and the highest occupied molecular orbital (HOMO, in blue) of the G bases of a polyG-polyC molecule in the dry conditions A form (Courtesy of Emilio Artacho) [23]. (c) Diagram of the dsDNA electronic model showing the π-π channel transfer integrals t XY n,n±1 (cylindrical rods), the interstrand transfer integrals between complementary bases t XY ⊥ (perpendicular bars), the glucosidic transfer integrals t P between nucleobases (spheres) and sugar-phosphate groups (cubes), and the on-site energies of these groups γ j .

Lattice Hamiltonian
In our lattice model, we treat each nucleotide (base + sugar + phosphate) as a point mass, helically arranged and mutually connected by means of elastic rods describing: (1) the sugar-phosphate backbone along a given strand and (2) the interstrand H-bonding between complementary bases (see Figure 1a) [24][25][26]. We explicitly take into consideration the mass difference among the four nucleobases, namely, m G = 347.05, m C = 307.05, m A = 331.06, and m T = 322.05 amu, and so, we realize that the mass of each bp as a whole is essentially the same, i.e., M ≡ m G + m C ∼ = m A + m T ∼ = 653.5 ± 0.5 amu. Adopting the reference frame indicated in Figure 1a, the position of the nth nucleobase can be expressed as x n = r n cos ϕ n , y n = r n sin ϕ n , and z n = cϕ n , where n labels the considered bp along the dsDNA, r n and ϕ n are the usual cylindrical coordinates, and c = h 0 /θ 0 , where h 0 0.34 nm is the equilibrium distance between two successive bp planes along the Z axis in the B-DNA form, while θ 0 is the equilibrium relative angular separation between neighboring bps. We note that in this Materials 2020, 13, 5119 4 of 19 model, the distance between successive bps along the Z direction is proportional to the twist angle. In this way, the helical structure is naturally preserved during the dynamical evolution, in agreement with dispersion relation data reported from inelastic X-ray scattering measurements [27]. Thus, we can express the Euclidean distance between adjacent bases along a given strand as: where we defined θ n,n±1 = ± (ϕ n±1 − ϕ n ) as the relative angle between two neighboring bps, and ρ n = r n − R 0 is the radial displacement about the equilibrium position (R 0 = 1 nm). We will further assume that the location of the dsDNA as a whole remains fixed, so that the center of mass is constant for each bp. Therefore, the radial displacements about the equilibrium position satisfy the relationshipρ n = λ n ρ n , where λ n ≡ m n /m n , and henceforth, the upper bar denotes the physical magnitudes of the complementary bases. Therefore,d n,n±1 for the opposite strand can be obtained by simply replacing ρ n withρ n in Equation (1). In equilibrium, both distances reduce to the value l 0 = d n,n±1 | eq. =d n,n±1 eq. = h 2 0 + 4R 2 0 sin 2 (θ 0 /2) 0.685 nm, where we adopted When describing the phonon dynamics of DNA at a molecular scale, one can disregard the inner degrees of freedom of the bases, since we can separate the fast vibrational motions of atoms about their equilibrium positions from the slower motions of molecular groups. Accordingly, we can write the dsDNA molecule lattice Hamiltonian as [16]: where n runs over the number N of bps, P ρ n and P ϕ n are the conjugate momenta of the nth bp radial and twist variables, respectively, and ξ = c 2 + R 2 0 1.147 nm is related to the helical geometry of the system, so that ξθ n,n±1 measures the helix arc length providing the shortest path between two points along a helical coil. In the limit of small radial and twist oscillations (r n R 0 , θ n,n+1 1), Equation (1) reads d n,n±1 = R 2 0 + c 2 θ n,n±1 ≡ ξθ n,n±1 , so that the Euclidean distance coincides with the helix arc length in this case [14].
The three elastic potential terms in Equation (2) describe the different interactions between the bases within the framework of the PDB model [18,19], namely: represents the radial stretching of the hydrogen bonds connecting complementary bases in the opposite strands of the double helix by means of Morse potentials of depth D n and width α n [19,28]. This potential term includes both the attraction due to the H-bonds forming the bps and the repulsion of the negatively charged phosphates in the backbone of the two strands, which is, in turn, screened by the surrounding solvent water molecules and positively charged counterions. The sequence dependence can be considered by adopting different values for the model parameters D n and α n , accounting for the different number of H-bonds in the G≡C and A=T bps [29]. The potential term: with u ± n,n+1 = (1 + λ n ) ρ n ± (1 + λ n+1 ) ρ n+1 , describes the stacking interaction between adjacent bps, whose role is to inhibit configurations with large relative radial displacements between neighboring pairs. This interaction is characterized by the exponential term that effectively modulates an otherwise harmonic radial oscillation. This term accounts for local constraints in nucleotide motions, measured in terms of the stacking stiffness k S n,n+1 and the interaction range b parameter, resulting in long-range cooperative elastic effects due to the distortion of hydrogen bonds and the overlap of the π-type orbitals [28]. The description of the radial degree of freedom via the non-linear potentials given by Equations (3) and (4) is more realistic than a purely harmonic approach and has been successful in capturing denaturation, as well as transcription initiation processes in several DNA model chains [19,[30][31][32]. Finally, the term: describes the harmonic interaction between neighboring bases along each backbone's strand.

The π-π Electronic Coupling
The radial and twist oscillations of bps have a significant impact on the molecular orbitals' overlap throughout the π-stacking, so that the resulting electronic transfer integrals' values explicitly depend on the dynamical degrees of freedom ρ n , ρ n±1 , and θ n,n±1 . The π and π * molecular orbitals are formed by the C, N, and O atomic p z orbitals perpendicular to the bps and pointing along the helical axis, as is illustrated in Figure 2. The p z orbitals from different bps couple by ppσ > 0 and ppπ < 0 hybridization, the different signs arising from the respective atomic orbital's lobe sign. According to the Slater-Koster theory, the transfer matrix element between two p z orbitals on neighboring bps is given by the combination of ppσ and ppπ hybridization contributions as V ij = V ppσ sin 2 ζ + V ppπ cos 2 ζ, where ζ measures the rise angle between successive bps (sin ζ = h ij /d ij ), and: where η ppπ and η ppσ describe the hybridization matrix elements, m is the electron mass, d ij = l 2 ij + h 2 ij is the Euclidean distance between atoms belonging to neighboring nucleobases, and R c describes the exponential tails of the atomic wave functions [33,34].
Version November 9, 2020 submitted to Materials 5 of 18 harmonic radial oscillation. This term accounts for local constraints in nucleotide motions, measured in terms of the stacking stiffness k S n,n+1 and the interaction range b parameter, resulting in long-range cooperative elastic effects due to the distortion of hydrogen bonds and the overlap of the π-type orbitals [28]. The description of the radial degree of freedom via the non-linear potentials given by Eqs. (3) and (4) is more realistic than a purely harmonic approach and has been successful in capturing denaturation as well as transcription initiation processes in several DNA model chains [19,[30][31][32]. Finally, the term describes the harmonic interaction between neighboring bases along each backbone's strand. The radial and twist oscillations of bps have a significant impact on the molecular orbitals overlap throughout the π-stacking, so that the resulting electronic transfer integrals values explicitly depend on the dynamical degrees of freedom ρ n , ρ n±1 and θ n,n±1 . The π and π * molecular orbitals are formed by the C, N and O atomic p z orbitals perpendicular to the bps and pointing along the helical axis, as it is illustrated in Figure 2. The p z orbitals from different bps couple by ppσ > 0 and ppπ < 0 hybridization, the different signs arising from the respective atomic orbital's lobes sign. According to the Slater-Koster theory the transfer matrix element between two p z orbitals on neighboring bps is given by the combination of ppσ and ppπ hybridization contributions as V ij = V ppσ sin 2 ζ + V ppπ cos 2 ζ, where ζ measures the rise angle between successive bps (sin ζ = h ij /d ij ), and where η ppπ and η ppσ describe the hybridization matrix elements, m is the electron mass, is the Euclidean distance between atoms belonging to neighboring nucleobases, and R c describes the 113 exponential tails of the atomic wave functions. [33,34]   The π-π coupling generally involves one base (say, Y) of the (n + 1)th bp and another base (say, X) of the nth bp, both belonging to the same strand (i.e., 5'-XY-3'), so that the π-π transfer integral between successive bps along the dsDNA helix is then given by: where N 1 and N 2 are the number of p z orbitals in the bps n + 1 and n, respectively, and c n j is the jth linear combination of atomic orbitals (LCAO) coefficient of the π molecular frontier orbital (HOMO or LUMO) of bp n. Making use of Equation (6) in Equation (7), assuming that the mean distances among atoms belonging to different nucleobases can be roughly approximated as d ij d and l 2 ij l 2 n,n±1 = r 2 n + r 2 n±1 − 2r n r n±1 cos θ n,n±1 , respectively, we obtain: [14] t XY n,n±1 (ρ, θ) = t XY where: is the transfer integral corresponding to the optimal face-to-face geometry (i.e., θ n,n±1 ≡ 0, ρ n = 0, ∀n) andη ≡ 1 + |η ppπ |/η ppσ = 1 + 2.26/5.27 1.429 [35]. In the B-DNA form equilibrium configuration (θ n,n±1 ≡ θ 0 , ρ n = 0, ∀n), Equation (8) reads: Sinceη > 0, we get t XY n,n±1 (0, θ 0 ) < t XY 0 , hence indicating that by explicitly considering helical geometry in the equilibrium configuration, the π − π base coupling strength is significantly reduced below that corresponding to the optimal face-to-face geometry. If we relax the equilibrium structure, allowing for the propagation of low frequency twist oscillations (acoustic modes), though keeping the radial variable describing H-bonding stretch oscillations fixed (no optical modes), Equation (8) can be approximated as: for small enough twists (cos θ n,n±1 1 − θ 2 n,n±1 /2), where the dimensionless parameter χ ≡η(R 0 /l 0 ) 2 2.92 measures the electron-phonon coupling strength. Despite its approximate nature, Equation (11) reasonably reproduces the main features of the transfer integral versus twist angle dependence derived from detailed quantum-chemistry calculations in the regime of low energies [12,14]. If we now allow for radial oscillations, still keeping within the small twist angle regime, Equation (8) adopts the form: where v n ≡ ρ n /R 0 . Finally, in the limit v n 1, Equation (12) can be approximated as: where only terms up to the second order are retained.

Electronic Hamiltonian
In order to obtain a realistic description of the rich dsDNA physico-chemistry, keeping at the same time the convenient mathematical simplicity, we will exploit the three-step renormalization approach sketched in Figure 3.
Version November 9, 2020 submitted to Materials 7 of 18

116
In order to obtain a realistic description of the rich dsDNA physico-chemistry, keeping at the same 117 time a convenient mathematical simplicity, we will exploit the three-step renormalization approach 118 sketched in Figure 3. In the fist step, the Watson-Crick bps present in the triplet codon shown in Figure 3a are renormalized to obtain the tight-binding model depicted in Figure 3b. The renormalized on-site energies and transfer integrals are respectively given by [36,37] : ε ij = t ij ⊥ = ε ji , which describes the charge carrier hopping from one base to its complementary base within the same bp [38], and are the effective transfer integrals between the bps and the sugar-phosphate groups, where t P is the effects [36,39]. 126 We note that the model depicted in Figure 3b can be properly regarded as a 3D generalization of the so-called fishbone model in 2D [40,41]. Accordingly, in the second renormalization step, the In the fist step, the Watson-Crick bps present in the triplet codon shown in Figure 3a are renormalized to obtain the tight-binding model depicted in Figure 3b. The renormalized on-site energies and transfer integrals are respectively given by [36,37]: ε ij = t ij ⊥ = ε ji , which describes the charge carrier hopping from one base to its complementary base within the same bp [38], and: are the effective transfer integrals between the bps and the sugar-phosphate groups, where t P is the glycosidic bond transfer integral, E is the charge carrier energy, γ j measure the sugar-phosphate groups' on-site energies, and ε j are the on-site energies of the corresponding nucleobases. In general, γ j will depend on the nature of the neighboring base, as well as the presence of water molecules and/or counterions attached to the backbone. Thus, the renormalized model parameters ε ij and τ j entail substantial physicochemical information concerning nucleotide interactions and backbone gating effects [36,39].
We note that the model depicted in Figure 3b can be properly regarded as a 3D generalization of the so-called fishbone model in 2D [40,41]. Accordingly, in the second renormalization step, the sugar-phosphate groups contribution is decimated [36,40,42], so that the original dsDNA molecule is mapped into the equivalent 1D binary lattice shown in Figure 3c, where the renormalized on-site energies labeled α and β correspond to the Watson-Crick complementary bps, and they are explicitly given by (E = γ j ): where t XY n,n±1 (ρ, θ) accounts for the aromatic base stacking between adjacent nucleotides, given by Equations (8)- (13). Therefore, the renormalized on-site potentialsε XY n (E) explicitly depend on the charge carrier energy, as well as on the angular and radial coordinates describing dsDNA oscillations, thereby enclosing all the relevant physicochemical information of the considered system. In this way, one obtains a realistic description, including 14 electronic model parameters, {ε j , t j , γ j , t GC ⊥ , t AT ⊥ }, fully describing CT throughout dsDNA molecules in terms of just two main functions, namely, XY n (E) and t XY n,n±1 (ρ, θ), in a unified way in terms of the effective 1D Hamiltonian: where c † n (c n ) is the creation (annihilation) operator for a charge at the nth site in the chain. In this way, the Hamiltonian given by Equation (16) provides a realistic treatment of CT mechanisms in dsDNA under physiological conditions, properly taking into account the influence of the dynamical state of the macromolecule on the CT efficiency. It is worth noting that epigenetic processes such as methylation (the addition of a methyl group (-CH 3 ) to one of the bases) will modify the on-site nucleobase energy and its effective mass alike, thereby changing both the mass ratio parameter λ among nucleotides and their γ parameter value. While the role of methylation-related on-site energy changes has been discussed in several recent works [43][44][45], the role of methylation-related dynamical effects in the CT efficiency has not. Due to the presence of both an on-site energy term ( XY n (E)) and a transfer integral term (t XY n,n±1 (ρ, θ)) in the diagonal term of the effective Hamiltonian given by Equation (16), one should expect the possible existence of resonance effects involving both electronic and dynamical physical parameters in an intertwined fashion.

General Expressions
From the lattice Hamiltonian given by Equation (2), we can straightforwardly obtain the canonical equations of motion: where P(ρ n ) = λ n ρ 2 n + 4m n M R 0 ρ n λ n ρ 2 n + 2R 0 ρ n and: along with: where: and as a first approximation, we assumed all the stacking stiffness parameters to take on the same value (i.e., k S n,n+1 = k S ∀n). In obtaining Equations (17) and (18), we neglected non-linear contributions related to theφ nρn andφ 2 n terms. Keeping only linear terms in the Taylor series of the functions appearing in the above expressions, we get the linearized equations of motion: where we reasonably assumed 2R 0 ρ n ξ 2 and introduced the twist frequency: along with the characteristic length l B ≡ 2 f 0 /g 0 , with f 0 = c 2 θ 0 + R 2 0 sin θ 0 0.759 nm 2 and g 0 = R 0 (1 − cos θ 0 ) 0.177 nm, so that l B 8.575 nm (i.e., about 25 bps) and: where a B ≡ g 0 ξ 2 / f 0 0.307 nm is a characteristic length whose value is comparable to the equilibrium bps separation h 0 . We introduced the coupled frequencies: and: where b B ≡ a B /ξ = g 0 ξ/ f 0 0.267 is a dimensionless factor. Therefore, ω ϕρ,n depends on the twist frequency, as well as the radial stretch H-bonding and lateral stacking oscillations of bps, whose characteristic frequencies are ω 2 H,n = D n α 2 n /(2M), and ω 2 S = k S /M, respectively. The site label in ω ϕρ,n arises from the presence of the λ n factor, as well as the fact that ω 2 H,n is site dependent, due to the different Morse potential parameters values for G≡C and A=T bps. On the other hand, ω ϕS,n±1 involves lateral stacking and twist oscillations. In this case, the site label dependence involves all the λ k terms. The set of coupled Equations (19) and (21) describes the dynamics of general dsDNA molecules, where two kinds of bps can be arranged either periodically or aperiodically [21,[46][47][48][49].

Dynamics of Homopolymer dsDNA Macromolecules
For the sake of simplicity, we will consider in this section the homopolymer case (i.e., polyA-polyT or polyG-polyC chains). In this case, the renormalized chain shown in Figure 3c becomes an effective monoatomic lattice (i.e., α ≡ β), where the renormalized on-site potentials depend on both the electron energy E and the phonon wavevector q, due to the presence of the π-π transfer integral in Equation (15). In addition, λ n = λ XY , D n = D XY ≡ D and α n = α XY ≡ α (see Table 1), so that the frequencies ω H,n ≡ ω H , ω ϕρ,n ≡ ω ϕρ and ω ϕS,n±1 ≡ ω ϕS are no longer site dependent. Hence, Equations (19)- (21) can be rewritten as (henceforth, we will drop the subscript XY in the λ parameters for the sake of clarity):φ where A λ ≡ (1 + λ)l −1 B and B λ ≡ a B (1 + λ −1 )/2 are constants, and we introduced the frequency ω 2 HS ≡ µ(ω 2 H + 2ω 2 S ), describing the coupling between the stacking and H-bond radial stretch oscillations, where µ ≡ λ −1 (1 + λ) 2 can take on two values. Making use of the model parameters listed in Table 1, we obtain the values listed in Table 2 for the characteristic frequencies just introduced, along with their related time scales. Table 1. Geometrical, dynamical, lattice model, and electronic model parameters adopted in the dsDNA homopolymer model studied in this work. The same effective Morse potential is used to describe H-bonding in both GC and AT bps. The spring constant k B is difficult to estimate, and different possible values, ranging from 0.04 to 0.5 eV Å −2 , have been reported in the literature [50][51][52][53][54][55].

Geometrical
Dynamical Lattice Electronic (eV) Electronic (eV) [25,56] t GC ⊥ = 0.01 [38,57,58]  From the data listed in Table 2, we see that the time scale of angular motions, determined by the twist frequency, amounts to ∼6 ps, which are an order of magnitude slower than those corresponding to the twist-radial coupled oscillations. The time scale related to bps' H-bond and stacking motions occupy an intermediate position, whereas coupled oscillations involving stretch and stacking motions are about 2.6 times quicker than the H-bond-mediated stretch oscillations alone. For the sake of comparison, the transition times reported for intrastrand hole transfer in ds-GT n GGG oligonucleotides range from τ = 0.5 ps for n = 1 to τ = 315 ps for n = 4 [63]. Quite interestingly, the electrical response of biological dsDNA chains to light irradiation has been recently investigated in order to engineer a DNA based molecular switch. In these experiments, it was observed that the electrical current turns on when the frequency of the incident light is above the 2 THz threshold [64], a value that coincides with that listed for ν ϕρ in Table 2. On the other hand, it is worth mentioning that the ν HS frequency value listed in Table 2 is smaller than the 2.83 and 3.04 THz frequencies experimentally observed by optical Kerr-effect spectroscopy in ds-GGCGGCCCGCGCGGGCCGCC and ds-ATTATTATTATATTA oligonucleotides, respectively. These oscillations are consistent with delocalized optical phonon modes with a wavelength extending throughout the molecule as a whole, and they are assumed to be related to the H-bonding dynamics between the two strands [65][66][67].
The mathematical structure of Equations (24) and (25) clearly indicates the correlated nature of next-neighboring bps dynamics. It is then convenient to zoom out our perspective and focus our attention on the dynamics of consecutive triplets of bps along the double-helix, which are closely related to the so-called codon units in genomics. To this end, we properly add up the dynamical equations corresponding to the consecutive bps corresponding to sites n − 1, n, and n + 1, grouping the resulting expression in terms of the collective variables x n ≡ 2ϕ n − ϕ n−1 − ϕ n+1 and y n ≡ 2ρ n − ρ n−1 − ρ n+1 , to obtain:ẍ Remarkably enough, we realize that Equations (26) and (27), describing the codon dynamics as a whole, are formally identical to Equations (24) and (25), which describe the motion of their constituent bps, since they are invariant upon the simultaneous variable exchange ϕ n ↔ x n and ρ n ↔ y n . This property can be regarded as expressing a renormalization of the dynamical equations when going from the bp local scale to the longer triplet codon scale. Quite interestingly, the very mathematical structure of the codon dynamics prescribed by Equations (26) and (27) guarantees that we will obtain similar dynamical equations by grouping codons in successive triplets of nested codon units recurrently, all the way up to the entire DNA molecule itself. Accordingly, the set of Equations (24) and (25) exhibits a self-similar symmetry upon triplet renormalization operation, so that by solving this fundamental dynamical equation set, we are actually disclosing the main features of the dynamics of the entire dsDNA macromolecule as a whole.
Inspired by previous results [14,15], we look for solutions to Equations (24) and (25) of the form ϕ n = ϕ 0 e i(ωt−nqξ) and ρ n = ρ 0 e i(ωt−nqξ) , describing a helical wave propagating throughout the dsDNA with frequency ω and wave vector q, where ϕ 0 8 • = 0.14 rad and ρ 0 0.05 nm are the twist and radial oscillation amplitudes at ambient temperature, respectively [68]. In so doing, Equations (24) and (25) can be expressed in the matrix form: The solution to Equation (28) requires the matrix determinant to identically vanish, thereby leading to a biquadratic equation whose solutions yield the dispersion relations for the acoustic and optical phonon branches given by: where G(q) ≡ 4ω 2 ϕ sin 2 (qξ/2), H(q) ≡ ω 2 ϕρ − ω 2 ϕS cos(qξ), and C λ ≡ A λ B λ . The exact dispersion relations given by Equation (29) are plotted in Figure 4, though they can be very well approximated by the simpler expressions: which extend previously reported results [25,30]. As we see, the acoustic branch is completely determined by twist oscillations, whereas the optical branch depends on both stretch and stacking oscillations. In particular, the q = 0 bandgap, ∆ν(0) ≡ ν + (0) − ν − (0) = √ µν H 1.70 THz ( 7 meV), is fixed by the H-bond frequency value. The maximum bandgap width occurs for q * = π/ξ 2.732 nm −1 , with ∆ν(q * ) = ν HS − 2ν ϕ 1.91 THz. The sound velocity obtained from the acoustic dispersion curve is 1.2 km s −1 , a figure smaller than the experimentally reported values ranging from 1.7 to 4.3 km s −1 , depending on the employed technique [69].
An interesting subject that can be addressed within the framework presented in this work refers to the behavior of the correlated oscillations when the DNA molecule is bonded with other small molecules. Generally speaking, the presence of ligands attached onto the sugar-phosphate backbone simultaneously affects the electronic properties and the mass distribution in a local region of the DNA molecule. This twofold effect can be properly accounted for in terms of a change in the γ value parameter, which modifies the electronic bandgap according to Equation (40) and a change in the values of parameters λ and M, thereby modifying the frequency values given by Equations (20)- (23). Thus, any effective bp mass increase leads to a slow down of frequencies ω ϕ , ω H , and ω S , along with their related coupled frequencies. In addition, the reduction of the twist frequency value makes the acoustic dispersion relation slope decline, so that the sound speed is reduced as well, ultimately leading to a lower thermal conductivity around the place where the small molecule has been bonded. On the other hand, a smaller ω ϕ value will widen the gap between the acoustic and optical branches (see Figure 4).  hand, a smalle ω ϕ value will wide the gap between the acoustic and optical branches (see Figure 4). Making use of the helical wave complex expressions for the variables ϕ n (t) and ρ n (t) into Eq. (13) we get t XY n,n±1 (ρ n , θ n,n±1 ) t XY where A 2 0 ≡ ϕ 2 0 + (ρ 0 /R 0 ) 2 0.02, so that 4χA 2 0 0.26 eV, and t XY n,n±1 > 0, ∀q. The resulting transfer integrals become site independent, which is a natural consequence of the synchronized motion mediated by helical waves propagating in the dsDNA molecule. Accordingly, we can write . Acoustic and optical phonon branches for polyG-polyC and polyA-polyT dsDNA homopolymers (their ν − (q) and ν + (q) curves respectively coincide up to the third decimal place), where q is measured in nm −1 . We used the characteristic and coupled frequency values listed in Table 2.

Charge Transfer through dsDNA Homopolymers
Making use of the helical wave complex expressions for the variables ϕ n (t) and ρ n (t) into Equation (13), we get: where A 2 0 ≡ ϕ 2 0 + (ρ 0 /R 0 ) 2 0.02, so that 4χA 2 0 0.26 eV and t XY n,n±1 > 0, ∀q. The resulting transfer integrals become site independent, which is a natural consequence of the synchronized motion mediated by helical waves propagating in the dsDNA molecule. Accordingly, we can write t XY n,n±1 (ρ n , θ n,n±1 ) = t XY (q) ≡ t XX (q) in the homopolymer case, and the Schrödinger nearest-neighbor tight-binding equation corresponding to the electronic Hamiltonian given by Equation (16) reads: where ψ n is the electronic wave function at site n. Equation (32) can be expressed in the matrix form: and within the framework of the transfer matrix formalism, the charge carrier dispersion relation of a dsDNA molecule containing N bps is given by (assuming periodic boundary conditions): where κ is the charge carrier wavevector. Since det M(E, q) = 1, we can use the Cayley-Hamilton theorem for unimodular matrices in order to calculate the required power matrix as M N (E, q) = U N−1 M(E, q) − U N−2 I, where I is the identity matrix and U m (x) ≡ sin [(m + 1)φ]/ sin φ, with x ≡ cos φ ≡ trM/2, are Chebyshev polynomials of the second kind, satisfying the recurrence relationship U m+1 − 2xU m + U m−1 = 0. In this way, we obtain [21]: Taking into account the relationship U m − U m−2 = 2 cos(mφ), the charge carrier dispersion relation can then be expressed as: According to Equation (15), the XY (E) function entails detailed information regarding the electronic structure model of the dsDNA molecule. X-ray experiments indicated that counterions condense around the nucleic acid chain in a tightly-bound layer, in agreement with early model calculations [70]. Therefore, a homogeneous charge distribution through the backbone can be assumed as a first approximation, that is γ j ≡ γ. In that case, we can write: where a XY ≡ ε X + ε Y and b XY ≡ (ε 2 X + ε 2 Y )/(2t 2 P ). Thus, plugging Equations (31) and (37) into Equation (36), we obtain the charge carrier dispersion relation in the explicit polynomial form: where ≡ γ(γb XY − a XY ) + t 2 P , and we introduced the auxiliary function: The resulting energy spectrum structure consists of two slightly asymmetric bands, E ± (κ, q), with relatively small widths (W ± ), separated by a gap, E g (0, q), whose value depends on the phonon wavevector q, as is illustrated in Figure 5 and Table 3 for the particular case γ = 0. By inspecting Figure 5, we can clearly appreciate that the phonon coupling gives rise to a systematic reduction of the bandgap width E g (0, 0) as q increases up to the value q * = π/ξ (see Figure 4), thereby enhancing the CT efficiency. Thus, the bandgap relative variation amounts to about 9% (6.5%) for polyG-polyC (polyA-polyT), respectively, as compared to the E g (0, 0) value. We note that the bandgap values listed in Table 3 are remarkably smaller than those usually reported in other studies [71][72][73][74][75][76]. Indeed, the dsDNA electronic structure is very sensitive to the precise value of the sugar-phosphate on-site energy. Thus, for the general case γ = 0, the bandgap width can be explicitly expressed as: so that a semiconductor-semimetal transition can be promoted by properly tuning the adopted γ value [36].
(polyA-polyT), respectively, as compared to the E g (0, 0) value. We note that the bandgap values listed in Table 3 are remarkably smaller than those usually reported in other studies [71][72][73][74][75][76]. Indeed, the dsDNA electronic structure is very sensitive to the precise value of the sugar-phosphate on-site energy. Thus, for the general case γ = 0 the bangap width can be explicitly expressed as so that a semiconductor-semimetal transition can be promoted by properly tuning the adopted γ value [36].  Table 1 and the choice γ = 0 eV. The wavevectors k and q are measured in nm −1 , and q * = π/ξ 2.732 nm −1 . The band structure of a polyA-polyT homopolymer is very similar (see Table 3) with the center of mass of the bands shifted upwards in energy by about 0.02. eV Table 3. PolyG-PolyC and PolyA-PolyT band structure properties (measured in meV). We have adopted the values γ = 0, a GC = 16.6 eV, a AT = 17.4 eV, b GC = 30.9, b AT = 33.6.

206
The physical picture inspiring this work relies on the fact that the presence of collective, 207 orchestrated oscillation motions of bps within the DNA double helix structure can efficiently enhance 208 the π − π orbital overlapping between bps along the backbone chain, hence promoting charge transfer 209 Figure 5. Electronic band structure for a polyG-polyC homopolymer. We used the electronic model parameter values listed in Table 1 and the choice γ = 0 eV. The wavevectors k and q are measured in nm −1 , and q * = π/ξ 2.732 nm −1 . The band structure of a polyA-polyT homopolymer is very similar (see Table 3) with the center of mass of the bands shifted upwards in energy by about 0.02 eV. Table 3. PolyG-polyC and polyA-polyT band structure properties (measured in meV). We adopted the values γ = 0, a GC = 16.6 eV, a AT = 17.4 eV, b GC = 30.9, and b AT = 33.6.

Conclusions
The physical picture inspiring this work relies on the fact that the presence of collective, orchestrated oscillation motions of bps within the DNA double helix structure can efficiently enhance the π − π orbital overlapping between bps along the backbone chain, hence promoting charge transfer via a long-range, phonon-correlated tunneling effect involving bps, which are relatively far apart. This property is intimately related to the helical geometry of the nucleobases' arrangement along the duplex chain, which makes possible that a local hopping process, involving a relatively small number of neighboring nucleotides, ultimately extends over the entire DNA chain as a consequence of the synchronized nature of the resulting helical wave. The possible presence of these helical waves may be relevant in the study of CT properties in dsDNA polymers exhibiting extensive chemically homogeneous regions along their strands, such as those reported in genomic studies of the telomere sequences of certain invertebrates [77] or in tandem repeats [78], involving mononucleotide triplet motifs (AAA, TTT, GGG, or CCC) [79]. Indeed, the existence of DNA-mediated charge migration has been related to the understanding of the damage recognition process, including the presence of lesions and mismatches. Furthermore, since CT dependence can be sensed electrically, it can be exploited for nanotechnological applications through base modifications or DNA-protein binding or with the task of designing nanoscale sensing of genomic mutations, opening new challenges for emerging nanobiotechnologies [80][81][82][83]. In addition, the fundamental dynamical mechanisms reported in this work are expected to also take place in other π − π molecular wires, such as G based quadruplexes [84].
Funding: This research received no external funding.
Acknowledgments: I sincerely thank Constantinos Simserides for his continued interest in my research works and Evgeni B. Starikov for illuminating conversations on the physico-chemical aspects of DNA molecules. I thank Victoria Hernández for a critical reading of the manuscript.

Conflicts of Interest:
The author declares no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: