A New Mathematical Model for Calculating the Electronic Coupling of a B-DNA Molecule
Dale J. Igram1, Jason W. Ribblett2, Eric R. Hedin1, Yong S. Joe1, *
1Center for Computational Nanoscience, Department of Physics and Astronomy, Ball State University, Muncie, IN, USA
2Department of Chemistry, Ball State University, Muncie, IN, USA
To cite this article:
Dale J. Igram, Jason W. Ribblett, Eric R. Hedin, Yong S. Joe. A New Mathematical Model for Calculating the Electronic Coupling of a B-DNA Molecule. American Journal of Physical Chemistry. Vol. 5, No. 2, 2016, pp. 17-25. doi: 10.11648/j.ajpc.20160502.11
Abstract: The charge transport properties of DNA have made this molecule very important for use in nanoscale electronics, molecular computing, and biosensoric devices. Early findings have suggested that DNA can behave as a conductor, semiconductor, or an insulator. This variation in electrical behavior is attributed to many factors such as environmental conditions, base sequence, DNA chain length, orientation, temperature, electrode contacts, and fluctuations. To better understand the charge transport characteristics of a DNA molecule, a more thorough understanding of the electronic coupling between base pairs is required. To achieve this goal, two mathematical methods for calculating the electronic interactions between base pairs of a DNA molecule have been developed, which utilize the concepts from Molecular Orbital Theory (MOT) and Electronic Band Structure Theory (EBST). The electronic coupling characteristics of a B-DNA molecule consisting of two Guanine-Cytosine base pairs have been examined for variation in the twist angle between the base pairs, the separation between base pairs, and the separation between base molecules in a given base pair, for both the HOMO and LUMO states. Comparison of results to published literature reveals similar outcomes. The electronic properties (metallic, semi-conducting, insulating) of a B-DNA molecule are also determined.
Keywords: Hückel Method, Slater-Koster Relations, Atomic Orbitals, Electronic States, Overlap Integral, Bond Integral
The B-DNA molecule is considered to be a potential building block for molecular electronics due to its self-assembly and self-recognition properties. The field of B-DNA electronics is highly interdisciplinary, merging physics, biology, chemistry, computer science, and engineering. Therefore, this new field shows great promise in regards to using individual DNA molecules for producing a new range of electronic devices such as nanoscale electronics [1, 2], molecular computing [3, 4], and biosensoric devices [5, 6], that are much smaller, faster, and more energy efficient than the present semiconductor-based electronic devices. Also, because of its self-assembly and self-recognition characteristics, DNA can easily adopt to various states and conformations, thereby providing the possibility of producing nanostructures with very high precision, beyond what is achievable with traditional silicon-based technologies . Equally important is the understanding of charge transport in DNA in relation to damage and mutation throughout the macromolecule [7-9], the detecting, manipulating, and sequencing of DNA [10-12], and the transport properties of other systems with interactions, such as molecular crystals and discotic materials [13, 14].
Several researchers have suggested theoretically or shown experimentally, that the B-DNA molecule has electrical conducting properties. As early as 1959, Eley  proposed that a DNA molecule might behave as a one-dimensional aromatic crystal and illustrated electron conductivity along the helical axis. His proposal was based on the results of the electrical conductivity of crystalline organic substances (- electron compounds). There are many substances that behave as semiconductors with an energy gap that decreases with the number of mobile -electrons in the molecule, thus suggesting that conductivity is associated with the intermolecular tunneling of thermally excited - electrons. By 1960, Ladik  showed that the sigma coupling could be sufficient to allow for conductivity if some of the bases were in an excited or ionized state where the one-electron orbital overlap integrals for the sigma coupling between DNA base molecules are perpendicular to the helical axis. About the same time as Ladik, Pullmann and Pullmann  demonstrated through the use of Molecular Orbital Theory calculations that the Guanine-Cytosine (GC) base pair would be a better electron donor and electron acceptor than the Adenine-Thymine (AT) base pair. In 1962, Eley and Spivey , experimentally showed that interactions of the stacked base pairs in DNA could lead to conducting behavior.
The interest in charge transfer in DNA took off in the early 1990s when Barton and Turro suggested that ultra-fast photo-induced charge transfer can occur over large distances between donors and acceptors that are inserted in the DNA [19-21]. Their hypothesis sparked a wide variety of experimental and theoretical studies into the nature of charge migration in DNA. Several charge transport experiments [22-31] have been performed on single DNA molecules revealing that different electrical characteristics exist: insulating [25, 26], semiconducting [27, 28], ohmic [29, 30], and superconducting . The variation in the charge transport properties may be due to the high sensitivity of charge propagation in DNA to extrinsic (interaction with hard substrates, metal-molecule contacts, aqueous environments) as well as intrinsic (dynamical structure fluctuations and base pair sequence) influences. This paper will be focusing only on the charge transfer in a DNA molecule which is influenced by intrinsic properties, in particular, the dynamical structure fluctuations and electronic coupling between base pairs.
Two mathematical methods required for the calculation of the electronic coupling between the base pairs of a DNA molecule were developed. The first method utilizes Molecular Orbital Theory, specifically, Linear Combination of Atomic Orbitals (LCAO) overlap integrals to calculate the bond integral parameter k values, which are then used in a highly modified version of the extended Hückel method to determine the molecular orbital wave function coefficients. This more generalized form of the extended Hückel method was developed by essentially eliminating most of its assumptions and some of its limitations, resulting in more accurate values for the coefficients. The second method utilizes an amended version of the Slater-Koster relations from Electronic Structure Theory to acquire the necessary interatomic matrix elements. This amendment was required due to the non-covalent interactions involved between the B-DNA base pairs. Only the LCAO coefficients for the frontier (HOMO and LUMO) molecular orbital wave functions were considered. With known LCAO coefficients and interatomic matrix elements, the electronic coupling parameter for a two base pair B-DNA molecular system was ascertained.
The electronic coupling parameter varies significantly as the DNA structure changes due to mechanical influences. These changes in are illustrated and discussed for variations in the twist angle between the base pairs, the separation between base pairs, displacement of one base pair with respect to the other in the and directions, and displacement between base molecules for one base pair with respect to the other in the and directions. There exist similarities between the results for the variations in the twist angle and separation of distance between base pairs and those of published literature . The reasons for discrepancies between the values presented in this paper and those of published literature  are discussed.
The electronic coupling between two successive base pairs  is a summation of the products of the LCAO molecular orbital coefficients, , of both base pairs and the interatomic matrix elements, , between the base pairs. The coefficients are calculated utilizing a newly developed method, which is a more generalized method than the ordinary extended Hückel method, due to the elimination of most of its assumptions and some of its limitations, and will be referred to as the modified-extended Hückel (meH) method. The interatomic matrix elements are typically determined from the Slater-Koster relations , but for this work these relations were amended to include both attractive and repulsive terms in order to more accurately calculate the non-covalent interactions that exists between B-DNA base pairs.
The B-DNA model used in this article consists of two stacked GC base pairs. Each base pair contains nineteen atomic locations, which is characterized by the general molecular orbital energy equation
where represent the Hamiltonian energy operator, E the total energy, the overlap integral, and the coefficients of the molecular orbital wave function. A secular determinant can be created by utilizing Eqn. (1). If the molecular system consists of other atoms besides carbon, then the secular determinant needs to be modified to reflect the fact that the different atoms will have a different number of valence electrons. Typically, adjustment parameters and are used, which are obtained empirically or theoretically. For the B-DNA structure, values from Table 8-3 in Ref.  were used; however, the parameter was determined using the meH method because its value depends upon the distance between two atomic orbitals.
By definition, where is the measured interaction energy between two atomic orbitals, and represents a standard for the benzene bond distance (1.397Å). values can be difficult to obtain; thus, the overlap integral can be used instead as it is determined theoretically. It has been proposed that , where is a non-energy quantity. The values of the parameter for all the combinations of paired atomic orbitals in the B-DNA molecule were ascertained using . The equations for for all Slater type atomic orbital pairs consisting of and atomic orbitals for n=1,2,3, and 5 have been formulated . It is well known that only the atomic orbitals contribute to the electronic coupling; thus, only the equations for and atomic orbital pairs are considered. A atomic orbital pair is realized when one atomic orbital is aligned parallel to but directly above the other atomic orbital, and a atomic orbital pair is created when one atomic orbital is parallel to but coplanar with the other atomic orbital (Figure 1).
The andatomic orbital pairs equations are expressed in terms of two parametersand, which are functions of the Slater values, inter-nuclear distance, and Bohr radius. Thus, we have for
where the parameters and are defined as
for and Equations (6) and (7) were obtained from published literature [37, 38].and are described by and , whereandare Slater values, inter-nuclear distance, and =Bohr radius (0.529Å). The coefficients were only determined for the frontier molecular orbitals (HOMO and LUMO).
Applying the LCAO method to solids in a rigorous manner is quite complicated. However, treating the LCAO method as an interpolation process allows simplifications to be made  which permit the potential energy in the Hamiltonian to be treated as a sum of spherical potentials located on the two atoms on which the atomic orbitals are located. The energy matrix components of the Hamiltonian operator are defined using
where represent the positions of the atoms for which orbitals are located, the positions of the atoms for which the orbitals are located, the displacement vector, and and as Löwdin orthogonalized functions (non-orthogonal functions that were made orthogonal by the Löwdin method) with the same symmetry properties as the atomic orbitals for a crystal. Each function can be expressed as a sum of the angular momentum components () with respect to , if the displacement vector defines an axis of a diatomic molecule. The quantities and represent the atomic orbitals on atom and at and , respectively. For atomic orbital functions, can be expressed as a linear combination of a and a function with respect to that axis.
The nature of these integrals can be determined by rotating axes and transforming spherical harmonic functions with one set of axes into spherical harmonic functions with another set of axes. Thus, to represent the integrals in Eq. (8), one needs to include contributions involving the product of an atomic orbital of this type on an atom at , another on an atom at , and the spherical potentials that are centered on these two atoms. By defining the direction cosines of the displacement vector as , where , and , Eq. (8) can be expressed in a form such as , similar to in Table I of Ref. . In this form, the function is a like function and another like function, thus allowing this new function to be written in terms of two integrals consisting of a orbital on the first atom and aorbital on the second, as well as a orbital on the first and a on the second. These integrals are symbolized by and the second by, where it is understood that the first and second indices represent the first and second orbitals, respectively. The third index represents the type of interaction or coupling that exists between the two orbitals. The integrals are functions of the distance between the atoms, thus the integrals have different values for different atomic pairs. As discussed previously, Eq. (8) can be expressed as
Because a B-DNA molecule consists of aromatic base molecules, there exists coupling or interactions between the atomic orbitals of one base pair with those of another base pair. When two orbitals are coupled, there can exist and interactions, which are characterized by the hybridization matrix element . The hybridization matrix elements used in Harrison’s model  were formulated to work well within a typical covalent bond distance but not for non-covalent interactions nor covalent interactions at distances greater than 2Å. This difficulty can be resolved by changing Harrison’s term to a standard term and by incorporating attractive and repulsive terms, which are used to scale the non-covalent interactions [33, 41, 42]. The attractive and repulsive terms [40, 41] are defined as and , respectively, where is the distance between two atomic orbitals, the typical covalent bond distance for a base molecule, an adjustment parameter, and the average of the van der Waals radii for two coupled atomic orbitals. Incorporating these terms changes the hybridization matrix elements from , to
Substituting Eqs. (10) and (11) into Eq. (9) yields
where is defined as the energy integral/interatomic matrix elements for two orbitals, the electron’s mass, the distance between two adjacent base pairs, and and are values obtained from a Density Functional Theory code like SIESTA.
The electronic coupling between two adjacent base pairs represents the strength of the coupling or the potential energy between two base pairs and is written as
where and represent the and orbitals for base pairs 1 and 2, respectively, the LCAO coefficient of base pair 1, the LCAO coefficient of base pair 2, the interatomic matrix elements for the and orbitals for base pairs 1 and 2, and and the total number of orbitals of base pairs 1 and 2, respectively. The difference between and is that represents the interaction energy between two atomic orbitals regardless of the type of molecular orbitals involved; whereas, represents the interaction energy between two atomic orbitals with the type of molecular orbitals considered.
Because a Guanine-Cytosine (GC) base pair contains 19orbitals, and there are two stacked base pairs, there are a total of 38orbitals creating a 19 x 19 matrix for . Each makes a small contribution to the electronic coupling . Equation (12) shows that will vary depending on how well the base pairs are stacked. In other words, by changing the B-DNA structure, such as by twisting the base pairs with respect to each other or varying the distance between the base pairs, may change significantly. Hence, it is very important to understand the dependency of regarding these conformational changes because it will affect the charge transport characteristics of the B-DNA molecule.
3. Results and Discussion
In this section, the results for a B-DNA molecule consisting of two GC base pairs will be presented along with a discussion. The data pertaining to the electronic coupling for the HOMO and LUMO energy states are provided as a function of the base pair’s twist angle, separation distance between base pairs, displacement of one base pair with respect to the other base pair in the x- and y- axes, and displacement of the base molecules within a given base pair with respect to the other base pair also in the x- and y- axes. In addition, the results representing only the electronic coupling as a function of the base pair’s twist angle and the separation distance between base pairs are compared to previous work , which reveals similar outcomes. The slight differences in the calculated versus the published data are discussed.
The electronic coupling is very dependent on the motion of the base pairs with respect to each other (Figure 2), where varies significantly for twist angle to. To determine the results, the following parameter values are used: the separation distance between the two base pairs of z = 3.375Å, a typical bond length of the B-DNA bases of Å, , an average of the van der Waals radii for C, N, and O atoms of Å, and adjustment parameter values C of eVÅ4 for HOMO and eVÅ4 for LUMO.
Figure 2 also reveals similar outcomes between the calculated and published values. The variation in is believed to be from the interaction matrix elements, because of their strong dependence on the geometry of the base pairs . In other words, as the twist angle changes, so do the and interactions. There are two basic reasons for the reduction of the electronic coupling: 1) there are positive and negative interactions occurring between two interacting atomic orbitals which could reduce or almost cancel each other, thus causing a smaller net atomic pair interaction and 2) there also are predominately large and pair interactions that are capable of reducing or canceling each other when summed up to calculate the total base pair coupling. Therefore, small individual atomic interactions, or an equal number of large positive and negative ones, can produce rather small base pair coupling. All the curves have one thing in common, which is that at only the interactions exist. As continues to increase, the interactions become less dominant with the introduction of some interactions, until goes to zero at specific values of as shown by the four curves (Figure 2). When becomes negative, then the interactions dominate with very little interactions occurring. An interesting note is that for the equilibrium twist angle of , the calculated LUMO state electronic coupling has a negative value as compared to the published value.
The electronic coupling also varies as a function of distance between the base pairs (Figure 3). Similar trends between the calculated and published data appear, with the exception of the published values for the LUMO state. The twist angle is held constant at with the other parameter values remaining unchanged, except for the separation distance.
There are numerous reasons for these discrepancies (Figures 2 and 3). Firstly, according to Ref. , a combination of a DFT code (SIESTA), which utilized a Double – Zeta Gaussian with Polarization (DZP) basis set, and a parameterized Slater-Koster model was utilized in generating the published values. The Slater-Koster parameters , and (0.87Å for B-DNA) were obtained by fitting to DFT information. Secondly, a Complete Neglect of Differential Overlap (CNDO) approximation was incorporated in order to simplify calculations by reducing the number of electron repulsive integrals to be calculated. Thirdly, a simple Hückel method was considered.
Whereas the model that was used to produce the calculated values consists of a combination of a generalized form of the extended Hückel method and a highly modified Slater-Koster method. This combination uses Coulomb and bond integral parameters, and attractive and repulsive exponential terms. In addition, no exponential cutoff distances and CNDO approximations were required due to the inclusion of the repulsive term used in Eqs. (10) and (11), which represents contributions from the electron-electron Coulomb repulsion and the ion-ion Coulomb repulsion. By not requiring an exponential cutoff distance, the exponential tails of the orbital wave functions extend to infinity, which means that the new model may provide more accurate approximations of the effects on electronic coupling due to additional interactions occurring at larger distances. In Ref. , the results pertaining to changes of the electronic coupling associated with rotation and separation between two base pairs only were considered, and any motion related to translations between two base pairs and between base molecules was ignored. In this paper, the additional degrees of freedom are contemplated and the data is discussed below.
It was found that the electronic coupling is affected in a somewhat unpredictable manner, by displacing the top base pair with respect to the stationary bottom base pair in the and directions (Figures 4 () and 5 ()). The axis is considered to be in the middle of the base pair but perpendicular to the line between the backbones, the axis in the middle of the base pair and parallel to the line between the backbones, and the axis along the helical axis. Because in Figure 4, the maximum peak for all the curves is at and , which is indicative of only interactions. The HOMO state reaches its maximum at and the LUMO state at Deviating from this point, the curves begin to decrease as a result of an increase in interactions. For the HOMO and LUMO states in the direction, does not vanish or go negative until the shift is beyond Å; however, this is not true for the HOMO and LUMO states in the direction. Also, for the HOMO and LUMO states in the direction is less sensitive to translations as compared to the HOMO and LUMO states in the direction (Figure 4). However, at the equilibrium conditions (separation distance of Å and twist angle of ) remains negative for all the curves, except for the LUMO state in the direction, where it becomes positive around Å and just slightly positive between and Å (Figure 5). This indicates that interactions dominate under these conditions. For and , the HOMO state has an electronic coupling value of and the LUMO state a value of . Another interesting item to note is that at about Å all four curves intersect at The variation of the electronic coupling for the HOMO state in the direction is less sensitive than for direction. In general, the results (Figures 4 and 5) correspond to interactions dominating and interactions dominating, respectively.
The effects on the electronic coupling associated with the base molecules in base pair 2 translating in opposite directions with respect to each other and to a stationary base pair 1 are presented (Figures 6 and 7). For the HOMO state at , has a value of and with respect to the LUMO state a value of (Figure 6). Similarly, the HOMO state has a value of and for the LUMO state a value of at (Figure 7).
There are similarities in the profile between the HOMO states for the and directions, and LUMO states for the and directions (Figures 4 and 6). Also, there exist similarities in profile when the twist angle is set to 36˚ (Figures 5 and 7). Additionally, the last four figures have revealed some intriguing results, which have provided further insight into the variation of the electronic coupling. The twist angle that exhibited minimum variation of the electronic coupling in regards to both the translation in the and directions of base pair 2 and between the base molecules in base pair 2 is (Figures 8 and 9). Because has mostly negative values it is indicative of interactions dominating, and that, relatively speaking, the electronic coupling between the two base pairs in the and range are the same (Figures 8 and 9).
In general, when the electronic coupling is at a maximum, whether due to or interactions, and a significant amount of wave function overlap exists, the final electronic state is a metallic electronic state. If the electronic coupling vanishes, then an insulating electronic state exists. When the electronic coupling is between these two conditions, there exists a semiconducting electronic state. Therefore, because a DNA molecule is very dynamic, it is conceivable that for a long DNA molecular chain, these three electronic states can occur simultaneously along the length of the chain.
The theories that were utilized for the calculation of the electronic coupling for holes and electrons that propagate from one base pair to another were presented. In particular, Molecular Orbital Theory (MOT) and Electronic Band Structure Theory (EBST) were explained. With regard to the MOT and EBST, two new mathematical methods were presented and shown to provide good values, which were in agreement with published literature [32, 33]. The slight differences in the calculated values versus the published, due to a difference in methodology, were discussed. Additionally, the equation that was utilized in the calculation of the electronic coupling between two adjacent base pairs was provided.
In this paper, the electronic coupling results for a B-DNA molecule consisting of two GC base pairs were discussed. These solutions were presented as functions of the base pairs’ twist angle, separation distance between base pairs, displacement of one base pair with respect to the other base pair in the x- and y- axes, and displacement of the base molecules within a given base pair with respect to the other base pair also in the x- and y- axes, for both the HOMO and LUMO states. The outcome of this work revealed several attributes: for , only the interactions exist; at , a mixture of and interactions occurred with the interactions being dominant; there are values where no electronic coupling can occur; for , as the distance between the base pairs decreases the interaction dominance increases; at , for all translations interactions prevailed; at , interactions dominated for all translations; and the interactions for have similar strength for all states. In summary, the electronic coupling is quite sensitive to changes in the twist angle and separation distance between the base pairs, but less sensitive to changes in the translation directions, especially at a twist angle of .