Transition between Random and Periodic Electron Currents on a DNA Chain

Abstract

By resorting to a model inspired to the standard Davydov and Holstein-Fröhlich models, in the present paper we study the motion of an electron along a chain of heavy particles modeling a sequence of nucleotides proper to a DNA fragment. Starting with a model Hamiltonian written in second quantization, we use the Time Dependent Variational Principle to work out the dynamical equations of the system. It can be found that, under the action of an external source of energy transferred to the electron, and according to the excitation site, the electron current can display either a broad frequency spectrum or a sharply peaked frequency spectrum. This sequence-dependent charge transfer phenomenology is suggestive of a potentially rich variety of electrodynamic interactions of DNA molecules under the action of electron excitation. This could imply the activation of interactions between DNA and transcription factors, or between DNA and external electromagnetic fields.

1. Introduction

A central question concerning the brain/mind activity and molecular dynamics in biological systems is whether in the warm, wet and noisy environment of living matter there are some processes for which quantum phenomena can play a relevant role. Some biological phenomena have been already ascribed to quantum effects, like in photosynthesis, bird orientation, and light response of opsins. In the search for relevant quantum effects in living matter, electrons are certainly good candidates, as is confirmed by the studies on the $\pi$-stacks of aromatic rings in proteins and in neuronal microtubules. Within this broad field of investigation, we focus on the behavior of electrons along DNA strands.

The possibility of charge transfer/transport along biomolecules, and notably along DNA, can have relevant biological consequences and has thus motivated several studies [1,2,3,4,5,6,7].

Following the quotations given in Ref. [7], electron transport on DNA can be involved in the action of DNA damage response enzymes, transcription factors, or polymerase co-factors, which are relevant processes in the cell [7,8]. For example, evidence has been given [7,9] that a DNA base excision repair enzyme enters the DNA repair process [10,11] through an electron transfer mechanism. Interestingly, electron transfer through damaged regions of DNA is markedly different with respect to electron transfer through healthy regions of DNA [10].

The aim of the present work is to investigate the spectral properties of an electron current generated along a segment of DNA. The reason for this study is two-fold, on the one side to investigate how DNA can respond to external electromagnetic excitations, on the other side with the perspective of a further investigation about the possible activation of selective attractive forces between selected sites of a DNA sequence and transcription factors.

Within a quantum mechanical framework, an electron transfer along some—essentially one dimensional—substrate is described by means of the probability current of the electron wave function $\psi \left(t\right)$ which, in general, reads as

$$\begin{array}{c}\hfill \overrightarrow{j}(\overrightarrow{x},t)=\frac{e\mathsf{\hslash }}{2m_{e}i}\left({\psi }^{\ast }\overrightarrow{\nabla }\psi -\psi \overrightarrow{\nabla }{\psi }^{\ast }\right)\end{array}$$

(1)

where e and $m_e$ are the electron charge and mass, respectively, and $\mathsf{\hslash }$ is the Planck’s constant; then, according to the D’Alambert equation ${\partial }^{\mu }{\partial }_{\mu }\overrightarrow{A}\left(t\right)={\mu }_{0}e\overrightarrow{j}(\overrightarrow{x},t)$, where ${\mu }_0$ is the vacuum permeability, this current can generate an electromagnetic field of components

$$\begin{array}{c}\hfill \overrightarrow{E}\left(t\right)=-\frac{\partial \overrightarrow{A}\left(t\right)}{\partial t}\\ \hfill \overrightarrow{B}\left(t\right)=\overrightarrow{\nabla }\times \overrightarrow{A}\left(t\right).\end{array}$$

(2)

Then, the effect of the currents flowing along two macromolecules (1 and 2) can be that of generating an intermolecular force of electrodynamic kind given by the formula

$${\overrightarrow{F}}_{12}\left(t\right)=-\frac{i_1\left(t\right)i_2\left(t\right)}{c^2}\oint \oint \frac{(d{\overrightarrow{l}}_1·d{\overrightarrow{l}}_2){\overrightarrow{x}}_{12}}{|{\overrightarrow{x}}_{12}{|}^3}$$

(3)

where the double integral is a geometric form factor, and the currents $i_{1,2}\left(t\right)$ are given by the mean values

$$i_{1,2}\left(t\right)=\frac{1}{l_{1,2}}{\int }_0^{l_{1,2}}{\overrightarrow{j}}_{1,2}(\overrightarrow{x},t)d\overrightarrow{x}$$

(4)

where $l_{1,2}$ stand for the linearized lengths of molecules 1 and 2, respectively.

According to the spectral properties of these currents, and in particular in case of the presence of co-resonances in their cross frequency spectrum ${\tilde{i}}_1^{\ast }\left(\omega \right){\tilde{i}}_2\left(\omega \right)$, the interaction force in Equation (3) could attain a strength of possible relevance in a biological context. Thus, the present work aims to contribute to an ongoing research field concerning the possible activation of intermolecular electrodynamic (resonant) forces [12,13,14].

The present work makes a first step in this direction by providing an investigation of the spectral properties of electron currents along DNA fragments.

2. Definition of the Model and Its Solution

There are several possible theoretical frameworks to model exciton propagation on a lattice, and the electron propagation on a lattice belongs to this class of phenomena. A paradigmatic model is the Haken–Strobl one [15] which describes different regimes of exciton–phonon interaction leading to the exciton damping through the scattering by lattice vibrations and thus describing the exciton motion by some kind of hopping process.

In the quantum biology literature, for example, the electronic energy transport to describe the coherent conveyance of electronic energy across chromophores of protein networks—via electrodynamic coupling of their transition electric dipole moments—is described by a tight-binding Hamiltonian for an interacting N-body system in the presence of a single excitation [16].

In the present work, in order to describe electronic motions along a DNA fragment, and from the perspective of the related electrodynamic interactions, we resort to a model partly borrowed from the standard Davydov and Holstein–Fröhlich models that have been originally introduced to account for electron–phonon interaction by explicitly also describing the dynamics of the underlying molecular lattice that, at room temperature, turns out to be essentially classical. Thus, to model the electrons moving along a given DNA sequence, the following Hamiltonian operator is assumed [17,18,19]

$$\begin{array}{c}\hfill \widehat{H}={\widehat{H}}_{el}+{\widehat{H}}_{ph}+{\widehat{H}}_{int},\end{array}$$

(5)

where the electronic and phononic parts of the Hamiltonian are given by

$${\widehat{H}}_{el}=\sum _{n=1}^N\left[E_0{\widehat{B}}_n^{†}{\widehat{B}}_n+ϵ\langle {\widehat{B}}_n^{†}{\widehat{B}}_n\rangle {\widehat{B}}_n^{†}{\widehat{B}}_n+J_n({\widehat{B}}_n^{†}{\widehat{B}}_{n+1}+{\widehat{B}}_n^{†}{\widehat{B}}_{n-1})\right],$$

(6)

$${\widehat{H}}_{ph}=\frac{1}{2}\sum _{n=1}^N\left[\frac{{\widehat{p}}_n^2}{M_n}+{\Omega }_n{({\widehat{u}}_{n+1}-{\widehat{u}}_n)}^2+\frac{1}{2}\mu {({\widehat{u}}_{n+1}-{\widehat{u}}_n)}^4\right],$$

(7)

and the electron–phonon Hamiltonian reads as

$${\widehat{H}}_{int}=\sum _{n=1}^N{\chi }_n({\widehat{u}}_{n+1}-{\widehat{u}}_n){\widehat{B}}_n^{†}{\widehat{B}}_n.$$

(8)

As shown in the following, the introduction of the site-dependent electron–phonon coupling constant ${\chi }_n$ brings about a definitely richer phenomenology with respect to the site-independent case [20]. Here, we considered only a longitudinal sequence of nucleotides where ${\widehat{B}}_n$ (${\widehat{B}}_n^{†}$) is the electronic annihilation (creation) operator relative to the lattice sites $n=1,2,...,N$. The term $E_0{\widehat{B}}_n^{†}{\widehat{B}}_n$ accounts for the initial “bare” electron energy distributed on several lattice sites according to initial shape of the electron wavefunction. The new constant $ϵ$ is the nonlinear electron–electron coupling energy due to the interaction of the moving electron along the DNA molecule with the electrons of the substrate of nucleotides and accounts for the Coulomb repulsion between the traveling electron and the charges localized on the nucleotides. The site-dependent parameter $J_n$ determines the strength of the nearest neighbor coupling energies of the electron tunneling across two neighboring nucleotides.

The vibronic part takes into account the longitudinal displacements of the nucleotides from their equilibrium positions, and each nucleotide at the n-th site along a DNA segment is described by the momentum and position operators, ${\widehat{p}}_n$ and ${\widehat{u}}_n$, by the mass $M_n$ and by the site-dependent ${\Omega }_n$, the spring parameter of two neighboring nucleotides. The parameter $\mu$ is the coupling constant of the nonlinear quartic term entailing the phonon–phonon interaction, of course absent in the harmonic approximation. The quartic term is introduced as a first (stable) term—beyond the harmonic approximation—in the power-law expansion around the minimum of typically nonlinear interparticle interaction potential (e.g., as is the case of Van der Waals potentials). The parameter ${\chi }_n$ is the site-dependent coupling parameter of the electron–lattice interaction.

Here we show how the quantum equations of motion for the Hamiltonian (5) can be derived using the time dependent variational principle (TDVP) in quantum mechanics. First, we work with the second of Davydov’s ansatz state vectors by assuming the following factorization

$$\begin{array}{ccc}\hfill \left|\psi \right(t)\rangle & =& |\Psi \left(t\right)\rangle |\Phi \left(t\right)\rangle ;|\Psi \left(t\right)\rangle =\sum _n{C}_n\left(t\right){\widehat{B}}_n^{†}{|0\rangle }_{el};\hfill \\ \hfill |\Phi (t)\rangle & =& e^{-\frac{i}{\mathsf{\hslash }}\sum [{\beta }_n\left(t\right){\widehat{p}}_n-{\pi }_n\left(t\right){\widehat{u}}_n]}{|0\rangle }_{ph}\hfill \end{array}$$

(9)

normalized to unity by the condition $\langle \psi \left(t\right)|\psi \left(t\right)\rangle ={\sum }_n{\left|C_n\left(t\right)\right|}^2=1$, where $|C_n{\left(t\right)|}^2$ is the probability for finding the electron at the n-th site. The state vector $|\Psi (t)\rangle$ describes an electron given a single quantum excitation propagating along the DNA sequence of N nucleotides and $|\Phi (t)\rangle$ describes the wave function from which the expectation values of ${\widehat{u}}_n$ and ${\widehat{p}}_n$ are obtained as

$$\begin{array}{ccc}\hfill \langle \Phi \left|{\widehat{u}}_n\right|\Phi \rangle & =& {\beta }_n\left(t\right),\hfill \\ \hfill \langle \Phi \left|{\widehat{p}}_n\right|\Phi \rangle & =& {\pi }_n\left(t\right).\hfill \end{array}$$

(10)

According to TDVP, which is equivalent to the least action principle, we introduce a new wave function $\left|\varphi \right(t)\rangle$—in terms of $\left|\psi \right(t)\rangle$ in Equation (9)—by defining a time-dependent phase factor $\left(S\right(t)\in \mathbb{R})$ so that

$$\begin{array}{c}\hfill |\varphi \left(t\right)\rangle =e^{iS\left(t\right)/\mathsf{\hslash }}|\psi \left(t\right)\rangle ,\end{array}$$

(11)

which implies the normalization $\langle \varphi \left(t\right)\left|\varphi \right(t)\rangle =1$. The wave function $\left|\varphi \right(t)\rangle$ satisfies the weaker form for the Schrödinger equation, $i\mathsf{\hslash }\langle \varphi \left(t\right)|{\partial }_t|\varphi \left(t\right)\rangle =\langle \varphi \left(t\right)|\widehat{H}|\varphi \left(t\right)\rangle$, giving

$$\begin{array}{c}\hfill -\dot{S}\left(t\right)+i\mathsf{\hslash }\langle \psi \left(t\right)|{\partial }_t|\psi \left(t\right)\rangle =\langle \psi \left(t\right)|\widehat{H}|\psi \left(t\right)\rangle .\end{array}$$

(12)

Now, the TDVP reads as

$$\begin{array}{c}\hfill \delta S\left(t\right)=0\mathrm{with}S\left(t\right)={\int }_0^{t}L\left(t^{\prime }\right)d{t}^{\prime }\end{array}$$

(13)

where $L\left(t\right)$ is the Lagrangian associated to the system

$$\begin{array}{c}\hfill L\left(t\right)=i\mathsf{\hslash }\langle \psi \left(t\right)|{\partial }_t|\psi \left(t\right)\rangle -\langle \psi \left(t\right)|\widehat{H}|\psi \left(t\right)\rangle ,\end{array}$$

(14)

hence and with Equation (9) we obtain

$$\begin{array}{c}\hfill L=\sum _n\left\{i\mathsf{\hslash }{\dot{C}}_n\left(t\right)C_n^{\ast }\left(t\right)+\frac{1}{2}\left({\pi }_n\left(t\right){\dot{\beta }}_n\left(t\right)-{\dot{\pi }}_n\left(t\right){\beta }_n\left(t\right)\right)-H(C_n,C_n^{\ast },{\beta }_n,{\pi }_n)\right\},\end{array}$$

(15)

where

$$\begin{array}{c}\hfill H(C_n,C_n^{\ast },{\beta }_n,{\pi }_n)=\langle \psi \left(t\right)|\widehat{H}|\psi \left(t\right)\rangle .\end{array}$$

(16)

By requiring the fulfilment of the stationary action condition of Equation (13), we obtain

$$\begin{array}{cc}\hfill \delta S\left(t\right)& ={\sum }_n\{i\mathsf{\hslash }\left(-{\dot{C}}_n^{\ast }\left(t\right)\delta C_n\left(t\right)+{\dot{C}}_n\left(t\right)\delta C_n^{\ast }\left(t\right)\right)+{\dot{\beta }}_n\left(t\right)\delta {\pi }_n\left(t\right)-{\dot{\pi }}_n\left(t\right)\delta {\beta }_n\left(t\right)\\ & -\left({\partial }_{C_n}H\right)\delta C_n-\left({\partial }_{C_n^{\ast }}H\right)\delta C_n^{\ast }-\left({\partial }_{{\beta }_n}H\right)\delta {\beta }_n-\left({\partial }_{{\pi }_n}H\right)\delta {\pi }_n\}=0,\end{array}$$

(17)

and arrive at the dynamical equations

$$\begin{array}{ccc}\hfill i\mathsf{\hslash }{\dot{C}}_n& =& {\partial }_{C_n^{\ast }}H\hfill \\ \hfill {\dot{\beta }}_n& =& {\partial }_{{\pi }_n}H\hfill \\ \hfill {\dot{\pi }}_n& =& -{\partial }_{{\beta }_n}H.\hfill \end{array}$$

(18)

The expectation value of the Hamiltonian is

$$\begin{array}{cc}\hfill \langle \psi \left|\widehat{H}\right|\psi \rangle & ={\sum }_n[E_0|C_n{|}^2+ϵ{\left|C_n\right|}^4+J_n(C_n^{\ast }{C}_{n+1}+C_{n+1}^{\ast }{C}_n)\\ & +\frac{1}{2}\left(\frac{1}{M_n}{\pi }_n^2+{\Omega }_n{({\beta }_{n+1}-{\beta }_n)}^2+\frac{1}{2}\mu {({\beta }_{n+1}-{\beta }_n)}^4\right)\\ & +{\chi }_n({\beta }_{n+1}-{\beta }_n){\left|C_n\right|}^2].\end{array}$$

(19)

Thus, from (18) and (19), we find the following explicit form of the equations of motion

$$\begin{array}{ccc}\hfill i\mathsf{\hslash }{\dot{C}}_n& =& \left(E_0+2ϵ{\left|C_n\right|}^2+{\chi }_n({\beta }_{n+1}-{\beta }_n)\right)C_n+J_n{C}_{n+1}+J_{n-1}{C}_{n-1},\hfill \\ \hfill M_n{\ddot{\beta }}_n& =& {\Omega }_n{\beta }_{n+1}+{\Omega }_{n-1}{\beta }_{n-1}-{\Omega }_{n-1}{\beta }_n-{\Omega }_n{\beta }_n+{\chi }_n|C_n{|}^2-{\chi }_{n-1}{\left|C_{n-1}\right|}^2\hfill \\ & +& \mu \left({({\beta }_{n+1}-{\beta }_n)}^3-{({\beta }_n-{\beta }_{n-1})}^3\right).\hfill \end{array}$$

(20)

3. Physical Parameters

In order to choose meaningful physical coupling parameters of the Hamiltonian, we borrow from Refs. [21,22,23] the values of the interaction energy between an electron and each of all the four nucleotides (reported in

Table 1. Electron–Ion interaction potential (EIIP) values for nucleotides adenine (A), thymine (T), guanine (G), and cytosine (C). From Refs. [21,22].

Nucleotide	EIIP Ry	EIIP eV	Nucleotide	EIIP Ry	EIIP eV
A	0.1260	1.7143	T	0.1335	1.8164
G	0.0806	1.0966	C	0.1340	1.8232

). We assume that the moving electron—of initial energy $E_0$—moves along the sequence of nucleotides constituting a given segment of DNA by tunneling across square potential barriers of variable height and of width $a=3.4 \mathring{A}$, the average distance between two nearest neighboring nucleotides [17]. We introduce the transmission coefficients $T_{n,n+1}$ from the probability $P(n\to n\pm 1)$ of tunneling from one potential well to the nearest one in order to set the electron hopping terms $J_n$ in (6)

$$T_{n,n+1}={\left[1+\frac{E_{n+1}^2{sinh}^2\left({\beta }_{n+1}a\right)}{4E_0(E_{n+1}-E_0)}\right]}^{-1},$$

(21)

where ${\beta }_{n+1}={[2m_e(E_{n+1}-E_0)/{\mathsf{\hslash }}^2]}^{1/2}$, $m_e$ is the mass of electron and $E_{n+1}$ are the interaction energies between the moving electron and the local nucleotide. Then we assume

$$\begin{array}{c}\hfill J_n=E_0{T}_{n,n+1}.\end{array}$$

(22)

For the Hamiltonian interaction (8) we can roughly estimate the electron–phonon coupling parameter ${\chi }_n$ as

$$\begin{array}{c}\hfill {\chi }_n=dE/dx=\frac{E_{n+1}-E_n}{a}.\end{array}$$

(23)

Numerical simulations are performed by adopting dimensionless physical parameters in the dimensionless expressions of the Hamitonian (19) and of the dynamical Equation (20). These are found by rescaling time and lengths as $t={\omega }^{-1}\tau$ and ${\beta }_n=L{b}_n$, respectively, where $L=\sqrt{\mathsf{\hslash }{\omega }^{-1}{M}_n^{-1}}$. The outcomes are

$$\begin{array}{cc}\hfill \langle \psi \left|\widehat{H}\right|\psi \rangle & =\sum _{n=1}^N[E^{\prime }|C_n{|}^2+{ϵ}^{\prime }{\left|C_n\right|}^4+J_n^{\prime }(C_n^{\ast }{C}_{n+1}+C_{n+1}^{\ast }{C}_n)\hfill \\ \hfill & +\frac{1}{2}\left({\dot{b}}_n^2+{\Omega }_n^{\prime }{(b_{n+1}-b_n)}^2+\frac{1}{2}{\mu }^{\prime }{(b_{n+1}-b_n)}^4\right)\hfill \\ \hfill & +{\chi }_n^{\prime }(b_{n+1}-b_n){\left|C_n\right|}^2],\hfill \end{array}$$

(24)

and

$$\begin{array}{ccc}\hfill i\frac{d{C}_n}{d\tau }& =& \left(E^{\prime }+2{ϵ}^{\prime }{\left|C_n\right|}^2+{\chi }_n^{\prime }(b_{n+1}-b_n)\right)C_n+J_n^{\prime }{C}_{n+1}+J_{n-1}^{\prime }{C}_{n-1},\hfill \\ \hfill \frac{d^2{b}_n}{d{\tau }^2}& =& {\Omega }_n^{\prime }{b}_{n+1}+{\Omega }_{n-1}^{\prime }{b}_{n-1}-{\Omega }_{n-1}^{\prime }{b}_n-{\Omega }_n^{\prime }{b}_n+{\chi }_n^{\prime }|C_n{|}^2-{\chi }_{n-1}^{\prime }{\left|C_{n-1}\right|}^2\hfill \\ & +& {\mu }^{\prime }\left[{(b_{n+1}-b_n)}^3-{(b_n-b_{n-1})}^3\right],\hfill \end{array}$$

(25)

where

$$\begin{array}{ccc}\hfill E^{\prime }& =& \frac{E_0}{\mathsf{\hslash }\omega };{ϵ}^{\prime }=\frac{ϵ}{\mathsf{\hslash }\omega };J_n^{\prime }=\frac{J_n}{\mathsf{\hslash }\omega };\hfill \\ \hfill {\chi }_n^{\prime }& =& \frac{{\chi }_n}{\sqrt{\mathsf{\hslash }{M}_n{\omega }^3}};{\Omega }_n^{\prime }=\frac{{\Omega }_n}{M_n{\omega }^2};{\mu }^{\prime }=\frac{\mu \mathsf{\hslash }}{M_n^2{\omega }^3}.\hfill \end{array}$$

(26)

In our simulations we resort to the known sound speed $V=a{({\Omega }_n/M_n)}^{1/2}$ on DNA; we borrow the value $V=1.69$ km/s from [24] and take it as an average constant (neglecting small local variations due to the different masses of the nucleotides). Thus, we obtain the constant dimensionless parameter ${\Omega }^{\prime }=V^2/a^2{\omega }^2$ from (26). In re-writing the dynamical equations we introduce the variables

$$\begin{array}{c}\hfill q_n=\frac{C_n+C_n^{\ast }}{\sqrt{2}},p_n=\frac{C_n-C_n^{\ast }}{i\sqrt{2}},\end{array}$$

(27)

so that Equation (25) become

$$\begin{array}{cc}\hfill {\dot{q}}_n& =\left[E^{\prime }+\frac{{ϵ}^{\prime }}{2}(q_n^2+p_n^2)+{\chi }_n^{\prime }(b_{n+1}-b_n)\right]p_n+J_n^{\prime }{p}_{n+1}+J_{n-1}^{\prime }{p}_{n-1},\hfill \end{array}$$

(28)

$$\begin{array}{cc}\hfill {\dot{p}}_n& =-\left[E^{\prime }+\frac{{ϵ}^{\prime }}{2}(q_n^2+p_n^2)+{\chi }_n^{\prime }(b_{n+1}-b_n)\right]q_n+J_n^{\prime }{q}_{n+1}+J_{n-1}^{\prime }{q}_{n-1}],\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\ddot{b}}_n& ={\Omega }^{\prime }(b_{n+1}+b_{n-1}-2b_n)+\frac{1}{2}\left({\chi }_n^{\prime }(q_n^2+p_n^2)-{\chi }_{n-1}^{\prime }(q_{n-1}^2+p_{n-1}^2)\right)\hfill \\ \hfill & +{\mu }^{\prime }\left[{(b_{n+1}-b_n)}^3-{(b_n-b_{n-1})}^3\right]\hfill \\ \hfill & ={\mathcal{B}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)].\hfill \end{array}$$

(30)

where the equation for ${\ddot{b}}_n$ can be also written as

$$\begin{array}{ccc}\hfill {\dot{b}}_n& =& {\pi }_n\hfill \\ \hfill {\dot{\pi }}_n& =& {\mathcal{B}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)].\hfill \end{array}$$

(31)

Finally, by combining a leap-frog scheme for (31) and a finite differences scheme for ${\dot{q}}_n$ and ${\dot{p}}_n$, Equations (28)–(30) are rewritten in a form suitable for numerical solution, that is

$$\begin{array}{ccc}\hfill q_n(t+\Delta t)& =& q_n\left(t\right)+\Delta t{\mathcal{Q}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)],\hfill \\ \hfill p_n(t+\Delta t)& =& p_n\left(t\right)+\Delta t{\mathcal{P}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)],\hfill \\ \hfill b_n(t+\Delta t)& =& b_n\left(t\right)+\Delta t{\pi }_n\left(t\right),\hfill \\ \hfill {\pi }_n(t+\Delta t)& =& {\pi }_n\left(t\right)+\Delta t{\mathcal{B}}_n[\mathbf{b}(t+\Delta t),\mathbf{q}(t+\Delta t),\mathbf{p}(t+\Delta t)],\hfill \end{array}$$

(32)

where ${\mathcal{Q}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)]$ and ${\mathcal{P}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)]$ are the r.h.s. of Equations (28) and (29), respectively. The integration scheme for $b_n\left(t\right)$ and ${\pi }_n\left(t\right)$ is a symplectic one, meaning that all the Poincaré invariants of the associated Hamiltonian flow are conserved, among these invariants there is energy. The simple leap-frog scheme cannot be applied to the equations for ${\dot{q}}_n\left(t\right)$ and ${\dot{p}}_n\left(t\right)$ because the r.h.s. of the equations for ${\dot{q}}_n\left(t\right)$ explicitly depend on $q_n\left(t\right)$ and $b_n\left(t\right)$; therefore, the first two equations in (32) are integrated with an Euler predictor–corrector to give

$$\begin{array}{ccc}\hfill q_n^{\left(0\right)}(t+\Delta t)& =& q_n\left(t\right)+\Delta t{\mathcal{Q}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)],\hfill \\ \hfill p_n^{\left(0\right)}(t+\Delta t)& =& p_n\left(t\right)+\Delta t{\mathcal{P}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)],\hfill \\ \hfill q_n^{\left(1\right)}(t+\Delta t)& =& q_n\left(t\right)+\frac{\Delta t}{2}\left\{{\mathcal{Q}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)]+{\mathcal{Q}}_n[\mathbf{b}\left(t\right),{\mathbf{q}}^{\left(0\right)}(t+\Delta t),{\mathbf{p}}^{\left(0\right)}(t+\Delta t)]\right\},\hfill \\ \hfill p_n^{\left(1\right)}(t+\Delta t)& =& p_n\left(t\right)+\frac{\Delta t}{2}\left\{{\mathcal{P}}_n[\mathbf{b}\left(t\right),\mathbf{q}\left(t\right),\mathbf{p}\left(t\right)]+{\mathcal{P}}_n[\mathbf{b}\left(t\right),{\mathbf{q}}^{\left(0\right)}(t+\Delta t)),{\mathbf{p}}^{\left(0\right)}(t+\Delta t)]\right\},\hfill \\ \hfill b_n(t+\Delta t)& =& b_n\left(t\right)+\Delta t{\pi }_n\left(t\right),\hfill \\ \hfill {\pi }_n(t+\Delta t)& =& {\pi }_n\left(t\right)+\Delta t{\mathcal{B}}_n[\mathbf{b}(t+\Delta t),{\mathbf{q}}^{\left(1\right)}(t+\Delta t),{\mathbf{p}}^{\left(1\right)}(t+\Delta t)],\hfill \end{array}$$

(33)

so that, thanks to fact that half of the set of the dynamical Equation (32) is integrated by means of a symplectic algorithm, and half of the equations are integrated by means of the Euler predictor–corrector, it turns out that by adopting sufficiently small integration time steps $\Delta t$ the total energy is very well conserved without any drift, just with zero-mean fluctuations around a given value fixed by the initial conditions. Regarding initial conditions, independently of the specific physical excitation mechanism, we assume an electron wavefunction described by the amplitudes $C_n(t=0)$ centered at the excitation site $n=n_0$ and distributed at time $t=0$ [17] as

$$\begin{array}{ccc}\hfill C_n(t=0)& =& \frac{1}{\sqrt{8{\sigma }_0}}\mathrm{sech}\left(\frac{n-n_0}{4{\sigma }_0}\right)\hfill \end{array}$$

(34)

where ${\sigma }_0$ specifies the width of the function. Periodic boundary conditions have been used.

Regarding the initial conditions of the phonon part of the system, we assume that the components of the DNA fragment under consideration are initially thermalized at the room temperature $T=310K$. At thermal equilibrium, the average kinetic and potential energy per degree of freedom are equal, therefore at $t=0$ the displacements and the associated velocities were initialized with random values of zero mean and according to the following prescription

$$\langle |b_n{\left(0\right)|\rangle }_n=\sqrt{\frac{k_{B}T}{\mathsf{\hslash }\omega {\Omega }^{\prime }}};\langle |{\dot{b}}_n\left(0\right){|\rangle }_n=\sqrt{\frac{k_{B}T}{\mathsf{\hslash }\omega }}.$$

(35)

expressed in dimensionless form. Periodic boundary conditions have been used also for the phonon part of the system.

4. Numerical Results

In numerical simulations, we adopted an integration time step $\Delta t=5\times {10}^{-6}$ (dimensionless units) entailing a very good energy conservation with typical fluctuations of relative amplitude $\Delta E/E\simeq {10}^{-6}$.

In what follows, we report the results obtained for two sequences of $N=66$ and $N=398$ nucleotides, respectively, for different initial electron energies $E_0$, and for initial excitation sites $n_0$ entering the initial wavefunction shape (34).

As described in the Introduction, the physical quantity of interest in what follows is the Fourier spectrum of the electron current activated on a segment of DNA. The average electron current is derived from the standard probability current $j(x_i,t)$ associated with the wave function $\left|\psi \right(t)\rangle$ in (9), which, in the discretized version along the chain of nucleotides and thus ready for its numerical computation, reads as

$$\begin{array}{c}\hfill j(x_i,t)=\frac{e\mathsf{\hslash }}{2m_{e}i}\left({\psi }^{\ast }(x_i,t)\frac{\psi (x_{i+1},t)-\psi (x_{i-1},t)}{2a}-\psi (x_i,t)\frac{{\psi }^{\ast }(x_{i+1},t)-{\psi }^{\ast }(x_{i-1},t)}{2a}\right)\end{array}$$

(36)

and hence the average current

$$i_{av}\left(t\right)=\frac{1}{l}{\int }_0^{l}j(x,t)dx=\frac{1}{Na}\sum _{i=1}^{N}j(x_i,t)a.$$

(37)

In Figure 1, Figure 2, Figure 3 and Figure 4 the outcomes are reported of the simulations performed for a DNA sequence coding for a subunit of the human haemoglobin molecule (HBB) consisting of $N=398$ nucleotides (see Appendix A). Figure 1 and Figure 2 display the behavior of the system when $E_0=0.658$ eV and for $n_0=N/2$ and $n_0=N/3$, respectively. Remarkably, the mean electron current is alternate in that it takes positive and negative values; however, the corresponding Fourier power spectra appear very different according to the initial excitation site. In fact, for $n_0=N/3$ it is well evident that the spectrum $|\tilde{i}{\left(\nu \right)|}^2$ peaks around the frequency $\nu \simeq 5$ THz, whereas for $n_0=N/2$ the spectrum appears much broader and noisy.

Then, by comparing the results obtained keeping the initial excitation around the site $n_0=N/3$ but increasing the energy to $E_0=0.79$ eV, an interesting and to some extent surprising result is found: the power spectrum of the current becomes sharper and peaks around $\nu \simeq 44$ THz, a much higher frequency indeed. This is shown in Figure 3.

On the other hand, with the initial excitation again localized around the site $n_0=N/2$ and with a lower initial activation energy, $E_0=0.46$ eV, Figure 4 shows that the electron wave function quickly spreads over the whole chain of nucleotides and the electron current power spectrum broadens with respect to the one found for $E_0=0.658$ eV. Of course, the parameter space of the system is very large and thus its systematic investigation is practically unfeasible. Nevertheless, the above reported results outline the existence of a richer phenomenology with respect to the excitation of just Davydov electro-solitons. This could be of interest in view of modeling specific processes involving electrodynamic interactions of DNA with other biomolecules or with external electromagnetic fields.

Panel (b) of the different figures shows the time evolution of random initial conditions for the displacements of the underlying sequences of masses modeling nucleotides of the DNA. According to the prescriptions of Equation (35), the random initial displacements and velocities are made at thermal equilibrium at ${310}^{\circ }$ K.

In Figure 5, Figure 6, Figure 7 and Figure 8 the results are reported as numerical simulations obtained for a shorter sequence of $N=66$ nucleotides containing the GAATTC recognition site of the restriction endonuclease enzyme EcoRI [25,26] (see Appendix B). These results confirm the richness of the phenomenology previously observed for the longer DNA sequence. The spatial distribution of the probability ${|\Psi (x,t)|}^2$ in Figure 5 where $E_0=0.2$ eV and $n_0=N/3$ seems completely frozen in time even though some low amplitude ripple exists which entail a non-vanishing electron current; the power spectrum $|\tilde{i}{\left(\nu \right)|}^2$ shows some peaks concentrated in the frequency interval 5–10 THz. By increasing the initial excitation energy of the electron above $0.6$ eV, and keeping the excitation centered at the same site $n_0=N/3$, we can see in Figure 6 a quick and complete spreading of the electron wave function ${|\Psi (x,t)|}^2$ and, correspondingly, a broad and noisy power spectrum of the electron current. A further increase in the electron excitation energy at $E_0=0.79$ eV, with the initial wavefunction centered at the site $n_0=N/2$, brings about a well evident ripple of ${|\Psi (x,t)|}^2$ around its initial peak, as can be seen in Figure 7. The associated electron current power spectrum turns out to peak drastically around 48 THz. Then, considering the initial electron energy at $E_0=0.9$ eV and displacing the initial excitation around the site $n_0=2N/3$ it is observed that the peak value of the electron wavefunction decreases much faster than in the previous case ($n_0=N/2$) and also the power spectrum of the current is very different from the previous case, displaying a broad and noisy pattern, as can be seen in Figure 8.

A comment about the possible biological significance of the energies adopted for the simulations is in order. The electron excitation values of $0.2$, $0.46$, and $0.6$ eV correspond to one, two and three DNA phosphodiester bonds, respectively, and are the same of those considered in [27] where an analysis was carried out for the same restriction enzyme DNA target sequence considered in the present work.

A complementary characterization of the electron current power spectra can be performed by means of spectral entropy. This is defined à la Shannon as follows

$$S\left(t\right)=-\sum _{m=1}^{\mathrm{M}/2}{w}_m\left(t\right)ln{w}_m\left(t\right);w_m\left(t\right)=\frac{E_m\left(t\right)}{E_T\left(t\right)}$$

(38)

the weights $w_m\left(t\right)$ are normalized by

$$E_T\left(t\right)=\sum _{m=1}^{\mathrm{M}/2}{E}_m\left(t\right);E_m\left(t\right)={\nu }_m^2{\left|\tilde{i}\left({\nu }_m\right)\right|}^2$$

(39)

where $\tilde{i}\left({\nu }_m\right)$ indicates the Fourier transform of the electron current. The frequencies are ${\nu }_m=m/T=m/(M\Delta T)$ where T is the length of time window which is Fourier analyzed, $\Delta T$ is a sampling time such that $T=M\Delta T$. As the input signals are real numbers and the Fourier spectra are computed through DFT algorithm we ignore the mirror part of the spectra thus ignoring the second half of the FFT. Then, the spectral entropy is normalized to give

$$\eta \left(t\right)=\frac{S_{max}\left(t\right)-S\left(t\right)}{S_{max}\left(t\right)}$$

(40)

so that it is $\eta =1$ when the power spectrum of the electron current is monochromatic, and $\eta =0$ for a flat spectrum such that $S\left(t\right)=S_{max}\left(t\right)$. Figure 9 shows the normalized entropy $\eta$ of the electrons versus the initial energy $E_0$. In panel (a) $\eta$ versus $E_0$ is reported for the longer DNA segment under different initial conditions: for $n_0=N/2$ the normalized entropy $\eta$ takes values approximately in the interval $0.1$–$0.15$ meaning that the corresponding power spectra are broad and noisy; for initial excitation site $n_0=N/3$ it can be found $\eta =0.45$ at $E_0=0.79$ eV and we recover what has been already displayed by Figure 3, that is, the power spectrum is narrow in frequency; intermediate values of $\eta$ can be found when the excitation site is $n_0=2N/3$ at energies $E_0=0.66,0.79,0.9$ eV.

In panel (b) $\eta$ versus $E_0$ is reported for the shorter DNA segment and synoptically shows the nontrivial appearance of some kind of, loosely speaking, “resonances” in which $\eta$ displays some bumps in its $E_0$-pattern which correspond to a significant narrowing of the electron current spectra, and thus to a less noisy and more coherent behavior of the electron current.

5. Concluding Remarks

The novelty of the model investigated in the present work concerns the introduction of site-dependent electron–phonon coupling constants. The sequence of values of these coupling constants follows the sequence of nucleotides along a given DNA segment. Moreover, the sequence of probabilities of electron jumping from one site to the next one depends on the specific sequence of nucleotides. In so doing, a richer phenomenology is found with respect to a similar model recently investigated [20]. Instead of observing the propagation of a standard Davydov electro-soliton, depending on the initial excitation energy, we observed localized periodic motions of the electrons—giving rise to a narrow frequency spectrum of the electron current—or, depending on the initial excitation site, to more extended motions associated with a broad noisy frequency spectrum. In both cases, the relevant spectral range belongs to the THz frequency domain. A qualitatively similar phenomenology was found by tackling two different DNA molecules, a subunit of the haemoglobin molecule (HBB) and an oligonucleotide with a specific recognition site of a restriction enzyme. Our findings suggest that the activation of periodic currents on specific sites could be at the origin of attractive forces between the DNA and a specific effector (transcription factor, enhancer, inhibitor, and so on). The prospective developments of the present work thus concern a new investigation/explanation of the physical grounds of the Resonant Recognition Model [21,22] by looking for co-resonances in the current frequency spectra of biochemical reaction partners. Moreover, electromagnetic signaling from electronic currents flowing along a DNA strand, and DNA response to externally applied electromagnetic fields could be further investigated in the light of the approach proposed in the present work.

Finally, further developments of the present model, aimed at describing co-resonance-activated intermolecular interactions, will have to take into account the role of hydration layers of the interacting bio-molecules. In fact, as put forward in [12], hydration layers made of spatially ordered water dipole moments can deeply affect the strength of the intermolecular electrodynamic interactions and the existence and relevance of these hydration layers is experimentally proved for DNA molecules [28,29] and for proteins [30].