Current-induced enhancement of DNA bubble creation

Current-induced heating of short double-stranded DNA chains is studied within a two-probe transport setup by using the Langevin approach. The electrons are modeled by a tight-binding Hamiltonian. The DNA atomic motion is described by the Peyrard–Bishop–Dauxois atomic potential, coupled with electrons through the Holstein interaction. The solvent environment is accounted for as a classical heat bath. Voltage biases of 0.1 ∼ 0.5 V can effectively break the base pairs and lead to the melting transition, which can be detected from the resulting significant reduction of the conductance. When the bias increases, the opening of base pairs near the leads with higher chemical potential is suppressed and bubble (localized separation of the double strand) formation becomes asymmetric. Our results suggest that the voltage bias can excite the base pairs, hence increases the chemical activity of DNA.


Introduction
The structure of double-stranded DNA (dsDNA) derives its stability from the hydrogen bonding between the base pairs and the stack interaction between neighboring bases. Under physiological conditions, the strength of the interactions for each base pair is of the order of k T B , where T is the bath temperature and k B is the Boltzmann constant, and thermal fluctuations can lead to local and transient separation of base pairs [1,2]. Breaking of a base pair weakens the stacking of the intrastrand adjacent bases and hence favors the breaking of the neighbouring base pairs. This cooperative opening of a sequence of consecutive base pairs ultimately leads to local denaturation in the closed double helix, i.e. DNA bubbles [3,4]. In solvent, the dsDNA is capable of closing the bubbles, and the fluctuating opening and reclosing of the double strand is so called DNA breathing. The breathing dynamics is dependent on sequence [5,6], length [7], bending rigidity [8] and supercoils [9] of the DNA chain. Besides solvent temperature, environmental factors such as salt concentration [5], pH value [10] and external stretching [11] can influence the bubble formation. While protein actions are usually discussed in terms of static structures, it is becoming increasingly clear that dynamics is essential for important biological functions [12][13][14], and it is meaningful to explore new ways of altering the dynamics of DNA.
Here we propose a new route to excite the base pairs and enhance the chemical activity of DNA. It is well known that Joule heating, i.e., the fluctuational part of current-induced forces in a molecule, is an important phenomenon in molecular electronics [15][16][17]. Another part of current-induced forces, i.e., the mean force, has been re-examined in recent years, and has a nonconservative component in out-of-equilibrium situations [18][19][20][21]. Its ability of doing net work along a closed path makes it a candidate for driving force of nanoscale engines. Nevertheless, the force can be destructive-it is able to continuously pump energy to the conductor and possibly leads to a runaway instability, resulting in damage and eventual failure of the conductor [22]. The nonconservative force may play a more important role in apparent heating of atomic nanowires [23,24], and is powerful enough to blow a molecular bridge connecting two metallic electrodes.
It is well known that electric current is an effective and controllable way to produce Joule heating in molecules and nanodevices, if the inelastic transport is taken into account. In this paper, we study the motion of the base pairs in some short DNA sequences under voltage biases to show how the creation of the bubbles in DNA molecules. Considering that specific solvent composition (water, ions) will influence the transport properties and dynamic process of the DNA sequences, a heat bath is included in our study to imitate the solvent environment. Our work is based on a framework of computing current-induced forces dealing with atomic motions which are slow compared to electronic time scale [21,22]. Electrical transport in DNA chain belongs to this regime because the bases are heavy groups (around 300 amu). We describe the nonlinear atomic dynamics by the simple potential model of Peyrard-Bishop-Dauxois (PBD) [25], which reduces the dsDNA to a onedimensional (1D) lattice with stretching of each base pair represented by a single coordinate. Its validity in characterizing the statistical and dynamical properties of DNA near denaturation transition has been demonstrated by several comparisons with experiments [26][27][28][29]. The pair stretching is coupled to the electrons through the Holstein interaction, where the on-site energy in charges is linearly dependent on the lattice displacement. The device setup is shown in figure 1.
There are two types of Waston-Crick base pairs, i.e., adenine-thymine (AT) and guanine-cytosine (GC). AT has two hydrogen bonds while GC has three, thus more energy is needed to open the GC pairs. Both experimental [30][31][32] and theoretical [33,34] works show that the sequences with high density of AT pairs have much higher probability to form the denaturation bubble. To have a conservative estimation of the breaking power of current-induced forces for general dsDNA chains, we consider the homogenous GC sequences at environmental temperature T = 300 K, for which the bubbles caused by thermal fluctuations are rare and transient presences [33]. Our results suggest that a bias of0.1 0.5V can effectively break the GC pairs and create denaturation bubbles. The result raises a proposal that voltage bias can be used as an effective route to excite the base pairs and enhance their chemical activity.
The paper is organized as follows. In section 2, we propose a Hamiltonian model to describe the DNA chain contacted with two leads, and then introduce the theoretical method to investigate the dynamic process of the DNA chain in a solvent environment with an applied voltage bias. The obtained results on the time proportion of a base pair, the conductance, the mean displacement and the time-averaged means forces to show the dynamic behaviors in the DNA chain are presented and discussed in section 3. Summary and concluding remarks are presented in section 4.

Model and method
We describe the electronic part of the DNA chain by using a 1D tight-binding model with each end coupled to a metal lead, as illustrated in figure 1(b represents the charge hopping between nearest neighboring sites, with = t 0.1 eV, a mediocre value within the range of the estimates for GC pairs (0.06 0.14 eV) [35,36] , where m a means the chemical potential in the lead α (=L, R). The tunneling between the chain and the leads is described by N is the length of the DNA chain. We assume the wide band   [37,38]. For the atomic motion, we focus on the transverse stretching of the pairs in the DNA chain. Thus, the atomic Hamiltonian is 2 . The potential energy part is described by the PBD model, including an on-site Morse potential and an anharmonic stacking interaction, Here, m is the reduced mass of the GC pair, y n represents the transverse stretching of the nth pair away from their equilibrium positions and the sum runs over all the base pairs. The Morse potential is a phenomenological description of hydrogen bonds in the same pair as well as other effective interactions between complementary nucleotides. The anharmonic interaction is essential for a better description of the dynamics of denaturation. When both y n and + y n 1 increase, the coupling between the nearest pairs becomes weaker and the energy decreases. This means that the interaction favors formation of denaturation bubbles. For homogenous chain of GC pairs, the parameters is set as 35 1 and the reduced mass m = 300 amu [25,27].
To study the dynamics of the GC pairs under the influence of the electrons, we adopt a simplified version of the semi-classical generalized Langevin equation (SGLE) approach [39][40][41][42][43][44][45][46]. Since the reduced mass of the GC pairs is much larger than that of the electrons, an adiabatic expansion of the electron influence functional is appropriate. The details of the derivation have been given elsewhere [17]. The SGLE, describing the dynamics of the DNA GC pairs, is written as˙( Here, the choosing value of γ is reasonable, since the preheating time allowing the DNA system to reach thermal equilibrium at each realization was 92 ps, and during this period the dissipation was 0.05 ps −1 [30]. Moreover, by adopting this parameter value, the thermodynamical properties from theoretical calculations are consistent with experimental observations [6,31]. Besides, each x n t is a white noise characterized by . The adiabatic force F n depends on the instantaneous atomic configuration of the DNA Hereafter, the limits of the energy integration is from -¥ to +¥ if not specified, Y is a vector made from the displacement of each pair in the system , with  = 1, and the average taken over the electronic density matrix at given configuration of the GC pairs Y. Green's functions in the time-domain and the frequency/energy domain are related through the Fourier transformation , ; e d t t i [47,48] (6), the average force on the n-th GC pair depends on the electron population, changing with the applied bias.
The electronic friction and noise are relevant beyond the adiabatic limit. To be consistent with the thermal bath, we take the zero frequency limit of the friction matrix, and ignore its frequency dependence. The friction matrix Here, we have dropped the arguments ( ) e Y ; , and ¶ e means derivative over energy ε. Given that the electronic fluctuations happen on much shorter time scales than that of the base pair motion, x n e is approximated by fluctuations locally correlated in time Here, the subscript s means taking the average over permutation of n and ¢ n , and m m = -eV R L is the applied bias. We have divided the noise into two parts: the first term in the square brackets of equation (9) is the equilibrium thermal noise, while the second is the correction due to the applied voltage bias. Thus, if eV = 0, Before closing this section, we make some comments on the approximation we use here. Firstly, when deriving the semi-classical generalized Langevin equation, the quantum effects of the thermal and electronic bath is taken into account, e.g., the zero-point motion, and the Bose distribution of phonons or electron-hole pairs [41]. But, to reduce the computation time, here we have taken the classical limit. Physically, this corresponds to take the time scale of the DNA atomic motion as the longest one in the whole problem. Thus, this approximation applies only when the applied voltage bias and the temperature are larger than the maximum frequency of the phonon system, well satisfied by the DNA atomic motion. Secondly, due to the simple Holstein type electron-vibration interaction we use here, the recently discussed current-induced nonconservative and effective magnetic force [19][20][21] are identically zero. Thus, we only look at the bias-induced renormalization of the atomic potential (equations (6), (5)) , and the stochastic Joule heating (equations (10), (11)).

Results and discussion
The Langevin equation (4) is numerically integrated using the stochastic Runge-Kutta algorithm of the second order [49], with a time step Dt of 1 fs and fixed boundary condition. The bubble formation is a fluctuational process of pair breaking and reclosing. To well characterize it, we compute the time proportion of a base pair staying in bubbles of different lengths, which is defined as The properties (i) and (iii) can be seen more quantitatively in figure 3, where we summed P n l over all bubble lengths ( =l 1 8) and also included the zero-bias result. When no bias is applied, even the center pairs have a total proportion less than 1.0% and the smaller proportion of the end pairs implies that most of the bubbles are not the full ones ( < l 8). For biases larger than 0.1 V, the values for the same pair are close to each other, except the difference near the left end (n = 1, 2), which clearly demonstrates the left-right asymmetry. Since at the biases of0.2 0.5 V the 8-bubble takes a large proportion and the close values of total bubble proportion indicate saturation of the bubble creation, we can infer that these biases have resulted in the melting transition, while the chain at = V 0.1 V is in a state between the closed strands and the melting phase.   figure 3, the second pair is the least likely to be broken and the sixth spends the longest time in denaturation. We show in figure 4 the stretching of these two pairs at V = 0.1 and 0.2 V as a direct checking of status of the states. It can be seen that the breaking of the pairs at 0.1 V exhibits an impulse pattern. The nearly synchronic opening of the two pairs implies that all the pairs tend to be broken at the same time and form an 8-bubble, which is consistent with the result in figure 2. So we find that the system at 0.1 V is indeed in a transition state: a denaturation phase can be reached but it is not sustainable. When the bias increases to 0.2 V, there are no clear interval and impulses, suggesting that the chain is in a steady denaturation phase. As the impulses are not regular and periodic, it is proper to assume that the bubble formation is due to stochastic fluctuation other than some ordered mechanisms like energy accumulation or periodic forces. This will become clear from the following analysis.

According to
It has been suggested that the mean stretching ( ) á ñ = å D y t yt Nt t N t n 1 1, , s s can well characterize the denaturation transition [50]. As shown in the inset of figure 5, the result for = V 0.1V is on the curve of sharp increase. When the biases are elevated to 0.2 V and above, a saturation is achieved, demonstrating that the states are in the denaturation phase. At = V 0.15 V there is a drop in the mean displacement, which appears to contradict with the expected monotonic increasing as biases are elevated. However, the monotonicity with voltage bias can not be assured since the current-induced mean force and the electronic fluctuation are two competing effects on the stretching of the base pairs. We can also see a smaller drop in mean displacement at = V 0.4 V. Since = V 0.15 V is in the critical region, the system is more sensitive to relative strength of the competing forces.
The transition can be detected by measuring change of the conductance, which has two orders of magnitude reduction from V = 0.05 to 0.2 V as shown in the inset of figure 5. The conductance is calculated by its definition of the current divided by the related voltage (I/V), in which the current is time average where G r and G a are temporal Green's functions calculated from configuration at time point t. Through the Holstein coupling, the stretching of base pairs affects the charge transport like fluctuating gates, which is positive most of the time because the PBD potential obviously favors positive displacement. Since on-site energy is zero in the electronic Hamiltonian, the transport channels lie around the Fermi energy (e = 0) when the displacements are zeros. In this case, there are always some channels within the bias window. In the denaturation phase, the effective gate voltage brings the channels upward, excludes them out of the bias window, and consequently leads to the dramatic decrease of conductance. As mentioned in section 2, the current-induced forces have two effects: (i) current-induced Joule heating and (ii), bias-induced modification of the effective atomic potential. In the following, we analyze these two effects and further understand the numerical results.
The bias-dependent electronic noise fluctuations, together with friction coefficient, give rise to the Joule heating effect. Actually, from equations (7), (10), we can define an effective temperature of the electronic bath at pair n as which is plotted in figure 6. We can see that T eff does not increase monotonically with the bias. From 0 to 0.1 V, the pair breaking is rare, and the displacements y n are around their equilibrium positions. The effective temperatures of all pairs increase due to the bias (see equation (10)). The hotter electronic bath is the energy source of bubble formation. But once the bubbles are formed at 0.2 V, T eff drops down. This is due to the fact that the electronic structure of the sequence is changed by the formation of bubbles. Further increasing of the bias leads to the increasing of T eff again. However, in this stage, the distribution of T eff within the chain changes drastically compared to the = V 0.1 V case. We see that T eff is largest at the site 6 and 7, which is well correlated with the distribution of time proportion in figure 2. It is noted that the distribution of T eff slightly breaks the left-right symmetry. This asymmetry, however, can not explain the asymmetry we observe in figure 2.
To shed light on the origination of the above asymmetry in the bubble formation, we turn to study the average force. Taking derivative of the electron Hamiltonian with respect to y n , it can be seen that the mean force is proportional to the mean occupation of electrons, † l = -á ñ F c c n n n , which always directs in the negative direction. Being position dependent, the mean force does no net work when the average is taken over an evolution long enough in time, and thus it functions like an effective potential bonding the base pairs and suppressing bubble creation. The average force may develop an asymmetry due to shifting of the left and right chemical potentials, as shown in figure 7. For the biases 0.4 V and 0.5 V, the asymmetry becomes violent: their values on the left and the right pairs can have a difference up to2 3orders of magnitude. Together with the consistency between its strengthening at the left end and the depressed breaking of the end pair, we believe that the asymmetry observed in figure 2 is ascribed to the average force.
Putting the above analysises together, we can get a consistent theoretical framework for the bias-induced bubble formation: the applied bias leads to an increase of electronic bath effective temperature T eff , which supplies the energy for the bubble formation in the DNA chain. Meanwhile, the average force develops an asymmetry due to shifting of chemical potentials. As a result, the probability of bubble formation at the left side of the DNA chain becomes smaller than that at the right.
From this picture, we may infer influence of altering the strength of the electron-phonon coupling and the thermal friction. As the dominant damping, it can be predicted that a weaker thermal friction leads to more effective bubble formation. Because both the mean force and the electronic fluctuation are proportional to the strength of the Holstein coupling, and the two forces play competing roles in the bubble creation, effect of changing l is less predictable. For a chain in the melting phase, the time proportion of a pair in broken states is nearly saturated and insensitive to variation of the parameters. To have a better comparison, we compute  Finally, to see the influence of the chain length on the bubble formation, we supplement the computations on the average stretching for the DNA sequence with 15 GC pairs in the inset of figure 8. For the case of N = 15, we find that the pairs near the two ends have a little smaller average stretching, which suggests that longer DNA chains are harder to separate. And meanwhile, no further suppression occurs in the center pairs, which actually have an average displacement much bigger than that of the chain with 10 pairs. This may be simply caused by the reduced anchor effect of the fixed boundary as the chain becomes longer.

Conclusions
We have simulated the dynamics of short DNA chains in a solvent environment with an applied voltage bias. A voltage bias as high as 0.2 V is capable of heating the DNA chain up to the melting phase, while for lower biases transient denaturation can occur but can not sustain. The current-induced force can break the base pairs and leads to the bubble creation in DNA chains. Due to the symmetry breaking in the chain by the unidirectional current, the bubble formation in DNA chains is not symmetric. Since the solvent environment is taken into account as a heat bath, the voltage bias can be applied as an effective route to create the bubbles in DNA chain, and to enhance the chemical activity of DNA molecules in experiments. Moreover, our numerical result of decreasing of the conductance versus elevated voltage biases agrees in order of magnitude with experimental measurements [51].
It is noted that our theoretical model is combination of the PBD potential and the tight-binding Hamiltonian to model charge transport along the π stack. For the PBD model itself, it does not take into account the helical structure of the DNA chain. If the model can be incorporated by an empirical potential, we can give a better description and the method can be applied to a wider range in low-dimensional physics. Moreover, for a real DNA chain, the conductivity of DNA chain is not only attributed by the π stack, but also influenced by the interface and solvent environment including water and ions. In our future study, we will develop the theoretical model and method towards real DNA chains, which helps us get more accurate dynamic behaviors of the DNA molecules.