Dynamics of the Davydov-Scott monomer in a thermal bath: comparison of the full quantum and semiclassical approaches

Dynamics of the Davydov-Scott monomer in a thermal bath: Comparison of the full quantum and semiclassical approaches J. Cuevas,1P. A. S. Silva,2F. R. Romero,1and L. Cruzeiro2 1Grupo de Física No Lineal, Universidad de Sevilla, Departamento de Física Aplicada I, Escuela Universitaria Politécnica, C/ Virgen de África, 7, 41011 Sevilla, Spain 2CCMAR and FCT, University of Algarve, Campus de Gambelas, 8000 Faro, Portugal 共Received 21 March 2007; published 12 July 2007兲 The quantum state diffusion equation is applied to the problem of energy transfer in proteins. Lindblad operators, capable of coupling the full quantum Davydov-Scott monomer to a thermal bath, are derived. Numerical simulations with the QSD equation show that the Lindblad operators derived do recreate the exact equilibrium ensemble for the Davydov-Scott monomer. Comparison of the results obtained with the full quantum and with the semiclassical systems shows that, at biological temperatures, the latter provides a good approximation of the former. DOI: 10.1103/PhysRevE.76.011907 PACS number共s兲: 87.14.Ee, 03.65.Yz I. INTRODUCTION The possibility that protein function involves the creation and transfer of vibrational excited states was first put forward in a classical paper by McClare 关1兴. This suggestion was then taken up by Davydov who was interested in muscle contraction 关2兴. In a series of studies, Davydov proposed that these vibrational excited states are amide I vibrations and that they can be stabilized by the interaction with local vibrations of the peptide group, thus affecting, and being affected by, the hydrogen bond length between peptide groups in a protein 关3兴. The analytical studies of Davydov on continuum models of ␣ helices were extended to numerical simulations of discrete models by Scott and co-workers 共for a review of the whole field until 1992 see Ref. 关4兴兲. The general conclusion from these studies was that, at low temperatures, the amide I excitation is localized and induces a shortening of the neighboring hydrogen bonds, a state that bears strong resemblances with the polaron 关7,8兴and that has been designated in the literature as the Davydov soliton. Proteins, however, must function at biological temperatures, and to be useful, the Davydov soliton must survive at these temperatures. The Davydov-Scott model includes two types of vibration: an intrasite vibration, the amide I vibration and intersite vibrations, the vibrations of the amide sites as a whole. Many of the studies assume that the latter vibrations can be treated classically, an approximation that we shall refer to here as semiclassical. The first simulations at finite temperature were performed within the semiclassical approximation and effected a coupling of the classical part of the system to a classical bath 关9兴. The result was that the Davydov soliton dispersed in a few picoseconds at biological temperatures. This result clashed with others, obtained later by Monte Carlo simulations of the full quantum system 共QMC兲, which showed that, at thermal equilibrium, the local distortion induced by the presence of an amide I excitation increases with temperature 关10兴. This conflict was resolved by showing that, in a semiclassical system, the coupling of the classical part to a classical bath leads to a classical behavior for the quantum part as well 关11兴. Thus, while the states of a classical amide I excitation, as described by the Davydov-Scott model, are predominantly delocalized at finite temperature, the states of a quantum amide I excitation in the same conditions are predominantly localized 关12兴. A set of dynamic equations, within the semiclassical approximation, was proposed 关11兴that leads to the same equilibrium averages as the QMC simulations at biological temperatures. While at low temperatures, these equations lead to localized amide I excitations that travel in an almost coherent fashion along protein ␣ helices, as in the low temperature studies 关3,4兴, at biological temperatures, although still localized, the amide I excitations jump stochastically from site to site, changing their shape and their velocity as they travel 关11,13兴. These finite temperature dynamical equations couple the stationary Schrödinger equation for the amide I excitations to a Langevin equation for the lattice and thus are only truly valid when the amide I excitations are much faster than the lattice motions 关13兴. In a recent study, the finite temperature dynamical equations in Refs. 关11,13兴were applied to amide I energy transfer in a real protein 关14兴. As it is not possible to treat all the degrees of freedom of a protein quantum mechanically, it is important to assess the accuracy of the dynamical equations in Refs. 关11,13兴. The main purpose of this study is to assess the degree of accuracy of the semiclassical model, and of the specific set of equations proposed in Refs. 关11,13兴, with respect to the full quantum model, particularly in what concerns dynamics at finite temperature. To simulate the quantum dynamics at finite temperature we will use the Quantum State Diffusion 共QSD兲theory 关15兴. QSD has been developed to describe open quantum systems and finds applications in phenomena such as the crossing of the quantum to the classical domain, decoherence and measurement processes. Here we will apply it to the study of the influence of a thermal bath on a quantum system. Since the number of variables needed to represent quantum states increase with the power of the number of amide sites, the comparison we make here is for the simplest possible Davydov-Scott Hamiltonian, representing just one peptide group 共one site for amide I excitations兲, hydrogen bonded to a nonexcitable group. The paper is organized as follows. In Sec. II the theoretical framework to deal with finite temperature dynamics in full quantum systems is developed. In Sec. III results ob- PHYSICAL REVIEW E 76, 011907 共2007兲 1539-3755/2007/76共1兲/011907共8兲©2007 The American Physical Society011907-1 tained in the semiclassical and in the full quantum approaches are presented and the paper ends with a discussion in Sec. IV. II. FINITE TEMPERATURE DYNAMICS OF THE DAVYDOV MONOMER The Hamiltonian for the full quantum Davydov monomer is H ˆq=ប ␻ 共 a ˆ†a ˆ+1 2 兲 + ␹ ⬘ua ˆ†a ˆ+1 2 ␬ u ˆ2+p ˆ2 2M, where a ˆ†共a ˆ兲is the creation 共annihilation兲operator for an amide excitation, u ˆis the operator for the displacement from the equilibrium position of the amide site, p ˆis the operator for the momentum of the amide site as a whole, ប ␻ is the energy of an amide I excitation, ␹ ⬘is the variation of the amide I energy with changes in the hydrogen bond length, ␬ is the elasticity of amide I site, and Mis the mass of the amide site. Here we adopt the values that have been most used in the literature, namely, ប ␻ =1660 cm−1, ␹ ⬘=62 pN, ␬ =13 N/m, and M =114 a.m.u. 关4兴. In terms of the operator b ˆ=冑M⍀ 2បu ˆ +ı冑1 2Mប⍀ p ˆ, the previous Hamiltonian can be written as H ˆq=ប ␻ 冉 a ˆ†a ˆ+1 2 冊 + ␹ 共b ˆ†+b ˆ兲a ˆ†a ˆ+ប⍀b ˆ†b ˆ,共1兲 where ␹ = ␹ ⬘冑ប/共2M⍀兲and ⍀=冑 ␬ /M. The Hamiltonian 共1兲only includes an extrinsic anharmonicity, due to the interaction of the amide I excitation with the hydrogen bond that connect the carbonyl groups, CvO to the amine groups H-N. Although multiquanta amide I states have a certain degree of intrinsic anharmonicity, as has been established experimentally 关5兴and considered theoretically 关6兴, we do not include it here because we restrict ourselves to one amide I quantum states only. This is because the sources for amide I vibrations we envisage here are biological processes such as the hydrolysis of ATP, which release, at most, energy for two amide I excitations. In future work we will consider also two-quanta energy states. We are mainly interested in what happens at finite temperature. Let us first consider the thermal equilibrium situation. Because at biological temperatures the amide I energy is approximately 8kBT, thermalizing the whole system, that is, both the amide I vibrations and the lattice displacements, at these temperatures or below, leads to zero amide I quanta, on average, and an amide site that is decoupled from the lattice, which fluctuates as a harmonic oscillator. We are not interested in this trivial total thermalization. What we want to know is the dynamics of the lattice when one amide I quantum is present. The Hamiltonian 共1兲conserves the number of amide I excitations and thus the thermal average of the number of amide I excitations is specified by the initial condition that we consider. On the other hand, the Hamiltonian 共1兲 does not conserve the number of quanta in the lattice, which fluctuates, at finite temperature. To determine the thermal equilibrium averages for the full quantum system, under the conditions just described, we note that the Hamiltonian 共1兲 can be diagonalized by the Lee, Low, and Pines transformation 关16兴,R ˆ=exp关 ␹ ប⍀共b ˆ†−b ˆ兲a ˆ†a ˆ兴, a procedure that leads to the transformed Hamiltonian H ˜ : H ˜ =ប ␻ 冉 A ˆ†A ˆ+1 2 冊 +ប⍀ 冉 B ˆ†B ˆ+1 2 冊 − ␹ 2 ប⍀共A ˆ†A ˆ兲2,共2兲 where the transformed operators A ˆand B ˆare given by A ˆ=R ˆ†a ˆR ˆ=a ˆexp 冋 ␹ ប⍀共b ˆ†−b ˆ兲 册 ,B ˆ=R ˆ†b ˆR ˆ=b ˆ+ ␹ ប⍀a ˆ†a ˆ. 共3兲 It is readily verified that a ˆ†a ˆ=A ˆ†A ˆ, so that these two operators have the same eigenvectors and eigenvalues. Using the diagonalized Hamiltonian 共2兲, it is easy to show that the thermal averages over lattice fluctuations, for the full quantum system, are 具具n ˆb典典 = ␯ BE + 冉 ␹ n ប⍀ 冊 2 ,具具u ˆ典典 =− ␹ ⬘n ␬ ,共4兲 where n ˆb=b ˆ†b ˆ,nis the number of amide I excitations, ␯ BE =1/关exp共ប⍀/kBT兲−1兴is the Bose-Einstein distribution, and 具具¯典典 stands for equilibrium averages over thermal fluctuations of the lattice. The corresponding thermal equilibrium values for a semiclassical Hamiltonian, in which the lattice displacement and momentum are real, are 具具nb典典 =kBT ប⍀ + 冉 ␹ n ប⍀ 冊 2 ,具具u典典 =− ␹ ⬘n ␬ ,共5兲 where kBis Boltzmann’s constant. Comparing Eq. 共4兲with Eq. 共5兲shows that the equilibrium average number of quanta in the semiclassical lattice only differs from the full quantum value because of the first term. However, as expected, the semiclassical thermal average number of lattice quanta 共5兲 tends to the quantum expression 共4兲as temperature increases. At biological temperatures, i.e., for T=37 °C, kBT/ប⍀ ⬇4.9, while the full quantum Bose-Einstein distribution ␯ BE⬇4.4, corresponding to a relative error of 11%. Moreover, the average displacements from the equilibrium position are exactly the same in the full quantum regime and in the semiclassical regime, at all temperatures. Let us now consider the non-equilibrium situation. For the full quantum system we resort to the QSD equation 关15兴 兩d ␺ 典=− ı បH ˆ兩 ␺ 典dt +兺 j 共L ˆj−具L ˆj典兲兩 ␺ 典d ␰ j +1 2兺 j 共2具Lj †典L ˆj−L ˆj †L ˆj−具L ˆj †典具L ˆj典兲兩 ␺ 典dt,共6兲 where H ˆis given by Eq. 共1兲,兩 ␺ 典is the wave function of the full quantum system, L ˆjare the Lindblad operators 关17兴that represent the influence of the environment on the system, and 具L ˆj典=具 ␺ 兩L ˆj兩 ␺ 典.InEq.共6兲the first term describes the deterministic evolution of the system, while the two other terms describe its interaction with the thermal bath. d ␰ jare independent normalized complex stochastic variables with a zero mean and a zero covariance 关15兴. The Lindblad operators L ˆjare not unique and can be defined by the Lindblad CUEVAS et al. PHYSICAL REVIEW E 76, 011907 共2007兲 011907-2 form of the master equation for the system 关17兴 d ␳ ˆ dt =− ı ប关H ˆ, ␳ ˆ兴+兺 j 冉 L ˆj ␳ ˆL ˆj †−1 2L ˆj †L ˆj ␳ ˆ−1 2 ␳ ˆL ˆj †L ˆj 冊 ,共7兲 where ␳ ˆis the density operator of the system. For a harmonic oscillator, the thermalization process can be achieved by the action of two Lindblad operators with the following form: L ˆ1=冑 ␥ 共1+ ␯ BE兲b ˆ,L ˆ2=冑 ␥ ␯ BEb ˆ†关18兴, where ␥ is an arbitrary parameter related to the intensity of the heat bath, i.e., the larger the ␥ , the stronger the interaction with the bath and the faster the relaxation to thermal equilibrium. As explained above, the thermalization we are interested in is that in which the number of amide I quanta is conserved, while the number of lattice quanta is allowed to fluctuate. To derive Lindblad operators capable of thermalizing the Davydov-Scott Hamiltonian 共1兲in this way we resort to a heuristic reasoning. The diagonalized Hamiltonian H ˜ 共2兲is linear in the lattice operators B ˆ, so, in terms of B ˆ, the Lindblad operators are the same as for the harmonic oscillator. Using the relationship between B ˆand b ˆgiven by Eq. 共3兲, we get the following Lindblad operators for our system: L ˆ1=冑 ␥ 共1+ ␯ BE兲B ˆ=冑 ␥ 共1+ ␯ BE兲 冉 b ˆ+ ␹ ប⍀a ˆ†a ˆ 冊 ,共8兲 L ˆ2=冑 ␥ ␯ BEB ˆ†=冑 ␥ ␯ BE 冉 b ˆ†+ ␹ ប⍀a ˆ†a ˆ 冊 .共9兲 With the latter Lindblad operators the bath is constituted by phonons with the same discrete energy levels as the lattice. Having obtained the Lindblad operators that thermalize the full quantum Davydov-Scott Hamiltonian, we can derive the quantum dynamical equations at finite temperature by making use of the fact that the time dependent expression for the thermal average of a given operator A ˆ, and of its derivative, are obtained from the time dependent density matrix using the relations 具A ˆ共t兲典=Tr关 ␳ ˆ共t兲A ˆ兴and 具dA ˆ共t兲/dt典 =Tr共d ␳ ˆ/dtA ˆ兲. Using the latter expressions together with Eqs. 共7兲–共9兲, leads to the following differential equations for the thermal averages of n ˆb共t兲and b ˆ共t兲: d具n ˆb典 dt =−i ␹ n ប共具b ˆ†典−具b ˆ典兲 − ␥ 具n ˆb典− ␥ ␹ n 2ប⍀共具b ˆ†典 +具b ˆ典兲 + ␥ ␯ BE,共10兲 d具b ˆ典 dt =i ␹ n ប− 冉 i⍀+ ␥ 2 冊 具b ˆ典− ␥ ␹ n 2ប⍀ .共11兲 In the latter equations and in what follows 具¯典denotes thermal average of the quantity in the square parenthesis and 具具¯典典, which was used in Eqs. 共4兲and 共5兲, denotes its thermal equilibrium value. III. RESULTS A. Averages One problem in comparing the full quantum approach with the semiclassical approach is that the states of the quantum lattice are represented by a wave function and there can be many wave functions that, when averaged over the displacement and momentum operators, lead to the same values of the displacement and momentum. Thus, one well-defined initial value problem in the semiclassical system can be mapped into a multitude of equivalent initial value problems in the corresponding quantum system. To minimize this problem we choose initial values of zero for the lattice displacement and momentum. With this choice, the solutions to the full quantum dynamical Eqs. 共10兲,共11兲are 具n ˆb共t兲典 =具具n ˆb典典 + 冋 2 冉 ␹ n ប⍀ 冊 2 −具具n ˆb典典 册 e− ␥ t/2 −2 冉 ␹ n ប⍀ 冊 2 cos共⍀t兲e− ␥ t/2,共12兲 具u ˆ共t兲典 =− ␹ ⬘n M⍀2关1 − cos共⍀t兲e− ␥ t/2兴,共13兲 具p ˆ共t兲典 =− ␹ ⬘n ⍀sin共⍀t兲e− ␥ t/2.共14兲 It is readily verified that equations 共12兲–共14兲lead to the correct thermal equilibrium averages 共4兲for the respective variables. Thus, the Lindblad operators 共8兲and 共9兲do allow for a correct thermalization of the quantum Davydov-Scott Hamiltonian 共1兲in that they lead to the exact quantum thermal equilibrium values. In the semiclassical approximation 关11,13兴, the nonequilibrium equations for the Davydov-Scott monomer are reduced to a single stochastic equation for the classical lattice displacement u, in which the number nof amide I excitations, appears as a factor in a constant driving force Mu ¨+M⍀2u+M ␥ u ˙=F共t兲− ␹ ⬘n,共15兲 where F共t兲is a stochastic variable that represents the classical bath and its time correlation obeys the fluctuationdissipation theorem 具F共t1兲F共t2兲典=2M ␥ KBT ␦ 共t1−t2兲. The ␦ function in time means that the stochastic forces F共t兲are Markovian. In this way, the nature of the bath is the same in the semiclassical and in the quantum systems. For ␹ ⬘=0, Eq. 共15兲describes a Brownian harmonic oscillator 关19兴. For ␹ ⬘ ⫽0 we have a shifted Brownian harmonic oscillator, for which the oscillations take place not around zero, but around the thermal equilibrium average 共4兲, as in the full quantum system. Two different dynamical regimes can be distinguished: overdamped 共for ␥ ⬎2⍀兲and underdamped 共for ␥ ⬍2⍀兲. In both the semiclassical and the quantum approaches, we concentrate on the underdamped regime 共for ␥ ⬍2⍀兲, more interesting because the system undergoes intrinsic oscillations. Also, a small ␥ corresponds to a weak coupling to the thermal bath, something that minimizes the effect of the Markovian nature of the stochastic forces. DYNAMICS OF THE DAVYDOV-SCOTT MONOMER IN A…PHYSICAL REVIEW E 76, 011907 共2007兲 011907-3 Another problem that arises when comparing the dynamics, at finite temperature, predicted by the two approximations, are the different representations of the thermal bath in each approximation. Indeed, inspection of Eqs. 共15兲and 共6兲 shows that while thermalization of the semiclassical system involves one real stochastic variable, the thermalization of the corresponding quantum system involves four independent real stochastic variables. For this reason, we first compare the time evolution of thermal averages in the two approximations and only consider individual trajectories later. Initial states for the classical lattice are specified by two real values: its displacement and its momentum. With the same initial conditions as before, integration of Eq. 共15兲leads to the following solutions: 具nb共t兲典 =kBT ប⍀ 冋 1− ␥ 2 4 ␤ 2cos共 ␤ t兲e− ␥ t/2 −⍀2 ␤ 2e− ␥ t 册 + 冉 ␹ n ប⍀ 冊 2 冋 1−2 冉 cos共 ␤ t兲+ ␥ 2 ␤ sin共 ␤ t兲 冊 e− ␥ t/2 冊冋 + 冉 1+ ␥ 2 2 ␤ 2sin2共 ␤ t兲+ ␥ 2 ␤ sin共2 ␤ t兲 冊 e− ␥ t 册 ,共16兲 具u共t兲典 =− ␹ ⬘n M⍀2 冋 1− 冉 cos共 ␤ t兲+ ␥ 2 ␤ sin共 ␤ t兲 冊 e− ␥ t/2 册 , 共17兲 具p共t兲典 =− ␹ ⬘n ⍀2 冉 ␤ + ␥ 2 4 ␤ 冊 sin共 ␤ t兲e− ␥ t/2,共18兲 where ␤ ⬅冑⍀2− ␥ 2/4. We see that, for t→⬁, all variables tend to their correct thermal equilibrium values, given by Eq. 共5兲. Also, in the limit ⍀ ␥ /2, the semiclassical displacement and momentum tend to the full quantum expressions 共12兲–共14兲. Thus, when comparing the two regimes we concentrate on the variable that differs most in the two regimes, namely, 具n ˆb共t兲典. Figure 1shows the dependence of the cumulative average of 具nb共t兲典, which is defined as C共具nbc共t兲典兲=1/t兰0 t具nb共 ␶ 兲典d ␶ , for T=10 and 310 K. This figure displays this dependence for the full quantum and semiclassical system together with the result of the integration of the QSD equation, using the program written by Schack and Brun 关20兴. Both at 10 and 310 K, we see that the averages obtained with the QSD equation 共dotted line兲are approximately on top of the exact analytical quantum results 共solid line兲. On the other hand, we conclude that, at T=10 K, quantum effects of the lattice are important and the semiclassical system is a poor approximation to the full quantum system. For the more biologically relevant temperature T=310 K, the quantum effects of the lattice are much less important and the semiclassical system leads approximately to the same time evolution as the full quantum system. As explained above, a direct comparison of single trajectories at finite temperature in the two regimes is precluded by the different ways in which the coupling to the thermal bath is made in each regime. However, it is of interest to know what the single trajectories look like, especially for the full quantum system. Such single trajectories are displayed in Fig. 2for the time evolution of the lattice displacement and for the number of quanta obtained in the full quantum system, from the numerical integration of the QSD equation 共6兲 共solid line兲, and in the semiclassical system, from the integration of Eq. 共15兲共dashed line兲,atT=10 K and T=310 K. Also shown in Fig. 2is a trajectory obtained with the QSD equation 共6兲for ␹ =0, in which case the Davydov-Scott monomer reduces to a simple harmonic oscillator with frequency ⍀, uncoupled to the amide I vibrations. As expected, while the Davydov-Scott monomer oscillates around a displaced equilibrium position, the harmonic oscillator oscillates around zero. The oscillations of all three systems have the frequency ⍀of the lattice and, in spite of the differences between the coupling to the thermal baths in the two regimes, the oscillations in the semiclassical system are approximately in phase with the oscillations in the quantum 0 100 200 300 400 500 60 0 0 0.1 0.2 0.3 0.4 0.5 C( <nb ( t ) > ) t (ps) Quantum Semiclassical QSD 0 100 200 300 400 500 600 0 1 2 3 4 5 6 C( <nb ( t ) > ) t (ps) Quantum Semiclassical QSD FIG. 1. 共Color online兲Time evolution of C共具n ˆb共t兲典兲 for the full quantum system 共solid line兲and for the semiclassical system 共dashed line兲, for ␥ =0.01 ⍀and n=1. The horizontal lines represent the corresponding thermal equilibrium averages. The dotted line was obtained numerically using the QSD program written by Schack and Brun 关20兴and represents an average over 500 individual trajectories, with a basis set with 200 levels for the lattice quanta and a time step of 0.01/⍀. Top panel corresponds to T=10 K, and bottom panel to T=310 K. The QSD line is not visible as it is almost coincident with the Quantum line. CUEVAS et al. PHYSICAL REVIEW E 76, 011907 共2007兲 011907-4 system for the first four cycles. At 310 K temperature, dephasing occurs faster, as expected, but recurrences are also visible. Comparing with the oscillations of the simple harmonic oscillator it becomes apparent, particularly at low temperatures, that the Davydov-Scott monomer extracts energy from the bath at a much higher rate. Indeed, the number of lattice quanta at T=310 K for ␹ =0 共dotted line兲is generally below that for ␹ ⫽0共solid line兲and for T=10 K, the number of lattice quanta for ␹ =0 is so small that it is practically on top of the horizontal axis, as can be seen in the inset at the bottom left of Fig. 2. All figures show that, even at biological temperatures, when the noise becomes significant, the quantum beats of the lattice, induced by the mere presence of an amide I excitation and enhanced by the thermal bath, are clearly visible. B. Standard deviations A variable of great interest when comparing quantum and classical approximations is the standard deviation. Let us first consider the standard deviations ␴ uand ␴ pof u ˆand p ˆ, respectively, where ␴ A=具A2典−具A典2is calculated in the full quantum regime. To this end we need the value of 具u ˆ2典and 具p ˆ2典, which are given by 具u ˆ2典=ប 2m⍀共具b ˆ†b ˆ†典+具b ˆb ˆ典+2具n ˆb典+1兲共19兲 具p ˆ2典=−mប⍀ 2共具b ˆ†b ˆ†典+具b ˆb ˆ典−2具n ˆb典−1兲,共20兲 where 具b ˆb ˆ典obeys the following differential equation: d具b ˆb ˆ典 dt =−共2i⍀+ ␥ 兲 冉 具b ˆb ˆ典− ␹ n ប⍀具b ˆ典 冊 ,共21兲 具b ˆ典being the solution to the differential equation 共11兲with initial condition 具b共0兲典=0, given by 具b ˆ共t兲典 = ␹ n ប⍀共e−i⍀t− ␥ t/2 −1兲.共22兲 Introducing this solution in Eq. 共21兲we get 0 2 4 6 8 10 12 −10 −8 −6 −4 −2 0 2 u ( t )( pm ) t (ps) QSD Semiclassical QSD (χ=0) 0 2 4 6 8 10 12 −30 −20 −10 0 10 20 30 u(t) (pm) t (ps) QSD Semiclassical QSD (χ=0) 0 2 4 6 8 10 12 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 nb (t) t (ps) 0 5 10 0 0.005 0.01 0 2 4 6 8 10 12 0 1 2 3 4 5 6 7 nb(t) t (ps) QSD Semiclassical QSD (χ=0) FIG. 2. 共Color online兲Time evolution of the lattice displacement 共upper row兲and the number of lattice quanta 共bottom row兲, for a single trajectory, in the full quantum system 共solid line兲and in the semiclassical system 共dashed line兲, as well as for the harmonic oscillator 共dotted line兲, when T=10 K 共left兲and when T=310 K 共right兲. The inset at the bottom left is a zoom up around n ˆb=0. DYNAMICS OF THE DAVYDOV-SCOTT MONOMER IN A…PHYSICAL REVIEW E 76, 011907 共2007兲 011907-5 具b ˆb ˆ典= 冉 ␹ n ប⍀ 冊 2 共1−e−i⍀t− ␥ t/2兲+Ce−2i⍀t− ␥ t共23兲 in which Cis a complex constant. In order to determine the value of Cwe need to know the initial values of ␴ uand ␴ p共which are the same as for the harmonic oscillator兲, together with their equilibrium values. They are as follows 关21兴: ␴ u共0兲=ប 2M⍀, ␴ p共0兲=បM⍀ 2,共24兲 ␴ u共⬁兲=ប 2M⍀coth ប⍀ 2KBT, ␴ p共⬁兲=បM⍀ 2coth ប⍀ 2KBT. 共25兲 With these values it can be shown that the complex constant C= ␹ n/共ប⍀兲and that the quantum standard deviations of u ˆand p ˆbecome ␴ u共t兲=ប 2M⍀关2 ␯ BE共1−e− ␥ t兲+1兴,共26兲 ␴ p共t兲=Mប⍀ 2关2 ␯ BE共1−e− ␥ t兲+1兴.共27兲 These results, that do not depend on ␹ , are the same as for the harmonic oscillator 关21兴. On the other hand, it can be shown that in the semiclassical approximation, the standard deviations ␴ u cand ␴ p chave the following expressions: ␴ u c共t兲= ␥ KBT 2M ␤ 2 冋 − 冉 ␥ 2⍀2cos共 ␤ t兲+ ␤ ⍀2sin共 ␤ t兲 冊 exp共− ␥ t/2兲 −2 ␥ exp共− ␥ t兲+2 ␤ 2 ␥ ⍀2 册 ,共28兲 ␴ p c共t兲=M ␥ KBT 2 ␤ 2 冋 − 冉 ␥ 2cos共 ␤ t兲− ␤ sin共 ␤ t兲 冊 exp共− ␥ t/2兲 −2⍀2 ␥ exp共− ␥ t兲+2 ␤ 2 ␥ 册 .共29兲 As was found in the full quantum regime, the standard deviations in the semiclassical regime do not depend on ␹ .In the limit ⍀ ␥ /2, their expressions can be approximated by ␴ u共t兲=ប ␯ MB M⍀共1−e− ␥ t兲,共30兲 ␴ p共t兲=Mប⍀ ␯ MB共1−e− ␥ t兲.共31兲 The latter equations only differ from those obtained in the full quantum regime Eqs. 共26兲and 共27兲for two reasons: first, the quantum expressions include a term associated with the uncertainty principle that states that the standard deviations cannot both be zero, even for well defined initial conditions, and, secondly, because where in the quantum expressions we have the Bose-Einstein distribution ␯ BE in the semiclassical expressions we have the classical number ␯ MB=kBT ប⍀ . The time evolution of the cumulative standard deviation is shown in Fig. 3. As was observed in Fig. 1, the semiclassical approximation is worse at T=10 K, where the relative error in the standard deviation of the displacement is approximately 67%, than at T=310 K, where that error is virtually zero, something that lends further support for the lack of importance of quantum effects of the lattice at biological temperatures. IV. DISCUSSION As far as we know, this is the first time that the QSD equation is applied to the problem of energy transfer in proteins and that Lindblad operators, capable of providing a correct thermalization of the Davydov-Scott system, are derived. Our main purpose is to make use of this full quantum formalism to evaluate the validity, at finite temperature, of the semiclassical approximation of the Davydov-Scott model. The relatively small system considered here allowed us to derive analytical expressions to important physical 0 100 200 300 400 500 600 33 33.05 33.1 33.15 33.2 C(σu(t)) (pm2) 0 100 200 300 400 500 600 0 5 10 C(σu(t)) (pm2) t (ps) 0 100 200 300 400 500 600 0 50 100 150 200 250 300 350 C(σu(t)) (pm2) t (p s ) Quantum Semiclassical QSD FIG. 3. 共Color online兲Time evolution of C关 ␴ u共t兲兴 for the full quantum system 共solid line兲and for the semiclassical system 共dashed line兲, for ␥ =0.01 ⍀and n=1. The horizontal lines represent the corresponding thermal equilibrium averages. The dotted line was obtained numerically using the QSD program 关20兴.Top panel corresponds to T=10 K, and bottom panel to T=310 K. CUEVAS et al. PHYSICAL REVIEW E 76, 011907 共2007兲 011907-6 variables, in an exact manner. A general conclusion is that, at T=10 K, quantum effects on the lattice cannot be neglected. On the other hand, at biological temperatures, the semiclassical Davydov-Scott monomer does provide a good approximation to the corresponding full quantum system. While the comparison that was made here is not possible to make for the large system that even a small real protein is, the fact that only a single site for amide I excitations was considered constitutes nevertheless a limitation of the present study. More extended systems may behave differently. Indeed, in a previous study of an extended Davydov-Scott system with 50 amide I sites 关12兴, in which only thermal equilibrium averages were considered, it was found that, at T =0.7 K, there was a 15% relative error in the semiclassical thermal equilibrium lattice displacement correlated with the position of the excitation, which disappeared above 10 K. On the other hand, for the one amide I site studied here the thermal equilibrium value of the lattice displacement is the same in the semiclassical and in the full quantum systems 关see Eqs. 共4兲and 共5兲兴, at all temperatures 共i.e., the relative error of the semiclassical thermal equilibrium lattice displacement is zero兲. In a forthcoming publication we will consider the system with two amide I sites, where the possibility of jumps of the amide I excitation from one site to another, which is likely to be responsible for the slightly different results obtained in more extended systems, is also included. It should also be said that the ultimate aim of this work is to understand how proteins work. Following McClare 关1兴, Davydov 关2,3兴, and Scott 关4兴, the main assumption is that protein function involves the storing of energy in the form of vibrational excited states, something that was recently designated as the VES hypothesis 关14兴. In recent years, experimental evidence for the feasibility of this hypothesis has been accumulating. Indeed, while an early objection to a role of vibrational excited states in bioenergetics was their short lifetimes, thought to be in the subpicosecond range 关1兴,itis now known that, even in isolated amino acids, the lifetime of the amide I mode can last approximately 1 ps and the first measurements of the lifetime of this mode in real proteins showed that, in myoglobin, the higher energy side of the amide I band has a lifetime of 15 ps 关22兴. More recently, a low-lying band of amide I has also been found in myoglobin 关23兴. This state, which is thought to arise from self-trapping, has a longer lifetime, of 30 ps at 50 K, than the normal amide I band, which is 5 ps 关23兴. On the other hand, Hamm and co-workers monitored NH excitations in a model ␣ -helix and not only find a self-trapped state for the NH stretch vibration in this system, but also that it can only arise when the helical structure is intact 关24兴, i.e., when the CvO groups are hydrogen bonded to each other, as is assumed in the Davydov-Scott model. Furthermore, experiments in the group of Hamm 关5,25–27兴also confirm the results obtained in previous computer simulations 关11,13兴according to which at low temperature vibrational excitations are self-trapped, while at biological temperatures they are localized because of static and dynamical disorder. Finally, experiments in the group of Hochstrasser 关28兴show that amide I vibrations can jump from one helix to another, as happens in computer simulations 关14兴, and as it must happen in real proteins if the VES hypothesis is valid. Thus, in spite of some perceived difficulties 关29兴it can be said that the recent experimental evidence supports the kind of interactions first postulated by Davydov 关3兴for amide I excitations in proteins and encourages the continuation of the work in this direction 关29兴.Itis thus important to make more accurate theoretical calculations that can lead to the quantitative predictions necessary to establish whether, in fact, proteins, in cells, do make use of vibrational excited states while performing their functions. The work in this paper constitutes a step in this direction. ACKNOWLEDGMENTS J.C. and F.R.R. acknowledge financial support from the MECD Project No. FIS2004-01183. The Fundação para a Ciência e a Tecnologia, Portugal, is also acknowledged for the financial support to P.A.S.S. 关1兴C. W. F. McClare, Ann. N.Y. Acad. Sci. 227,74共1974兲. 关2兴A. S. Davydov, J. Theor. Biol. 38, 559 共1973兲. 关3兴A. S. Davydov, J. Theor. Biol. 66, 379 共1977兲; Int. J. Quantum Chem. 16,5共1979兲;Biology and Quantum Mechanics 共Pergamon, New York, 1982兲. 关4兴A. Scott, Phys. Rep. 217,1共1992兲. 关5兴J. Edler and P. Hamm, J. Chem. Phys. 117, 2415 共2002兲. 关6兴V. Pouthier, Phys. Rev. E 68, 021909 共2003兲. 关7兴T. Holstein, Ann. Phys. 共N.Y.兲8, 325 共1959兲. 关8兴T. Holstein, Ann. Phys. 共N.Y.兲8, 343 共1959兲. 关9兴P. S. Lomdahl and W. C. Kerr, Phys. Rev. Lett. 55, 1235 共1985兲. 关10兴X. Wang, D. W. Brown, and K. Lindenberg, Phys. Rev. Lett. 62, 1796 共1989兲. 关11兴L. Cruzeiro-Hansson, Europhys. Lett. 33, 655 共1996兲. 关12兴L. Cruzeiro-Hansson and V. M. Kenkre, Phys. Lett. A 203, 362 共1995兲. 关13兴L. Cruzeiro-Hansson and S. Takeno, Phys. Rev. E 56, 894 共1997兲. 关14兴L. Cruzeiro, J. Phys.: Condens. Matter 17, 7833 共2005兲. 关15兴I. Percival, Quantum State Diffusion 共Cambridge University Press, Cambridge, 1998兲. 关16兴T. D. Lee, F. E. Low, and D. Pines, Phys. Rev. 90, 297 共1953兲. 关17兴G. Lindblad, Commun. Math. Phys. 48,119共1976兲. 关18兴V. V. Dodonov, S. S. Mizrahi, and A. L. de S. Silva, J. Opt. B: Quantum Semiclassical Opt. 2, 271 共2000兲. 关19兴S. Chandrasekhar, Rev. Mod. Phys. 15,1共1943兲; G. E. Uhlenbeck and L. S. Ornstein, Phys. Rev. 36, 823 共1930兲. 关20兴R. Schack and T. A. Brun, Comput. Phys. Commun. 102, 210 共1997兲. 关21兴A. Isar, A. Sandulescu, and W. Scheid, Physica A 322, 233 共2003兲. 关22兴A. Xie, L. van Der Meer, W. Hoff, and R. H. Austin, Phys. Rev. Lett. 84, 5435 共2000兲. DYNAMICS OF THE DAVYDOV-SCOTT MONOMER IN A…PHYSICAL REVIEW E 76, 011907 共2007兲 011907-7 关23兴R. H. Austin, A. Xie, L. V. Meer, M. Shinn, and G. Neil, J. Phys.: Condens. Matter 15, S1693 共2003兲. 关24兴J. Edler, R. Pfister, V. Pouthier, C. Falvo, and P. Hamm, Phys. Rev. Lett. 93, 106405 共2004兲. 关25兴J. Edler, P. Hamm, and A. C. Scott, Phys. Rev. Lett. 88, 067403 共2002兲. 关26兴J. Edler and P. Hamm, J. Chem. Phys. 119, 2709 共2003兲. 关27兴J. Edler and P. Hamm, Phys. Rev. B 69, 214301 共2004兲. 关28兴C. Fang, A. Senes, L. Cristian, W. F. DeGrado, and R. M. Hochstrasser, Proc. Natl. Acad. Sci. U.S.A. 103, 16740 共2006兲. 关29兴P. A. Lindgard and A. M. Stoneham, J. Phys.: Condens. Matter 15,V5共2003兲. CUEVAS et al. PHYSICAL REVIEW E 76, 011907 共2007兲 011907-8