Directed polaron propagation in linear polypeptides induced by intramolecular vibrations and external electric pulses
DDirected polaron propagation in linear polypeptides induced by intramolecularvibrations and external electric pulses
DOI: 10.1103/PhysRevE.98.012401
J. Luo * B. M. A. G. Piette † *[email protected] † [email protected] of Mathematical Sciences, Durham UniversityDurham, DH1 3LE, United Kingdom a r X i v : . [ c ond - m a t . o t h e r] J u l bstract We study the propagation of α -helix polarons in a model describing the non-adiabatic inter-action between an electron and a lattice of quantum mechanical oscillators at physiologicaltemperature. We show that when excited by a sub-picosecond electric pulse, as induced byexperimentally observed sub-picosecond charge separation, the polaron is displaced by up tohundreds of lattice sites before the electron becomes delocalised. We discuss biophysical im-plications of our results. Introduction
Directed motion of electrons in proteins is important to a wide range of biological processes,from cellular respiration to photosynthesis [1–4]. Past studies have tended to examine thetransport process in the framework of quantum tunnelling [5–8]. More recently, reports haveemerged in which electron transport in either linear or helical polypeptides is modelled as apolaron effect [9–13]. That is, the electron potential is modulated by the vibrational dynamicsof the protein molecule, and this interaction induces a deep potential well, trapping the electron.The “self-trapped” electron and the local molecular distortion form a quasi-particle compoundknown as a polaron, which is a familiar concept in solid state physics [14–17]. If the electronis to propagate in a lossless solitonic fashion, as opposed to dispersing and becoming a freeparticle, it must move in sync with the local molecular distortion. Thus, the integrity of thepolaron remains intact.It has been shown that, once a static polaron is formed on a one-dimensional molecular lattice,transport is possible after a kick in the lattice pinning mode [9], or under the influence of abiharmonic electric field with zero mean [10–12] or harmonic electric field with non-zero mean[13]. In the current study, we consider polarons formed in α -helicies by electrons interactingwith intramolecular vibrations. We consider both thermalised and non-thermalised lattices.Our aim is to find a suitable external forcing which facilitates directed electron motion alongthe lattice, and to describe the characteristics of such motion. We discover that an electricpulse with appropriate amplitude and hundred-femtosecond timespan can be used to displacea polaron by tens of lattice sites. Such electric pulses match exactly, both in amplitude andtimespan, those induced by the charge separation observed in biological complexes reportedin [18, 19]. When these pulses are repeated periodically in time, we find that the polaron canremain intact for several periods, during which it can be displaced by hundreds of sites. Asimilar excitation was described in [20] for DNA polarons but without thermal effects.In Section 2, we outline our mathematical model for the non-adiabatic interaction betweenan electron and lattice vibrations, and describe the physical parameters, and we derive aset of dynamical equations from the model. In Section 3 we present solutions representingstationary polaron states, comparing numerical and analytical results. Section 4 concernspropagating solutions, where we find suitable electric fields capable of displacing a polaronfrom its stationary state and sustaining its motion. We examine the effect of a single electricpulse as well as periodically repeated pulses, and we characterise the resulting polaron motionin terms of velocity and stability. By taking thermal effects into account, we investigate thestability of the polaron with respect to random forces due to temperature in the environment.At physiological temperature, not only is the polaron dynamics stable, but also the randomforces promote directed transport in the sense that, when thermal effects are present, thepolaron has a stronger tendency to move in one direction over the other.Throughout this study, values of physical parameters are chosen to correspond to an electroninteracting with intramolecular oscillators in a hydrogen-bonded peptide chain in an α -helix.The most common secondary protein structure, the α -helix consists of amino acid residueslinked by hydrogen bonds [21, 22], and it has various intermolecular and intramolecular vibra-tional degrees of freedom, most notably the stretching of the hydrogen bonds [23], and variousintramolecular vibrational modes such as amide-I vibrations in C=O double bonds [24]. Wehave chosen to consider the C=O double oscillator as the unit of our lattice, and this determinesthe key parameter values in our model. 1 The model and dynamical equations
We consider a linear molecular lattice with unit mass M and equilibrium spacing R . In theabsence of extraneous electrons, every unit cell independently oscillates with natural angularfrequency Ω. We model the lattice by quantum mechanical operators as opposed to classicalvariables, because we wish to assign a specific frequecy Ω corresponding to the amide-I mode,and it has been shown that the absorption band of a theoretical classical amide-I oscillator is 40times wider than the quantum mechanical one [25]. We consider amide-I vibration excitationson the lattice points, which could be induced by a passing electron for example, and describethe system using a Fr¨ohlich-Holstein Hamiltonian, (cid:98) H = (cid:98) H e + (cid:98) H l + (cid:98) H int , where (cid:98) H e = N (cid:88) n =0 J (cid:98) A † n (cid:98) A n − N − (cid:88) n =0 J (cid:16) (cid:98) A † n +1 (cid:98) A n + (cid:98) A † n (cid:98) A n +1 (cid:17) − N (cid:88) n =0 eE ( t ) R ( n − n ) (cid:98) A † n (cid:98) A n , (1a) (cid:98) H l = N (cid:88) n =0 (cid:32) (cid:98) P n M + M Ω (cid:98) Q n (cid:33) , (1b) (cid:98) H int = N (cid:88) n =0 χ (cid:98) Q n (cid:98) A † n (cid:98) A n , (1c)In eq. (1), n = 0 , . . . , N labels the lattice sites. (cid:98) H e describes a tight-binding electron, where (cid:98) A † n and (cid:98) A n are, respectively, the operators of electron creation and annihilation at the n th latticesite. J is the potential energy of a stationary electron in a transfer-free and distortion-freelattice, and − J is the electron exchange energy. − e is the electron charge, E ( t ) the amplitudeof an external electric field, and n an arbitrary lattice site at which the potential energy dueto E ( t ) is set to zero. The last term in eq. (1a) represents the modification of on-site electronenergies due to the presence of the electric field [12]. (cid:98) H l corresponds to the energy contributionfrom the lattice, where (cid:98) Q n and (cid:98) P n are, respectively, the displacement and conjugate momentumoperators for the n th oscillator. The form of the interaction Hamiltonian, (cid:98) H int , is derived fromthe assumption that on-site electron energies are modulated by the displacement field, (cid:98) Q n , andwe have retained only the linear term, involving coupling constant χ , in the Taylor expansionfor this modified energy [9]. The operators satisfy the commutation relations (cid:104) (cid:98) Q m , (cid:98) P n (cid:105) = − (cid:104) (cid:98) Q m , (cid:98) P † n (cid:105) = i (cid:126) δ mn , (2a) (cid:104) (cid:98) Q m , (cid:98) A n (cid:105) = (cid:104) (cid:98) Q m , (cid:98) A † n (cid:105) = 0 = (cid:104) (cid:98) P m , (cid:98) A n (cid:105) = (cid:104) (cid:98) P m , (cid:98) A † n (cid:105) , (2b)and the fermionic anti-commutation relation (cid:98) A m (cid:98) A † n + (cid:98) A † n (cid:98) A m = δ mn . (3)Denoting the vacuum states of (cid:98) H e and (cid:98) H l by, respectively, | e (cid:105) and | a (cid:105) , we have (cid:98) A n | e (cid:105) = 0and (cid:98) Q n | a (cid:105) = (cid:98) P n | a (cid:105) = 0. Then at time t the electronic state is a superposition of singleexcitations, | Ψ e ( t ) (cid:105) = (cid:80) Nn =0 α n ( t ) (cid:98) A † n | e (cid:105) for some complex coefficients α n . We are assumingthat the intramolecular oscillators are in a Glauber state [26, 27], | Ψ a ( t ) (cid:105) = exp( (cid:98) σ ( t )) | a (cid:105) ,where (cid:98) σ ( t ) = i (cid:126) N (cid:88) n =0 (cid:16) p n ( t ) (cid:98) Q n − q n ( t ) (cid:98) P n (cid:17) , (4)for some real coefficients p n , q n . The Glauber state is a pure state and therefore does not accountfor entanglement effects, so our model is semi-classical despite the appearance of the (cid:98) P n , (cid:98) Q n | Ψ( t ) (cid:105) = N (cid:88) n =0 α n ( t ) exp( (cid:98) σ ( t )) (cid:98) A † n | e (cid:105) | a (cid:105) , (5)with the normalisation condition, N (cid:88) n =0 | α n | = 1 . (6)Using eq. (2), and the fact that (cid:98) σ † = − (cid:98) σ , as well as the Baker-Hausdorff identity for quantumoperators, we deriveexp( (cid:98) σ † ) (cid:98) Q n exp( (cid:98) σ ) = (cid:98) Q n + q n , exp( (cid:98) σ † ) (cid:98) Q n exp( (cid:98) σ ) = (cid:98) Q n + 2 q n (cid:98) Q n + q n , (7a)exp( (cid:98) σ † ) (cid:98) P n exp( (cid:98) σ ) = (cid:98) P n + p n , exp( (cid:98) σ † ) (cid:98) P n exp( (cid:98) σ ) = (cid:98) P n + 2 p n (cid:98) P n + p n . (7b)It therefore follows that the expected value of the total energy of the system in state | Ψ (cid:105) , (cid:104) H (cid:105) := (cid:104) Ψ | (cid:98) H | Ψ (cid:105) , is (cid:104) H (cid:105) = N (cid:88) j =0 N (cid:88) k =0 α ∗ j α k (cid:68) (cid:12)(cid:12)(cid:12) (cid:98) A j exp( (cid:98) σ † ) (cid:98) H exp( (cid:98) σ ) (cid:98) A † k (cid:12)(cid:12)(cid:12) (cid:69) (8a)= N (cid:88) j =0 N (cid:88) k =0 α ∗ j α k (cid:42) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:98) A j (cid:98) H e + N (cid:88) m =0 (cid:16) p m M + M Ω q m χq m (cid:98) A † m (cid:98) A m (cid:17) (cid:98) A † k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:43) , (8b)where | (cid:105) = | e (cid:105) | a (cid:105) . To derive dynamical equations for α n and q n , we proceed as follows. For q n , we have d q n / d t = ∂ (cid:104) H (cid:105) /∂p n and d p n / d t = − ∂ (cid:104) H (cid:105) /∂q n , the combination of which gives M d q n d t = − (cid:16) M Ω q n + χ | α n | (cid:17) . (9)For α n , we deduce from eq. (8a), making use of eq. (5) and the Schr¨odinger equation (cid:98) H | Ψ (cid:105) = i (cid:126) ∂ | Ψ (cid:105) /∂t , that ∂ (cid:104) H (cid:105) ∂α ∗ n = i (cid:126) N (cid:88) k =0 (cid:42) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:98) A n exp( (cid:98) σ † ) (cid:18) d α k d t + α k d (cid:98) σ d t (cid:19) exp( (cid:98) σ ) (cid:98) A † k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:43) , (10)At this point, a standard treatment is to invoke the adiabatic approximation, that d (cid:98) σ/ d t isnegligible compared to d α/ d t [28]. Here we do not make such an assumption, therefore ∂ (cid:104) H (cid:105) ∂α ∗ n = i (cid:126) N (cid:88) k =0 d α k d t δ nk − N (cid:88) k =0 α k (cid:42) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:98) A n N (cid:88) m =0 (cid:18) d p m d t (cid:16) (cid:98) Q m + q m (cid:17) − d q m d t (cid:16) (cid:98) P m + p m (cid:17)(cid:19) (cid:98) A † k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:43) = i (cid:126) d α n d t + (2 W + I ) α n , (11)where W ( t ) = (cid:104) (cid:98) H l (cid:105) = 12 N (cid:88) m =0 (cid:18) M (cid:16) d q m d t (cid:17) + M Ω q m (cid:19) , I ( t ) = (cid:104) (cid:98) H int (cid:105) = N (cid:88) m =0 χq m | α m | (12)3re the expected energy contributions from the lattice and electron-lattice interaction, respec-tively. Meanwhile, from eq. (8b), ∂ (cid:104) H (cid:105) ∂α ∗ n = N (cid:88) k =0 α k (cid:42) (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:98) A n (cid:32) N (cid:88) m =0 J (cid:98) A † m (cid:98) A m − N − (cid:88) m =0 J (cid:16) (cid:98) A † m +1 (cid:98) A m + (cid:98) A † m (cid:98) A m +1 (cid:17) − N (cid:88) m =0 eE ( t ) R ( m − n ) (cid:98) A † m (cid:98) A m + N (cid:88) m =0 (cid:16) p m M + M Ω q m χq m (cid:98) A † m (cid:98) A m (cid:17)(cid:33) (cid:98) A † k (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) (cid:43) = α n (cid:2) J + χq n − eER ( n − n ) + W (cid:3) − J ( α n − + α n +1 ) . (13)Comparing eqs. (11) and (13), we obtain the following equation for α n . i (cid:126) d α n d t = α n (cid:2) J + χq n − eER ( n − n ) − ( W + I ) (cid:3) − J ( α n − + α n +1 ) . (14)Equations (9) and (14) are the coupled dynamical equations for our system. We have defined α − = α N +1 = 0, so that eq. (14) holds at the boundaries. We note that, had we used theadiabatic approximation, the − ( W + I ) in eq. (14) would have been replaced by + W . Next,we define a variable ψ n ( t ) with | ψ n | = | α n | by α n ( t ) = ψ n ( t ) exp (cid:20) it (cid:126) ( − J + 2 J ) (cid:21) , (15)which redefines the zero of energy measurements so that J = 0 in eq. (14), and the J termbecomes a discrete Laplacian, − J ( α n − + α n +1 − α n ). To account for the effect on the latticeof its thermal environment, we add two Langevin terms, − Γ d q n / d t + F n ( t ), to the r.h.s. ofeq. (9) [29, 30]. Here, Γ is a viscous damping coefficient, which depends on temperature of theenvironment, and F n ( t ) are Gaussian stochastic terms describing thermal fluctuations, withzero mean and satisfying the correlation relation, (cid:104) F m ( t ) , F n ( t (cid:48) ) (cid:105) = 2Γ k B T δ m,n δ ( t − t (cid:48) ), where k B is the Boltzmann constant and T is the temperature. Scaling time by Ω − and lengthby √ (cid:126) M − Ω − , we have the following non-dimensionalised dynamical equations for ψ n anddimensionless lattice displacements u n . i ˙ ψ n = (cid:0) κu n − w ( τ ) − η ( τ ) (cid:1) ψ n − ρ ( ψ n − + ψ n +1 − ψ n ) − (cid:15) ( τ )( n − n ) ψ n , (16a)¨ u n = − u n − κ | ψ n | − γ ˙ u n + f n ( τ ) , (16b)where τ = Ω t, u n = q n / √ (cid:126) M − Ω − , w = W/ ( (cid:126) Ω) , η = I/ ( (cid:126) Ω), and κ = χ √ (cid:126) M Ω , ρ = J (cid:126) Ω , (cid:15) = eER (cid:126) Ω , γ = Γ M Ω , f n = F n √ (cid:126) M Ω . (17)At the boundary we have ψ − = ψ N +1 = 0. Physically, ρ is the characteristic timescaleseparation between the electron and lattice, which is why it is known in the literature asthe adiabaticity parameter [9], and κ is a measure of the coupling strength between electronand lattice. We take the following parameter values appropriate for a single hydrogen-bondedpeptide chain in an α -helix [9, 21, 22, 24, 31–33]. In cases where no exact value is known, asuitable range is given. M = 1 . × − kg, R = 4 . . × s − , J (cid:47) χ and E ( t ) as adjustable parameters in the model, so that κ and (cid:15) ( τ ) are adjustable.Appropriate estimates of Γ can be obtained by approximating the lattice units as spheres,for lack of a better method, and using Stokes’ Law [34]. Its value decreases as temperatureis raised, but for the remainder of this study we fix the dimensionless parameter γ = 0 . T = 310K. We also have ρ (cid:47)
5, and4he thermal energy θ = k B T / ( (cid:126) Ω) = 0 .
13 enters the system via the stochastic forcing, f n ( τ ),which satisfies (cid:104) f n ( τ ) (cid:105) = 0 , (cid:104) f m ( τ ) , f n ( τ + ∆ τ ) (cid:105) = 2 γθδ m,n / ∆ τ. (18)We investigate the effect of f n in Sections 3 and 4 by comparing results of solving the systemneglecting f n and solving the system with f n taken into account. With f n ( τ ) = 0, the stationary solution to eq. (16b) is u n = − κ | ψ n | . (19)This reflects the fact that lattice excitations and electronic excitations in our model are directlyrelated, one arising from the presence of the other. Using a gauge transformation ψ n (cid:55)→ ψ n exp[ i ( w + η ) τ ], where w and η are constants in the steady state, we set w + η = 0 ineq. (16a). Putting eq. (19), as well as (cid:15) ( τ ) = 0, into eq. (16a), we obtain iρ − ˙ ψ n + ( ψ n − + ψ n +1 − ψ n ) + λ | ψ n | ψ n = 0 , (20)where λ := κ ρ = χ M Ω J (21)is known as the effective coupling parameter [35]. Since ψ n is stationary, the time-evolution of ψ n manifests only as a variation in its phase. That is, for some eigenvalue H , we have ψ n ( τ ) = φ n exp ( − iρH τ ) , (22)where φ n are time-independent. Now ρ appears only as a phase factor, it is immediately clearthat the stationary | ψ n | solutions to our system are fully characterised by λ [36]. In the limit N (cid:29)
1, eq. (20) is the spatially-discretised nonlinear Schr¨odinger equation (NLSE).That is, from the NLSE i (cid:126) ∂ t (cid:101) ψ + R J ∂ x (cid:101) ψ + λJ | (cid:101) ψ | (cid:101) ψ = 0 on a domain of size RN with N (cid:29) ∂ x (cid:101) ψ = [ (cid:101) ψ ( x − R ) + (cid:101) ψ ( x + R ) − (cid:101) ψ ( x )] /R and identifying (cid:101) ψ ( nR ) with ψ n . One can therefore obtain an approximate solution to eq. (20) by discretising thestationary solution to the NLSE, namely (cid:101) ψ ( x ) = (cid:112) λ/ (cid:2) λ ( x − x ) / (4 R ) (cid:3) exp( − iJ H t/ (cid:126) )with eigenvalue H = − λ /
16 [37–39]. Indeed, the approximate stationary solution obtainedvia the continuum limit is | ψ n | = | φ n | = ( λ/
8) sech [ λ ( n − n ) / λ , whilst the profile height is proportional to λ .We note that, as well as stationary solutions, the NLSE also admits mobile solutions in theform of travelling waves. But they are of no interest to us because eq. (20) is derived from thesystem only if the u n field - an integral part of the polaron - is stationary.5 .2 Numerical solutions with f n = 0 Without invoking the continuum limit, we must compute the stationary solutions numerically.To achieve this, we use a numerical shooting method with the φ n ∼ sech( n ) approximations asinitial guess. This gives us φ n solutions, from which we can derive all the important informationabout the stationary polaron state, namely the electron probability density | ψ n | = | φ n | , thestationary lattice configuration u n = − κ | ψ n | , as well as the binding energy of the polaron, E b ,which for a general state Ψ( t ) is given w.r.t. J and in units of (cid:126) Ω by E b = (cid:104) Ψ | (cid:98) H | Ψ (cid:105) − J (cid:126) Ω = − N − (cid:88) n =0 ρ (cid:0) ψ ∗ n +1 ψ n + ψ ∗ n ψ n +1 (cid:1) + N (cid:88) n =0 (cid:32) u n u n (cid:33) + N (cid:88) n =0 κu n | ψ n | , (23)and in the stationary state the binding energy is E = − N − (cid:88) n =0 ρ (cid:0) ψ ∗ n +1 ψ n + ψ ∗ n ψ n +1 (cid:1) − N (cid:88) n =0 κ | ψ n | . (24)We note that the gauge transformation given by eq. (15) effectively shifts the energy spectrumof the system by the constant 2 J − J , which is why we measure energy from J instead of, asis commonly done, the lowest energy in the electron band, − J .
85 90 95 100 105 110 11500.20.40.60.8 n | ψ n | ρ = 4 , κ = 3 . ρ = 4 , κ = 3 . ρ = 4 , κ = 3 . ρ = 5 , κ = 3 . ρ = 3 , κ = 3 . (a)
85 90 95 100 105 110 115−0.04−0.03−0.02−0.010 q n R κ max | ψ n | ρ = 3 ρ = 5 ρ = 3 ρ = 5 (b) E b Figure 1: Stationary solutions to eq. (16) with (cid:15) ( τ ) = f n ( τ ) = 0. (a) Some | ψ n | solutions (leftaxis), and corresponding q n /R which is lattice distortion in units of equilibrium spacing (right axis),computed using various combinations of parameters ρ and κ . (b) Dependence on ρ and κ of twokey characteristics of stationary polaron states: the maximum localisation probability, max | ψ n | (leftaxis), and the binding energy, E (right axis). Both are expressed as families of functions of κ ,parametrised by ρ = 3 . , . , . , . , . We have obtained stationary solutions to eq. (16) with (cid:15) = f n = 0 on a grid of size N = 200. InFigure 1(a) we see some | ψ n | profiles, the height of which increases with κ and decreases with ρ . This is as we expected, since the stationary solutions are fully characterised by λ = κ /ρ .It also shows q n being proportional to | ψ n | , as it should be. Figure 1(b) serves to explain thedependency of the polaron state on ρ and κ . For fixed ρ , max | ψ n | increases with κ whilst E decreases with κ . The negative sign of E signifies that energy needs to be put into the systemin order to break up the polaron; thus, a decrease in E is an indication that the electronis more strongly bound to the lattice. Since a larger κ indicates a stronger electron-latticeinteraction, we do expect that it results in a more strongly bound polaron. Meanwhile, forfixed κ , max | ψ n | decreases with ρ and E decreases with ρ . An extension of Figure 1(b) isFigure 2(a), where we see max | ψ n | and E as surfaces over the parameter space ( ρ, κ ). Bydrawing contour lines of the max | ψ n | surface, we obtain Figure 2(b). Indeed, these contourlines are the parabolae κ /ρ = constant, i.e. lines of constant λ . This means that the shape ofa | ψ n | profile depends solely on λ , again as we would expect.6 (a) κρ max | ψ n | E b ρκ λ = 2 . λ = 2 . λ = 3 . (b) Contour where max | ψ n | = 0 . | ψ n | = 0 . | ψ n | = 0 . λ = κ / ρ = constant Figure 2: Dependence of max | ψ n | and E on the parameter space ( ρ, κ ). (a) max | ψ n | (positive z -axis) and E (negative z -axis) as surfaces over the ( ρ, κ ) plane. (b) Some contour lines of the max | ψ n | surface, projected onto the ( ρ, κ ) plane. Lines of κ /ρ = constant are included for comparison. τ E b Figure 3: Under non-zero stochastic forcing at θ = 0 . E b of the polaron fluctuates around astationary mean, after a rapid initial increase. λ = 2 . We thermalise the stationary po-laron by time-evolving the systemin the presence of random forces f n ( τ ), which take values according toeq. (18) with θ = 0 .
13 (310 Kelvin).For numerical stability, we use a smalltime-step of ∆ τ = 0 .
01. Regard-less of the value of λ , the polaronalways settles in a quasi-stationarystate , where its binding energy, aftera very small initial increase, oscillatesabout a steady value, and the electronprobability density fluctuates arounda steady configuration. This is a stateof thermal equilibrium. The energy variation is always very small in magnitude, as Figure 3illustrates. The changes in | ψ n | is also small, with max | ψ n | never deviating by more than0 .
01 times its initial value. These results show that our polaron is stable against thermalfluctuations acting on the lattice at physiological temperature.
With stationary polarons as initial conditions, we impose external electric fields, represented bythe (cid:15) ( τ ) term in eq. (16), and integrate our system forward in time using the RK4 method. Thetime-step remains ∆ τ = 0 .
01. We also set in eq. (16a) n = ¯ n where the stationary | ψ n | profileis maximum. Neglecting f n ( τ ), we look for suitable (cid:15) which can facilitate polaron propagation;then in Section 4.3 we investigate thermal effects by turning on f n ( τ ). We fix λ = κ /ρ = 2 . ρ we also change κ accordingly. This is because as we saw inSection 3 the electron probability distribution of stationary polarons is parametrised only by λ . We would like to study the motion of polarons with a moderate maximum localisationprobability, and indeed for λ = 2 .
80 we have max | ψ n | = 0 . λ , κ (or equivalently ρ ) is also important to the dynamics of anon-stationary polaron, in that ( λ = 2 . , κ = 3 .
35) and ( λ = 2 . , κ = 3 .
00) produce very7ifferent results. We explain this further in Sections 4.1 to 4.3.Our numerical solutions show that several natural choices of (cid:15) ( τ ) produce negative results. Asinusoidal electric field causes the polaron simply to oscillate, regardless of the field’s amplitudeand period. Using constant (cid:15) , we find that for every initial condition there exists a thresholdamplitude (cid:15) such that no polaron displacement occurs if (cid:15) < (cid:15) , and if (cid:15) ≥ (cid:15) then theelectron delocalises within several hundred units of time. Delocalisation is the phenomenonwhere the electron “escapes” the local potential well, thus destroying the polaron as its twoconstituent parts become decoupled. This can occur when excessive energy is imparted to theelectron. Although a theoretical delocalised state is represented by a probability density with | ψ n | ∼ O (1 /N ), for practical purposes we consider delocalisation to have occurred wheneverthe weaker condition max | ψ n | < . | ψ n | havethe same order of magnitude as max | ψ n | , and this is always accompanied by a significantdecay in the polaron’s binding energy. Our understanding is that the constant electric field isprone to destroying the polaron because it raises the electron energy in a sudden and continualmanner. In an attempt to counter these issues, we have used a period of linear increase in (cid:15) ( τ ) which brings it slowly to a constant value over O (10 ) units of time. However, the resultremains that as soon as (cid:15) reaches the threshold value then the electron delocalises. This callsfor a pulse-like electric potential, which peaks at a certain amplitude before resetting to zero,in theory allowing the polaron to regain stability once the peak has passed. Consider an electric pulse of the form (cid:15) ( τ ) = (cid:40) A sin ( πτ / ∆ T ) , , if τ < ∆ T,τ ≥ ∆ T, (25)where ∆ T is the timespan of the pulse and A is the amplitude. For every ∆ T we find thatthere is some critical pulse amplitude A c with the following property. If A < A c , the pulse τ P o l a r o np o s i t i o n (a) ∆ T = 3 ,A = A c = 0 . ∆ T = 30 ,A = A c = 0 . ∆ T = 300 ,A = A c = 0 . ∆ T = 300 ,A = 0 . ∆ T = 3000 ,A = A c = 0 . τ B i nd i n g e n e r g y (b) ∆ T = 3 ,A = A c = 0 . ∆ T = 3000 ,A = A c = 0 . Figure 4: Polaron motion under the electric pulse with timespan ∆ T and amplitude A . λ = 2 . κ = 3 .
35. (a) Some polaron trajectories. (b) Evolution of polaron binding energy. causes the polaron to move away and then back to the vicinity of the initial position, beforesettling in a quasi-steady state of small oscillations about the initial position. The energy of thepolaron is raised slightly by the pulse. If A ≥ A c , the polaron moves away during the pulse butdoes not return, and instead settles in an oscillatory quasi-steady state some lattice sites awayfrom its starting position. Some examples of such trajectories are shown in Figure 4(a). A c isnegatively correlated with ∆ T , which is to be expected as a longer pulse need not have as highan amplitude as a shorter one in order to impart the same amount of energy to the electron.8 .26 0.27 0.28 0.29 0.30 0.31−1012 ∆ T = 3 ∆ T = 30 D ∆ T = 300 A ∆ T = 3000 Figure 5: Polaron displacement, D , as function of∆ T and A [cf. eq. (25)]. λ = 2 . κ = 3 . A being too small. Triangles (red) indicate zerodisplacement due to delocalisation before end ofpulse. When there is displacement, the energyof the polaron is significantly higher inthe new quasi-steady state than in theinitial steady state, as Figure 4(b) illus-trates. This means that the electron haspicked up a large amount of energy fromthe pulse. Figure 5 captures the way inwhich polaron displacement, D , varies aswe increase A beyond the critical value.We see that the relationship between D and A > A c is fairly erratic, with noclear indication of positive correlation.However, we can discern the qualitativecharacteristic that, the longer the times-pan ∆ T , the larger the polaron displace-ment. This is because the system is sub-ject to the extra energy from the pulsefor longer. We note two more features ofthe polaron dynamics. Firstly, some combinations of ∆ T and A cause displacement of thepolaron in the direction of higher electric potential, i.e. smaller n . This occurs when, for in-stance, ∆ T = 30 and A = 0 .
080 or 0 . . < A < .
087 (at our resolutionof ∆ A = 0 . A causethe electron to delocalise before the end of the pulse, due to the excessive energy input. Forexample, when ∆ T = 3000, if A ≥ .
105 then delocalisation always occurs before τ = 3000,and some smaller values of A such as 0.085 produces the same effect. If we simply repeat the propagation-inducing single pulse over time, so that (cid:15) ( τ ) is periodicand takes the form (cid:15) ( τ ) = A sin ( πτ / ∆ T ) for τ ≥
0, we obtain polaron motion which isunsustainable, in the sense that delocalisation occurs within two or three periods. Whereas forthe single pulse the polaron is permanent, i.e. it remains quasi-stable once the electric field isreset to zero (as long as delocalisation has not occurred during the pulse), under the repeatedpulse the polaron is transient, as it has a finite lifetime τ at which the electron delocalises.Even though the polaron can move by several hundred lattice sites during its lifetime, itsbinding energy increases so rapidly that this type of polaron propagation cannot be consideredan efficient transport mechanism. We understand the cause of the polaron’s short lifetime tobe as follows. We saw in Section 4.1 that as a pulse hits the polaron, it raises the polaron’senergy, making the electron less bound to the lattice. Repeated applications of the same pulsetherefore eventually decouples the electron from the lattice. Crucially, the time it takes thepolaron’s binding energy to re-settle at a quasi-steady value can be significantly longer thanthe timespan of the pulse [cf. Figure 4(b)]. This means that right at the end of the first pulsethe binding energy is much higher than it would be in its quasi-steady state, and hitting thesystem with a second pulse straight away would raise the energy even higher. Thus, to prolongthe polaron lifetime, we set (cid:15) ( τ ) to zero after each pulse for an amount of time equal to some S ∆ T , which we call the relaxation period, allowing the system the necessary time to settle in9 new quasi-steady state before another pulse hits. We write this periodic forcing as (cid:15) ( τ ) = (cid:40) A sin ( πτ / ∆ T ) , , if τ − c (cid:2) (1 + S ) ∆ T (cid:3) < ∆ T,τ − c (cid:2) (1 + S ) ∆ T (cid:3) ≥ ∆ T, (26)where c is the largest integer such that τ − c (cid:2) (1 + S ) ∆ T (cid:3) ≥
0. We find that, for ∆ T = 3000,regardless of S , delocalisation always occurs before the end of the second pulse. For this reason,we restrict ourselves to ∆ T = 3 ,
30 and 300, for which S = 10 gives a long enough relaxationperiod for our purposes. Figure 6(a) contains examples of polaron trajectories under the τ P o l a r o np o s i t i o n (a) ∆ T = 3 , A = A c = 0 . ∆ T = 30 , A = A c = 0 . ∆ T = 300 , A = A c = 0 . ∆ T = 300 , A = 0 . τ B i nd i n g e n e r g y (b) ∆ T = 3 , A = A c = 0 . ∆ T = 30 , A = A c = 0 . ∆ T = 300 , A = A c = 0 . ∆ T = 300 , A = 0 . Figure 6: Polaron motion under periodic pulses with timespan ∆ T and amplitude A , andrelaxation period 10∆ T . λ = 2 . κ = 3 .
35. (a) Some polaron trajectories. (b) Evolutionof polaron binding energy. periodic forcing, and Figure 6(b) shows the corresponding evolutions of the polaron’s bindingenergy. We see that by adding a relaxation period after each pulse we allow time for the polaronto settle into a quasi-stationary state, hence the periodic lowering of binding energy. As thepolaron stabilises, its movement stalls, hence the ladder-like trajectories featuring jumps oftens of lattice sites followed by plateaus. Compared to the single pulse of Section 4.1, periodicpulses cause much larger polaron displacements. Moreover, we saw in Section 4.1 that a pulsealways raises the polaron energy (even if it does not cause a displacement), making the polaronmore susceptible to moving under further pulses, which is why the critical pulse amplitude A c for periodic pulses is lower than the A c we saw for the single pulse. It is also why, as we see inFigure 6(a), at A = A c the polaron does not begin to move until several pulses have hit.Figure 7 is a visualisation of the way a polaron moves during a pulse and settles afterwards.The polaron’s size, represented by the breadth of the electron probability function, oscillatesduring the pulse, and after each pulse the probability density always becomes broader than itwas before. Examining Figure 8, we recognise a clear negative correlation between the polaron’sdisplacement D and the pulse amplitude A for A > A c , as opposed to the erratic relationshipbetween D and A in the case of a single pulse [cf. Figure 5]. The negative correlation can beexplained as follows. The polaron’s lifetime is negatively correlated with A , as a stronger pulseraises the polaron energy by a larger amount and its repeated application causes delocalisationmore quickly. Meanwhile, the displacement per pulse tends to increase with A , as we see in the3 rd column of subfigures in Figure 8, but this increase is small compared to the decay in polaronlifetime. As a result, total displacement over the polaron’s lifetime is a decreasing function of A , for A > A c . This has the implication that A c is not only the critical amplitude, but also in asense the optimal amplitude , and it induces the largest displacement. We note in addition that10 τ n Figure 7: Evolution of the electronprobability density, | ψ n | , under pe-riodic pulses with ∆ T = 300 , A = A c = 0 . λ = 2 . , κ = 3 . | ψ n | profile at some τ . (The broader the density, thedarker the shade.) ∆ T = 3 D τ V ∆ T = 30 A ∆ T = 300 Figure 8: Polaron displacement, D , lifetime, τ , anddisplacement per pulse, V , as functions of ∆ T and A [cf. eq. (26)]. λ = 2 . κ = 3 . S = 10. under the periodic forcing the polaron propagation is directed , meaning all combinations of ∆ T and A cause displacements in the same direction, which was not the case under single-pulseforcing. Our results show that using different combinations of ρ and κ , keeping λ = κ /ρ fixed, κ A c (a) ∆ T = 3 ∆ T = 30 ∆ T = 300 ∆ T = 3000 κ A c (b) ∆ T = 3 ∆ T = 30 ∆ T = 300 Figure 9: A c as a function of κ , parametrised by ∆ T . λ = 2 .
80. (a) Single-pulseforcing by eq. (25). (b) Periodic forcing by eq. (26). produces figures which are characteristically similar to Figure 8, showing polaron displacementand lifetime as decreasing functions of A for A > A c . However, the value of A c is dependenton more than just λ . Figure 9 exhibits the dependence of A c on the electron-lattice couplingstrength κ , with λ = 2 .
80 fixed. It shows that the more strongly coupled our system is, themore difficult it becomes to displace a polaron using pulse-like electric fields, in the sense thata larger amplitude is required to cause the onset of polaron propagation. Comparing Figure9(a) and Figure 9(b), we see that periodic forcing requires a much lower pulse amplitude toachieve polaron displacement than single-pulse forcing, particularly when ∆ T is small. Indeed,when ∆ T = 3, A c for the periodic forcing is an order of magnitude smaller than that for thesingle pulse. 11 .3 Thermal stability We study the effect of stochastic forcing from the environment by first evolving the system ofeq. (16) under a non-zero f n ( τ ) and (cid:15) ( τ ) = 0 until it reaches thermal equilibrium as we describedin Section 3.3, and then turning on (cid:15) ( τ ) as per Sections 4.1 and 4.2. The stochastic term exists ∆ T = 3 ∆ T = 30 D ∆ T = 300 A ∆ T = 3000 Figure 10: Mean displacement D as functionof ∆ T and A ; single pulse, θ = 0 . λ = 2 . κ = 3 .
35. Squares (black) indicate zero dis-placement due to A being too small. Triangles(red) indicate zero displacement due to delo-calisation before end of pulse. ∆ T = 3 D τ V ∆ T = 30 A ∆ T = 300 Figure 11: Mean values of displacement, D ,lifetime, τ , and displacement per pulse, V ,as functions of ∆ T and A ; periodic pulses, θ = 0 . λ = 2 . κ = 3 . S = 10. due to thermal energy θ := k B T / ( (cid:126) Ω) and satisfies eq. (18). We fix θ = 0 .
13 which correspondsto physiological temperature (310K). For each set of parameter values ( λ, κ, ∆ T, A ), we haverun 100 numerical simulations and taken the mean values of key scalar quantities associatedwith the polaron motion, namely its displacement and, in the case of periodic pulse-like electricfields, lifetime and displacement per pulse. Figure 10 presents the way in which the polarondisplacement D depends on the amplitude A of the single-pulse forcing, under random thermalfluctuations, and it may be compared directly to Figure 5, for which f n = 0. One of the κ A c (a) ∆ T = 3 ∆ T = 30 ∆ T = 300 ∆ T = 3000 κ A c (b) ∆ T = 3 ∆ T = 30 ∆ T = 300 Figure 12: A c as a function of κ , parametrised by ∆ T . λ = 2 .
80, thermal energy θ = 0 .
13. (a) Single-pulse forcing by eq. (25). (b) Periodic forcing by eq. (26). notable effects of the stochastic forcing is making the polaron motion more directed, as Figure10 shows few combinations of ∆ T and A causing negative displacements. Figure 11 showsresults of combining the periodic excitation of eq. (26) with stochastic forcing, and we cancompare it with Figure 8 for which f n = 0. The phenomenon that polaron displacement perpulse is an increasing function of A is more clearly seen when stochastic forces are in play.The correlation between total displacement and A for A > A c , and between polaron lifetime12nd A , both of which are negative, are also stronger under thermal fluctuations. Due to thisnegative correlation, the critical pulse amplitude A c serves also as optimal amplitude, beingthe pulse strength that induces the largest displacement. Polaron propagation is stronglydirected, in that the mean value over 100 numerical simulations of total displacement is severalhundred lattice sites in the positive direction. We also note the stabilising effect of the thermalfluctuations, which is evidenced by the smoothness of the displacement function D of A , asopposed to the jagged D -versus- A curves for f n = 0 [cf. Figure 8] which exhibit significantdips at some values of A . Crucially, the polaron’s lifetime is not at all reduced in the thermalenvironment compared to its lifetime under zero random forcing. The stochastic forcing alsohas a small effect on the critical pulse amplitudes, A c , namely that it slightly decreases A c ,meaning it is easier to displace a polaron when thermal fluctuations take place. ComparingFigure 12 to Figure 9 reveals this difference. Figure 12 also reveals that the stochastic forcinghas little bearing on the way in which A c increases with the electron-lattice coupling strength, κ .Finally, we note that A c increases with ∆ T when the forcing is periodic, whilst the correlationis negative under the single-pulse forcing. We have put forward a model for the interaction between an electron and amide-I vibrations in alinear polypeptide. We have shown that the interaction can result in stationary polarons whoseprobability density and binding energy depend on an effective coupling parameter λ = κ /ρ ,where ρ and κ represent respectively the adiabaticity and coupling strength in the electron-amide system. In particular, the maximum value of the probability density function | ψ n | isproportional to λ , and a more localised | ψ n | represents a state with lower energy. To inducethe propagation of a polaron with a moderate size from its stationary state, we have usedan external excitation in the form of a squared-sinusoidal electric pulse, after constant andpermanently sinusoidal electric fields produced negative results. The timespan of a pulse rangedbetween 3 and 3000 units, which correspond respectively to 0.01ps and 10ps, ensuring that thedynamical evolution of our polaron under a single pulse takes place within the picosecondtimescale reported in [40] as the lifetime of amide-I excitations. We have discovered that forevery pulse timespan ∆ T there exists a pulse amplitude A = A c which is critical, in the sensethat an excitation displaces the polaron if and only if A ≥ A c . When displacement occurs,the polaron typically moves along the polypeptide by a distance which is positively correlatedwith ∆ T , before settling in a quasi-stationary state, with its energy raised compared to itspre-excitation level. By repeating the electric pulse periodically in time, with a sufficientlylong relaxation period between pulses to allow the polaron to settle, we have found that thepolaron can remain intact for up to tens of pulses, as each pulse causes a displacement in thesame direction along the polypeptide. For sufficiently small ∆ T , there can be many pulseswithin the characteristic amide-I lifetime, causing a total displacement by tens of peptideunits. For A > A c , the total displacement and polaron lifetime are both decreasing functionsof A , even though the displacement per pulse increases with A . Moreover, while fixing λ wehave varied the coupling strength κ in order to investigate its effect on the polaron dynamics,and we have found that A c is positively correlated with κ . Our results show that polaronpropagation induced by pulse-like electric potentials can occur at physiological temperatures.Indeed, thermally-induced stochastic forcing on the peptide units have a stabilising effect onthe system, and it promotes directed transport towards one end of the polypeptide over theother. Over the lifetime of a polaron, it can move along the polypeptide by up to hundreds ofunits. 13ur model puts the electron-lattice system in a coherent state characterised by electron cre-ation and annihilation operators, as well as position and momentum operators for the latticepoints. We have studied an alternative model where the lattice points are described by classicalvariables, and the dynamical equations for that system are the same as eq. (16) except that ω ( τ ) and η ( τ ) do not appear. The results we obtained from the alternative model were almostindistinguishable from that which we have presented in this paper, except for the fact that thequantum correction terms ω ( τ ) and η ( τ ) appear to have the effect of lowering the critical am-plitude of electric pulses by up to 5%, when thermal fluctuations are neglected. With thermalforcing taken into account, results under the two models are essentially identical. This agreeswith findings in [41].In our model we have thermalised the system classically using a Langevin term for the latticefield. F¨orner [42] argued that one should use a quantum thermalisation of the system andshowed that when doing so the thermal stability of the Davydov soliton is smaller than whatis predicted by a classical thermalisation. While it would be interesting to do a full quantumthermalisation for our model, we believe it would make very little difference as in our casethe exchange energy J is three orders of magnitude larger than in the Davydov model andapproximately 30 times larger than the thermal energy at 300 Kelvin.The strengths of electric fields involved in our model ranged from 0.5 to 20 millivolts perangstrom. In the biological cell, if a pair of opposite charges on either side of the plasmamembrane spontaneously localised, the resulting dipole would create an electric field in thecentre of the membrane with strength ∼ / ( (cid:15) r d ), where (cid:15) r and d are respectively the relativepermittivity and width of the membrane [43, 44]. Taking (cid:15) r = 5 [45] and d = 80˚A [46, 47], thiselectric field has strength 3.6 millivolts per angstrom, which is within the range of values wehave used in our model. Moreover, observations of hundred-femtosecond charge separation inbiological complexes have been reported [18, 19], the timescale of which matches the timespanof electric pulses we have considered. It is therefore conceivable that the pulse-like electricpotentials in our model could naturally occur and induce directed electron transport alongpolypeptides across the cell membrane. References [1] G. Feher, J. P. Allen, M. Y. Okamura, and D. C. Rees.
Nature , 339:111, 1989.[2] A. E. Senior, S. Nadanaciva, and J. Weber.
Biochim. Biophys. Acta , 1553:188, 2002.[3] G. Karp.
Cell and Molecular Biology . Wiley, 5th edition, 2008.[4] M. Kaucikas, K. Maghlaoui, J. Barber, T. Renger, and J. J. van Thor.
Nat. Commun ,7:13977, 2016.[5] D. N. Beratan, J. N. Onuchic, and J. J. Hopfield.
J. Chem. Phys. , 86:4488, 1987.[6] D. N. Beratan, J. N. Betts, and J. N. Onuchic.
Science , 252:1285, 1991.[7] D. N. Beratan, J. N. Betts, and J. N. Onuchic.
J. Phys. Chem. , 96:2852, 1992.[8] M. Hasegawa, H. Nagashima, R. Minobe, T. Tachikawa, H. Mino, and Y. Kobori.
J. Phys.Chem. Lett. , 8:1179, 2017.[9] D Hennig.
Phys. Rev. E , 64:041908, 2001.1410] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.
J. Phys. Condens.Matter , 20:255242, 2008.[11] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.
J. Phys. Condens.Matter , 22:155105, 2010.[12] L. S. Brizhik, A. A. Eremko, B. M. A. G. Piette, and W. J. Zakrzewski.
Phys. Rev. E ,89:062905, 2014.[13] J. Luo and B. M. A. G. Piette.
Eur. Phys. J. B , 90:155, 2017.[14] L. D. Landau.
Phys. Z. Sowjet. , 3:664, 1933.[15] S. Pekar.
Zh. Eksp. Teor. Fiz. , 16:341, 1946.[16] H. Fr¨ohlich.
Proc. R. Soc. Lond. A , 215:291, 1952.[17] T. Holstein.
Ann. Phys. , 8:325, 1959.[18] Y. Gauduel, S. Berrod, A. Migus, N. Yamada, and A. Antonetti.
Biochemistry , 27:2509,1987.[19] D. Zhong and A. H. Zewail.
Proc. Natl. Acad. Sci. U.S.A. , 98:11867, 2001.[20] V. Lakhno and A. Korshunova.
Eur. Phys. J. B , 79:147, 2011.[21] L. Pauling, R. B. Corey, and H. R. Branson.
Proc. Natl. Acad. Sci. USA , 37:205, 1951.[22] J. D. Dunitz.
Angew. Chem. Int. Ed. , 40:4167, 2001.[23] K. G. Brown, S. C. Erfurth, E. W. Small, and W. L. Peticolas.
Proc. Natl. Acad. Sci.USA , 69:1467, 1972.[24] N. A. Nevskaya and Y. N. Chirgadze.
Biopolymers , 15:637, 1976.[25] L. Cruzeiro-Hansson and S. Takeno.
Phys. Rev. E , 56:894, 1997.[26] J. C. Eilbeck, P. S. Lomdahl, and A. C. Scott.
Phys. Rev. B , 30:4703, 1984.[27] W. C. Kerr and P. S. Lomdahl.
Phys. Rev. B , 35:3629, 1987.[28] A. Scott.
Phys. Rep. , 217:1, 1992.[29] D. S. Lemons and A. Gythiel.
Am. J. Phys. , 65:1079, 1997.[30] T. Schlick.
Molecular Modeling and Simulation . Springer, 2nd edition, 2010.[31] C. Chothia, M. Levitt, and D. Richardson.
Proc. Natl. Acad. Sci. USA , 74:4130, 1977.[32] D. J. Barlow and J. M. Thornton.
J. Mol. Biol. , 201:601, 1988.[33] A. A. Zavitsas.
J. Phys. Chem. , 91:5573, 1987.[34] K. J. Laidler and J. H. Meiser.
Physical Chemistry . Benjamin/Cummings, Menlo Park,CA, 1982.[35] A. S. Alexandrov and N. Mott.
Polarons and Bipolarons . World Scientific, Singapore,1995.[36] L. Cruzeiro-Hansson, L. C. Eilbeck, J. L. Mar´ın, and F. M. Russell.
Phys. Lett. A , 266:160,2000.[37] V. E. Zakharov and A. B. Shabat.
Sov. Phys., JETP , 34:62, 1972.1538] R. Hirota.
J. Math. Phys. , 14:805, 1973.[39] Z.-D. Li, Q.-Y. Li, X.-H. Hu, Z.-X. Zheng, and Y. Sun.
Ann. Phys. , 322:2545, 2007.[40] J. Edler and P. Hamm.
J. Chem. Phys. , 117:2415, 2002.[41] L. Cruzeiro-Hansson and V. M. Kenkre.
Phys. Lett. A , 203:362, 1995.[42] W. F¨orner.
J. Phys.: Condens. Matter , 5:823, 1993.[43] J. D. Jackson
Classical Electrodynamics . Wiley, 3rd edition, 1999.[44] L. Li, C. Li, Z. Zhang, and E. Alexov.
J. Chem. Theory Comput. , 9:2126, 2013.[45] J. C. Weaver and K. H. Schoenbach.
IEEE Transactions on Dielectrics and ElectricalInsulation , 10:715, 2003.[46] L. McCaughan and S. Krimm.
Science , 207:1481, 1980.[47] R. M. Hochmuth, C. A. Evans, H. C. Wiles, and J. T. McCown.