A quantum mechanical model for hole transport through DNA: predicting conditions for oscillatory/non-oscillatory behavior

As the potential energy diagram depicted in figure 1 is a simplified 1D potential well, we applied the following Schrodinger Equation in 1D formulation:

$\frac{\partial }{\partial x^2}\psi =\frac{2m}{{\hslash }^2}(E-V)\psi$ (1a)

The corresponding equation for the probability current density is given by

$S=\frac{\hslash }{2im}(\psi *\frac{\partial \psi }{\partial x}-\psi \frac{\partial \psi *}{\partial x})$ (1b)

Here, m is the effective mass of the hole, ψ is its wavefunction for the energy state E, and V is the potential energy. The above equation is a steady state equation when there is no change in particle energy, and under a bound condition of figure 1 potential there cannot be any probability current S in any direction. If we assume that after the formation of a hole at the lone G site its wavefunction interacts with the A bases in the barrier and the remote GGG base cluster, the solution will no longer remain as a time invariant one. A time dependent Schrodinger Equation is a complicated one to solve, for which we applied the following simple mathematical treatment after the consideration of three different hole states (initial, intermediate and final states, discussed in section 2.4). If any portion of the wavefunction interacts with the vibrating atoms of bases (referred to as inelastic scattering) and changes its original energy state, the probability of the existence of the hole in that particular state will gradually decrease. In order to include this effect, an imaginary quantity is introduced in the (E-V) term of equation (1a) to yield eqn (2).

$\frac{{\partial }^2}{\partial x^2}\psi =\frac{2m}{{\hslash }^2}(E-V+i{V}_i)\psi$ (2)

Here, i is the imaginary quantity and Vi is the newly introduced term in the equation. The Vi term can be treated either as a component of an imaginary energy or an imaginary potential. Whatever the choice is, the final equation remains the same. It may be noted that the imaginary potential term Vi is a function of temperature [28,29] as the molecular vibrations or even chemical reaction rates can vary with temperature. So, it is possible to address the temperature variation by making Vi as a function of temperature. For any given energy E of the hole, the functional forms of the wavefunction in different regions (Figure 1) are given below:

$\psi =C{e}^{\gamma x}+ D{e}^{-\gamma x}$ (3)

Where γ = $i\sqrt{2m_1^{ }{}^{ }(E+i{V}_{i1})/{ћ}^2}$ for region PQ

= $i\sqrt{2m_2^{ }(E-V+i{V}_{i2})/{ћ}^2}$ for region QR

= $\tau i\frac{1}{2}$ for region RS

Here V is the height of the potential barrier and m1, m2 and m3 are the effective mass in PQ, QR and RS region (Figure 1) respectively. Vi2 and Vi3 are the imaginary potential in QR and RS regions respectively. Inelastic scattering in region PQ is ignored (Vi1 = 0). The usual solution technique of Schrodinger’s equation for such a potential structure involves application of boundary conditions at each potential discontinuity and evaluation of ‘C’s and ‘D’s for each region of a normalized wavefunction distribution. We used the technique presented in [30,31] to find out the eigenenergy and normalized wavefunction. The technique presented in [31] is particularly useful to evaluate normalized wave functions in presence of imaginary potentials. Once the wavefunction is known, the probability current density S (eqn.1b) can be easily calculated at any position using the form of solution presented in eqn. (3).

After the imaginary potential term is introduced within the A and GGG regions, the values of γ (eq.3) become complex resulting in gradual decay of the wave function in those regions due to inelastic scattering. If we represent the time constant of wave function decay as τ_i, then we can write

$\tau i=\frac{\underset{G}{\int }|\psi |{}^{2}dx}{S_Q}$ (4)

Here, S_Q is the probability current density at location Q (Figure 1) and has been defined by the number of holes being transferred per unit time. If such a process continues, after a while, the total wave function from the well at G will vanish and the hole will disappear from G and make a transition to GGG. No backward transition is considered as the initial energy at G is much higher than that at GGG [32].

Deriving equations for product concentrations at different base sites

Whenever a hole is formed at the lone G base, it has a very high energy barrier on its left hand side and the probability that the hole will cross that barrier is almost nil (Figure 1). The barrier on the right hand side is much smaller with a height of around 0.73 ev [26], and the hole may tunnel through that barrier or may have sufficient energy exchange with its neighboring bases due to thermal vibration or phononic interaction to overcome the barrier, and move to the GGG region. The tunneling through the barrier is a coherent process and there is no change in energy. On the other hand, if we consider an energy interchange with the surrounding molecules, there is a change in energy of the hole and a coherent process cannot be assumed. However, due to finite probability of both the processes, they would occur simultaneously and the ultimate result is the combined effect of both the tunneling and inelastic scattering. Hence, it can be concluded that a hole generated at the G site will die out from that state and have the probability of reappearing in some other state at some other location. It is a well-known QM phenomenon that the energy states go lower as the width of a quantum well gets larger. This situation is depicted in figure 1 by indicating E_G to be much higher in energy than E_GGG. Once the hole leaves its original energy state E_G, it will have a high probability of appearing in a lower energy state like E_GGG. As mentioned above, in calculating product concentrations at different base sites, we considered the following three hole-states that contribute to product formation.

Product formation during initial hole state

As soon as a hole is formed at the G site, its wave function (ψ₀) is primarily localized onto it with some penetration into the barrier on the right as shown in figure 2a. This is a metastable condition, and it is expected that the metastable wave function (ψ₀) will tunnel through the barrier (Figure 2b) and elastically transform to the wave function ψ₁ (Figure 3), which the wave function is corresponding to the whole potential structure. The time required for such a transformation is the time required for the accumulation of wavefunction ψ₁ in GGG. If the probability current density at R (Figure 1) for this state (ψ₀) is S_R0 then, the time constant (τ₁) for such a process to buildup ψ₁ in GGG is given by

${\tau }_1= \frac{{{\displaystyle \int }}_{GGG}^{ }|{\psi }_1{|}^{2}dx}{S_{R0}}$ (5)

If we assume that the probability of product formation rate at any location is directly proportional to the probability of hole existence (includes hole formation and decay) in that location and inversely proportional to the reaction time then

P_G0 = B α_G0 τ_C0 (1- exp(-τ₁ /τ_C0))/τ_GC (6a)

P_A0 = B α_A0 τ_C0 (1-exp(-τ₁ /τ_C0))/τ_AC (6b)

P_GGG0 = B α_GGG0 τ_C0 (1-exp(-τ₁ /τ_C0))/τ_GGGC (6c)

Here τ_XC is the reaction time in the region X (region for G, A or GGG in figure 1) and τ_C0 is the mean reaction time for the entire structure (Figure 1) calculated by taking the weighted average of reaction rates (inverse of reaction time) in different regions (please see Appendix 1 for a detailed derivation). The proportionality constant B takes into account any variation in the probability of existence of the wave function in the energy state under consideration and has the value of unity in this case as the hole is originated at this state. The α_Xj represents the probability of hole formation at a particular site and is calculated from the normalized wavefunction over the region X by applying the following eqn:

${\alpha }_{Xj}=\underset{X}{\overset{ }{{\displaystyle \int }}}|{\psi }_j{|}^{2}dx$ (6d)

Here, X represents the region (G, A or GGG) and j is the subscript for the wave function under consideration. If we consider X to be the region for G and set j = 0, then eqn. 6d gives a value for α_G0 obtained by integrating the square of the wave function ψ₀ over the region of G (PQ in figure 1). Other values of α, defined later in eqns 8 and 10, are also calculated in a similar manner. For a quick check of the formula 6a-c, when τ₁ or life-time of the hole is zero, P_G0= 0, that is, life time of the hole is too short for that energy state to form any product. On the other hand, when τ_GC (reaction time) is infinity or when the reaction takes place very slowly compared to other regions, again P_G0= 0. As expected, in this situation the hole does not exist at the base site long enough to form any product.

When the product formation probability due to ψ₀ at different base sites are known from eqns 6a-c, the total probability (for initial hole state) is given by

P₀= P_G0+ P_A0+ P_GGG0 (7)

Product formation during intermediate hole state

Now, the product formation probabilities at G, A and GGG sites together with the total probability (P₁) are given as follows, similar to eqn. 6a, 6b and 6c, for the hole state defined by ψ₁:

P_G1 = (1-P₀) α_G1 τ_C1 (1- exp(-τ₂ /τ_C1))/τ_GC (8a)

P_A1 = (1-P₀) αA1 τ_C1 (1-exp(-τ₂ /τ_C1))/τ_AC (8b)

P_GGG1 = (1-P₀) αGGG1 τ_C1 (1-exp(-τ₂ /τ_C1))/τ_GGGC (8c)

P₁= P_G1+ P_A1+ P_GGG1 (9)

As described earlier, α values are calculated by applying eq.6d with j = 1 for different regions of bases. Here (1-P₀) = B, as defined in eqns 6a-c and eqns in Appendix 1, and P₀₀.

Product formation during final hole state and total quantities

As a result of inelastic scattering the hole finally makes a transition to the lowest energy state, EGGG. Figure 4 shows the schematic potential energy diagram along with the probability distribution function for the final hole state represented by ψ₂.

P_G2 = (1-P₀-P₁) α_G2 τ_C2 (1- exp(-τ₃/τ_C2))/τ_GC (10a)

P_A2 = (1-P₀-P₁) α_A2 τ_C2 (1-exp(-τ₃ /τ_C2))/τ_AC (10b)

P_GGG2 = (1-P₀-P₁) α_GGG2 τ_C2 (1-exp(-τ₃ /τ_C2))/τ_GGGC (10c)

This is the lower most state of the potential energy structure, and the backward transition to higher states by inelastic scattering has been ignored as the next higher state is about 0.5 eV higher than this lowest state, and has a very small probability of occupancy. Hence, the value of τ₃ (time constant for hole transition to a higher state) is considered to be infinitely large in eqns. 10a-c. As before, α values are obtained by using eqn. 6d with j = 2. Hence, the total probability for product formation at each base site is given by

P_G = P_G0 + P_G1 + P_G2 (11)

P_A = P_A0 + P_A1 + P_A2 (12)

P_GGG = P_GGG0 + P_GGG1 + P_GGG2 (13)

In addition to giving the total amounts of products at G, A or GGG sites (eqns 11-13), this model also allows one to calculate product values within any specific location by taking into account the wavefunction distribution at that location. For example, we can calculate the distribution of products within the individual G locations of the GGG cluster. If GGG(1), GGG(2) (middle G) and GGG(3) are the three G locations within the GGG cluster (region RS as shown in figure 1), then the amounts of products formed at these individual locations are proportional to the ratio of the probability at that location (G) to the probability at the region encompassing the GGG cluster corresponding to each of the wave functions, ψ₀, ψ₁ and ψ₂, and can be expressed as

${\text{P}}_{\text{GGG}}\left(\text{r}\right){\displaystyle \sum }_j ={\text{P}}_{\text{GGGj}}\times ( \alpha { }_{GGGj}(r)/\alpha { }_{GGGj})$ (14)

where j = 0, 1 or 2 (corresponding to the subscripts of the wave functions ψ₀, ψ₁ and ψ₂) and ${\alpha }_{GGGj}(r)=\underset{GGG(r)}{\overset{ }{{\displaystyle \int }}}|{\psi }_j{|}^{2}dx$ , where r is the G location within GGG (Figure 1) and assumes values 1,2 or 3 and ${\alpha }_{GGGj}=\underset{GGG}{\overset{ }{{\displaystyle \int }}}|{\psi }_j{|}^{2}dx$ . Since the calculation of α requires ψ₀, ψ₁, and ψ₂ functions, in next section we discuss how to calculate these functions and their associated energies.

Applying QM model for calculating ψ₀, ψ₁, and ψ₂

As mentioned above, these calculations were done numerically by using a method described in [30,31] and using the MATLAB series of programs [33]. In the calculation of ψ₀ we assume that soon after the formation of the hole at G, it will not experience the potential due to GGG, for which we consider a potential profile with an infinitely extended A barriers to the right as shown in figure 2a. With the barrier height [22] of 0.73 eV, the lowest eigenstate for the G well (Figure 2a) has the energy of 0.64 eV. Then, we assume a potential structure like that shown in figure 2b and calculated the rate of tunneling through the barrier of A bases. After tunneling, the wave function is distributed as shown schematically in figure 3, and is represented by ψ₁. Since tunneling is a coherent process both ψ₀ and ψ₁ will have the same energy (0.64 eV). Although ψ₀ belongs to the lowest energy state for the potential structure of figure 2a, ψ₁ does not for the potential structure of figure 3 or 4. The lowest eigenstate (ψ₂) for the potential structure of figure 3 or 4 is substantially lower in energy with a value of 0.16 eV (E_GGG). For this large energy gap, the possibility of back transfer of hole from the low energy GGG site to high energy G site has been neglected in our numerical calculation. Even though for illustration purposes we provided schematic diagrams of wavefunctions (Figures 2-4), some of the actual functions, as calculated from the potential structure are presented later (Figure 7a-c).

We also examined other immediately higher energy states (which ranges from 0.565 to 0.612 eV) for different barrier widths, and find that they contribute very little toward the hole population or product concentration due to their lower probability of occupancy, and hence, these excited states were ignored. It is important to know that the hole populations at different eigenstates are determined by the Maxwell Boltzmann probability distribution function.

Calculating current components for tunneling and inelastic scattering

In calculating tunneling and scattered current components, we first assumed V_i = 0 in the barrier region (but not in GGG) so that there is no inelastic scattering due to the barrier, and all the current that is observed at the boundary of G and (A)_n is due to tunneling current (S_T) only. Then, we introduced V_i in the barrier and calculated the current at G. This current includes both tunneling (S_T) and inelastically scattered (S_I) components. Thus, we separated the tunneling and scattered components from the total current and present results in table 1. Figure 5 shows how the ratio of tunneling to inelastic current density (S_T/S_I) decreases with an increased barrier width.

Effective hole mass, imaginary potential and reaction time

Because of the lack of reported results on effective mass or imaginary potential in the literature, we obtained these values by fitting the experimental data [13] points. These parameters were needed for the predictions of oscillatory and non-oscillatory nature of the product yield. The best fit curve and the data points are shown in figure 6. The effective masses of the hole at the G (or GGG) and A sites are determined to be around 1.42 m₀ and 1.50 m₀ respectively, where m₀ is the rest mass of an electron. These masses are in line with the typical hole masses that are considered in different charge transport phenomena in physics and electronic engineering. Similarly, the imaginary potential (V_i) value has been determined to be 0.25 × 10^-10 eV.

As discussed above, in the process of hole transfer from G to GGG, if the hole remains at a particular base site longer than the time required for its reaction with water, products will form. For the water trapping reaction with a hole at the G site (forming G+), the experiment by Giese and Spichty [34] provides a rate constant (k) value of around 6.0 × 10⁴/s and hence, the time of reaction (τ_GC) of 1.67 × 10^-5 sec (1/k). The same experiment also suggests reaction at the GGG site to be 3 times faster, and hence, τ_GGGC = τ_GC/3. While the experiment suggests a much slower reaction at the A site, no rate constant value has yet been reported for the water trapping reaction at this site. In the absence of such experimental data, we examined a number of values by considering rate at A to be 5, 10 or 15 times slower than that at the G site. These results suggest that 10 times slower rate or τ_AC = τ_GC × 10 gives a product concentration at A in better agreement with that of the experiment [13]. It is to be mentioned that only the percent product formation inside the barrier is dependent upon this value and has no effect on the overall findings of this study.