Charge diffusion in homogeneous molecular chains based on the analysis of generalized frequency spectra in the framework of the Holstein model

Full Text

1. Introduction Kinetics of charge transfer in polarized one-dimensional chains at finite temperature is an attractive and pressing theoretical problem by itself [1]-[4]. Its particular importance is related to the discovery of conducting properties in DNA [5]-[7]. On the first hand, this transfer may determine the biological functions of DNA. On the second hand, this opens up fresh opportunities in nanobioelectronics [8]-[10]. A basic qualitative simulation of charge transfer in DNA is a pioneered Holstein model for describing the charge transfer process in one-dimensional chains [11]. In the semi-classical approximation, DNA nucleotide pairs are modeled by a sequence of unconnected oscillators arranged in a chain [12]. A charge travels along a chain in the strong coupling approximation. Motions of the sites influence the charge propagation, and conversely, the charge probability density affects the charge motion. The Holstein model is rather simple and can hardly describe the charge propagation in DNA realistically. Nevertheless, it is nonlinear and demonstrates different regimes. The study of this model can give qualitative insight into the nature of charge transfer in one-dimensional chains, including biological polymers. The Holstein model is thoroughly investigated. Numerous papers are devoted to the analytical and numerical analysis of the original quantum model and its semi-classical approximation. However, the charge transfer kinetics at finite temperature is still to be understood. Of importance in this respect is the diffusion coefficient of a charge. Earlier, we studied its temperature dependence inhomogeneous chains of different nucleotide structures and found the scales which mainly determine differences in the charge diffusion in different chains [13]. This work continues the earlier started analysis. In order to better appreciate the kinetics of an added charge propagation along an infinite homogeneous chain, find all the motions which contribute into its mobility and determine the diffusion coefficient, we calculate and analyze autocorrelation functions of the charge propagation velocities and a relevant generalized frequency spectrum which is a frequency-dependent diffusion coefficient. We show that the earlier suggested scale is natural for the model, and the diffusion coefficients differ only in the low-frequency range for different sequences. In the limit of high and low temperatures, the nucleotide structure is irrelevant. Analysis of autocorrelation functions suggests that the charge propagation is contributed by dissipation and reflection at long distances as well as reflections from neighboring sites. As the temperature grows, these motions arise, coexist, and alternate. In the limit of finite temperature, the charge diffusion ceases. A detailed straightforward analysis of the Holstein model kinetics at finite temperature requires direct numerical simulation concerned with a vast amount of computations. Determinate time dependencies should be identified from chaotic trajectories by averaging a large number of such trajectories. An essential tool to accomplish these ends is our unique method [14], [15], which enables speeding up calculations by three orders of magnitude as compared to difference schemes. The approach implies combining Magnus expansion methods used to reveal a charge evolution with a stochastic difference scheme applied to calculate motion trajectories of classical sites. This computational method enabled us to model reasonably sized samples so that to calculate a velocity autocorrelation function smooth enough to be subsequently analyzed. 2. Theory 1. Discrete Holstein model in semi-classical approximation. Modelling a thermostat The charge transfer in a one-dimensional double-strand chain will be considered in the framework of the Holstein model. An electron travels along a chain whose sites are modeled by unconnected classical harmonical oscillators. The model Hamiltonian reads: ̂ = ∑ |⟩ ⟨| + ∑ |⟩ ⟨| + ∑ 1 2 1 + ∑ ̇ 2, (1) , 2 2 where are matrix elements or energies of electron transitions between the sites and , is a constant of electron coupling with displacement of the -th site , and are mass and elasticity coefficient of the oscillators, respectively. We will deal with homogeneous chains in the nearest neighbor approximation: = ,±1. If we seek a solution corresponding to Hamiltonian (1), in the form |Ψ⟩ = ∑ () |⟩ , then the relevant Schroedinger equation for a charge and the equation for the sites motion at finite temperature is written as: ℏ = (-1 + +1) + , (2a) 2 2 2 = - - || + Γ + () . (2b) To model the temperature of the surrounding medium we place the chain in a Langevin thermostat. For this purpose equation for classical sites (2b) will include a term with friction (where Γ is a friction coefficient) and a random force () with a normal distribution and the autocorrelation function: ⟨ () (′)⟩ = 2Γ , ( - ′) , where is temperature and is the Bolzmann constant. 2. Change-over to dimensionless parameters. Scaling of autocorrelation functions If we neglect the influence of quantum equation (2a) on classical equation (2b), the latter becomes a Langevin-type equation, and the solutions do not depend on the charge distribution over the chain. In paper [16], we called this variant of the system an “adiabatic approximation”. In the complete system, the quantum subsystem affects classical displacements but to a limited extent. Therefore we explicitly separate a temperature multiplier in front of them, proceeding from the expression for dispersion of the oscillator coordinate with the elasticity coefficient : ⟨2 ⟩ = . This expression follows from virial theorem and determines the scale of variation of a classical oscillator amplitude as a function of temperature: = √ . The time scale specifies the matrix element = ℏ/. Let us rewrite equations (2) with regard to the chosen scales in dimensionless form: = -1 + +1 + , (3a) 2 2 2 2 = -0 - || + + () , (3b) where 0 ⟨ () ( ′)⟩ = 22 ( - ′) , ℏ ℏ Γ 2 ℏ2 0 = √ , = , = 3, = √ . The quantum equation contains only one dimensionless parameter , which is the amplitude multiplier of classical displacements. It is mainly responsible for the influence of the classical subsystem on the quantum one, depending on temperature. This fact suggests that the electron distribution kinetics may mainly depend on for all chain sequences. In particular, the influence of the classical system becomes negligible for → 0, and the dependence on the nucleotide sequence completely disappears. Our results demonstrate that this choice of the scale is appropriate. Apart from frequency 0 and friction , responsible for oscillator motion characteristics, the classical equation involves a parameter /, which determines a reverse influence of the quantum subsystem on the classical one. Analysis of / suggests that this influence is neglectable for high temperatures. Besides, we also can neglect it in the case of nucleotide chains with a large matrix element . If we set a constant = 0 for these cases, we obtain an adiabatic approximation considered earlier in [16]. 3. Charge mean-square displacement and generalized frequency spectrum The kinetic properties of stochastic system (3) are described by the velocity autocorrelation function and the generalized frequency spectrum. These functions can be calculated from the density of probability of finding an electron at site in a chain of length 2 + 1: (, ) = () ∗ () . In the semi-classical approximation the root-mean-square displacements in dimensionless form are calculated by the formula: ⟨2⟩ = ( ) = ∑ =- 2 (, ). (4) Then the diffusion coefficient can be found as a slope of the curve () on its linear segment or as an asymptotic value of the derivative of the displacement with respect to time () = ()′, or via the integral of the second derivative of the displacement with respect to time () = ()″. We will call this coefficient a static one and denote it as (0): (0) = 1 ∞ ∫ () . 2 0 The natural generalization of the static diffusion coefficient is a frequencydependent diffusion coefficient (), which can be found as a generalized frequency spectrum of the second derivative of the displacement function with respect to time: () = ∞ 1 ∫ () cos (2) , (5) 2 0 where is dimensionless frequency which is related to frequency in hertz () as: = = 2 ℏ . ℎ Function () in semi-classical approximation is a velocity autocorrelation function. 3. Results and discussion At the initial moment a charge was inserted in the center of a homogeneous chain of length 2 + 1 with temperature distribution of the degrees of freedom of harmonical oscillators: (0) = ,0, (0) = (0, 1) , (0) = (0, 1) , (6) where , is Kronecker delta, (, ) is a normal distribution with the mean and dispersion 2. The trajectory of the system motion was calculated from initial states (6), each time with new values of and ̇. The calculations were carried out by the method, which combines Magnus expansion with a stochastic difference scheme [14], [15]. For the trajectories to correspond to charge propagation along an infinite chain, the chain length was chosen such that during the time of modeling the front of propagation of the charge probability density would not come up to its end. Then we calculated the second derivative of the root-mean-square displacement with respect to (3), (4) by averaging over samples. After that, we found the frequency spectrum by formula (5) using the fast cosine transform. We chose the parameters listed in Table 1 as in our previous work [13]. These parameters are typical for DNA. Dimensional parameters of the Holstein model for dsDNA Table 1 = 10-21 g = 0.062 eV/Å2 Γ = 6 ⋅ 10-10 g/s = 0.13 eV/Å = 0.030 eV = 0.041 eV = 0.084 eV = 0.158 eV We repeated the calculations for a set of temperatures and homogeneous dsDNA chains composed of different nucleotides. Table 2 lists dimensionless parameters. Dimensionless parameters for dsDNA as a function of a nucleotide sequence Table 2 parameter poly A poly C poly G poly T 0 0.02193 0.01605 0.007837 0.004167 0.01316 0.009631 0.004702 0.0025 0.004219 0.001654 1.925 × 10-4 2.894 × 10-5 A common property of the velocity autocorrelation functions at any temperature is that they are equal to 4 in the small area after an initial time instant. This value corresponds to the velocity of the ballistic motion at zero temperature, which is equal to 2. A comparison of autocorrelation functions for chains composed of different nucleotides suggests that for the same value of , the functions differ only in the attenuation asymptotics. The spectra demonstrate differences at low frequencies for finite values of and, accordingly, different diffusion coefficients (0) (see Figure 1, dashed lines). Here we use the scale were these differences are insignificant. In particular, they lack in the limit of low and high temperatures. Earlier in work [13], we carried out a more detailed analysis of the temperature dependencies of the static diffusion coefficient for homogeneous dsDNA chains with different nucleotide sequences and suggested the scales in which the temperature dependencies of (0) are close for chains with different sequences. Analyzing the spectra in this scale here, we see that differences in the frequency-dependent diffusion coefficient blend as the frequency is increasing, and virtually disappear at the frequency, where an absolute maximum of the spectrum is observed (see Figure 1, solid lines). 103 102 101 B 2D 100 A C A D B 10-1 C D 10-2 10-3 10-1 100 101 102 3 Figure 1. Temperature dependencies of the frequency-dependent coefficient; A, Poly A dsDNA; B, Poly C dsDNA; C, Poly G dsDNA; D, Poly T dsDNA. Dashed lines for the dependencies of 2 (0) on ; solid lines for the dependencies of (max 2 ()) on As the temperature changes, the form of autocorrelation functions and spectra change too. Their analysis enables us to identify four different regimes. Since we have shown that any differences between the chains with different nucleotide structures are insignificant in our scale, and their character is determined, we will deal only with Poly A dsDNA chains in what follows. In the limit of infinitely small temperatures, the autocorrelation function tends to attenuate exponent: () = (0) -/0 . Its generalized frequency spectrum has an analytical form similar to that of the real part of the frequency-dependent conductivity, according to Drude formula: () = (0) 0 0 1 + 422 2 . In accordance with fluctuation-dissipative theorem, the mobility () is related to () as [17]: (2) = (2) . The dependence of (0) on temperature is given in [18]. Notice also that the maximum of the frequency-dependent diffusion coefficient and, accordingly, the mobility at the relevant frequency is greater than the static diffusion coefficient (see Figure 1). It suggests that DNA can serve as a conductor of alternating current. In the range of extra-low temperatures for near-zero , the autocorrelation functions and the spectra keep monotonously attenuating (see Figure 2) until the dimensionless parameter reaches the value = 0.14 in Poly T dsDNA or = 0.3 in Poly A dsDNA. This boundary is slightly different for chains with different nucleotide sequences (see Figure 1). For the first regime, scattering at long distances exceeding one site on the average is typical. 1. 4 2. A 3.5 3 102 B 2.5 A 2D 2 D C B A C 1.5 101 1 D 0.5 0 100 -0.5 100 101 102 103 = 10-4 10-3 10-2 10-1 100 ? Figure 2. (a) Velocity autocorrelation functions and (b) generalized frequency spectrum of the charge propagation along Poly A dsDNA in the range of low and extra-low temperatures; A, = 0.1; B, = 0.2; C, = 0.4; D, = 0.8 In the range of low temperatures, until the dimensionless parameter becomes = 1, reflection at long distances arises. In this regime, the velocity autocorrelation function (see Figure 2(a)) is not monotone and has a single minimum in the negative range. The absolute value of the minimum increases with growing , while its coordinate decreases. Now the spectra demonstrate a low-frequency maximum whose frequency increases with growing , while the value decreases (see Figure 2(b)). Then, up to the value of = 3.5, a range of moderate temperatures is observed. In this temperature range a high-frequency maximum caused by oscillations of the autocorrelation function () occurs after the first minimum (see Figure 3(a)). At the same time, the absolute value of the first minimum of the autocorrelation function keeps growing. High-frequency oscillations arise while the function attenuates after the first minimum. Thus, in the range of moderate temperatures, reflection at long distances, and reflection from neighboring sites coexist. Looking at the spectrum 2 (), we see that the low-frequency maximum keeps decreasing, while its frequency decreases and approaches the frequency of the second maximum (Figure 3(b)). The position of the high-frequency peak on the spectrum is independent of temperature. The first and the second maxima of the spectrum are well seen and have nearly the same value for = 1.9 (see the inset in Figure 3(b)). 1. 4 (b) 3.5 3 A 2.5 B 2 C 100 A B C D A 1.5 D E 1 0.5 0 -0.5 2D 10-1 10-2 E 2.3 2.2 C 2.1 2 1_ -1 10-3 0.1 0.2 0.3 : 10-1 100 101 102 = 10-3 10-2 10-1 100 ? Figure 3. (a) Velocity autocorrelation functions and (b) generalized frequency spectrum of the charge propagation along Poly A dsDNA in the range moderate temperatures; A, = 1.2; B, = 1.6; C, = 1.9; D, = 2.4; E, = 3.5 Notice that the value of the parameter = 1.9 for which the high-frequency peak is approximately equal in value to the low-frequency one is the same for different types of chains, which is confirmed by the curves of the maxima shown in Figure 1. The value of the parameter for which the low-frequency maximum of the spectrum disappears while the high-frequency one emerges can reasonably be called a critical value. Dimensional critical temperatures corresponding to the critical parameter are listed in Table 3 for different types of chains. Table 3 Critical temperature for = 1.9 and position of high-frequency peak for different nucleotide sequences units poly A poly C poly G poly T THz 14.5 19.8 40.6 76.4 ∘ C -134.76 -14.68 811.80 3565.37 Regularity and homogeneity of the chain cause the high-frequency peak. It occurs at a frequency of = 1/. The dimension frequency depends only on matrix element: = 2 ℎ. The values of this peak frequency for the different nucleotide sequences are presented in Table 3. Obviously, in the case of inhomogeneous regular chains, these frequencies will be more than one, and the spectra will contain several relevant peaks. In the range of high temperatures, motion mainly fulfills as hopping between neighboring sites. As distinct from the previous regimes, the absolute value of the first minimum of the autocorrelation function decreases here with growing (see Figure 4(a)). A decrease of its coordinate is considerably retarded and tends to 1. At the same time, the initial attenuation occurs more sharply. Initially, the function sharply decreases to zero and then slowly goes to the first minimum. Oscillations of the autocorrelation function in this regime occur about the abscissa axis rather than about the low-frequency envelope. (a) 4 3.5 (b) 100 A 3 2.5 2 E D C B A B 10-1 2D C A 1.5 1 0.5 0 -0.5 D 10-2 E 1_ : 10-3 10-3 10-2 10-1 100 101 = 10-2 10-1 100 101 102 ? Figure 4. (a) Velocity autocorrelation functions and (b) generalized frequency spectrum of the charge propagation along Poly A dsDNA in the range of high and ultrahigh temperatures; A, = 4.4; B, = 10; C, = 40; D, = 160; E, = 640 The spectrum demonstrates only a high-frequency peak, which decreases with growing (see Figure 4(b)). The spectrum smoothly widens. Obviously, in the limit of high temperatures, the spectrum widens and becomes monotone. The charge propagation ceases. As is seen from Figures 2-4, for rather high , at any temperatures, except for extra-low and superhigh ones, () increases with growing temperature. This effect is known for static disorder [17]. 4. Conclusion We considered the frequency spectra of the diffusion coefficient (mobility) of an excess charge in a Holstein molecular chain for different temperatures. We revealed the character of motions contributing to the charge diffusion in the Holstein model. The presence of temperature causes a dynamical disorder even in an ideal chain. Let us call the main differences between the frequency dependencies of chains with static disorder and those with the dynamic disorder caused by temperature fluctuations. In the static case, the diffusion coefficient increases with growing frequency [17], while in the dynamic case, more complicated behavior is observed (see Figures 2-4). We demonstrated the existence of asymptotics at low and high temperatures. For finite temperatures, we identified four temperature ranges in which autocorrelation functions and their spectra have similar peculiarities/features. Accordingly, the charge diffusion inside these ranges demonstrates the same character. The character is changing in passing to a different temperature range. Low temperatures cause weak scattering without reflections, as a result of which the ballistic motion changes for the diffusion one. As the temperature increases, at first, a reflection at long distances emerges, then it is added with reflections from neighboring sites. Both these motions slow down the charge diffusion. In the limit of high temperatures, the motion degenerates into hopping between neighboring sites, and the charge diffusion ceases. We revealed the character of motions contributing to the charge diffusion in the Holstein model. All the results are presented in dimensionless values, which enabled us to identify the main parameters and recognize the different regimes. The autocorrelation functions in themselves, as well as their spectra, are also a significant result. Their calculation was made possible only with the use of special unique methods of numerical integration of differential equations. Presently, DNA is considered to be a promising material for constructing molecular wires. However, even homogeneous synthetic nucleotide chains at room temperatures demonstrate deficient mobility in the case of direct current [19], [20]. A possible solution to this problem, according to the results obtained, is the use of not direct but alternating current of rather high frequency.

About the authors

Dmitry A Tikhonov

Institute of Mathematical Problems of Biology Branch of Keldysh Institute of Applied Mathematics of RAS; Institute of Theoretical and Experimental Biophysics of RAS

Email: dmitry.tikhonov@gmail.com
Candidate of Physical and Mathematical Sciences, Senior researcher, Institute of Mathematical Problems of Biology Branch of Keldysh Institute of Applied Mathematics Russian Academy of Sciences 1 Professor Vitkevich St., Pushchino, Moscow Region 142290, Russian Federation; 3 Institutskaya St., Pushchino, Moscow Region 142290, Russian Federation

Egor V Sobolev

Institute of Mathematical Problems of Biology Branch of Keldysh Institute of Applied Mathematics of RAS; European Molecular Biology Laboratory, Hamburg Unit

Email: egor@embl-hamburg.de
Candidate of Physical and Mathematical Sciences, Postdoctoral fellow, European Molecular Biology Laboratory, Hamburg Unit 1 Professor Vitkevich St., Pushchino, Moscow Region 142290, Russian Federation; c/o DESY, Building 25A, Notkestraße 85, 22607 Hamburg, Germany

Victor D Lakhno

Institute of Mathematical Problems of Biology Branch of Keldysh Institute of Applied Mathematics of RAS

Email: lak@impb.ru
Doctor of Physical and Mathematical Sciences, Scientific Director, Institute of Mathematical Problems of Biology Branch of Keldysh Institute of Applied Mathematics Russian Academy of Sciences 1 Professor Vitkevich St., Pushchino, Moscow Region 142290, Russian Federation

References