Energy transfer to the phonons of a macromolecule through light pumping

In the present paper we address the problem of the energy downconversion of the light absorbed by a protein into its internal vibrational modes. We consider the case in which the light receptors are fluorophores either naturally co-expressed with the protein or artificially covalently bound to some of its amino acids. In a recent work [Phys. Rev. X 8, 031061 (2018)], it has been experimentally found that by shining a laser light on the fluorophores attached to a protein the energy fed to it can be channeled into the normal mode of lowest frequency of vibration thus making the subunits of the protein coherently oscillate. Even if the phonon condensation phenomenon has been theoretically explained, the first step - the energy transfer from electronic excitation into phonon excitation - has been left open. The present work is aimed at filling this gap.

In order to sketch the conceptual framework of the present work and what has motivated it, let us begin by quoting a recent work 1 where the activation of out-of-equilibrium collective intramolecular vibrations of a model protein has been reported. This phenomenon has been induced by light pumping, realised by shining a laser light on an aqueous solution of BSA (Bovine Serum Albumin) protein molecules each one carrying a few fluorophores covalently attached to their Lysine residues. The fluorophores were excited with a blue light at 4880Å and then they re-emitted a broadband fluorescence radiation peaked at 5190Å , thus the difference between the absorbed and re-emitted photon energies resulted in a concentration of an average energy of 0.19 eV at the fluorophores sites which thus became "hot points" on each protein. A continuous energy supply of this kind was experimentally found effective to excite the vibrational modes of the proteins and, with an energy supply rate exceeding a suitable threshold, this eventually led to a phonon condensation phenomenon into the lowest vibrational frequency. The relevance of this out-of-equilibrium collective molecular vibrations consists in the possibility of activating long-range electrodynamic interactions between bio-macromolecules 2 . The reason is that, at thermal equilibrium, a macromolecule vibrates incoherently with a broad spectrum of modes, whereas the action of an external source of energy promoting a phenomenon of phonon condensation can induce the coherent motion of the molecular subunits, so that, the resulting collective vibration can bring about a large oscillating dipole moment. Under this condition long-range and resonant (thus selective) electrodynamic forces can be activated. In turn, these electrodynamic forces could help explaining the astonishing efficiency of the impressively complex biochemical machinery at work in living cells 3 , where the different actors (proteins, DNA and RNA) find their cognate partners and targets in the right place, at the right time and in the right sequence in an overcrowded environment (the cytosol). Electrodynamic resonant/selective forces are the only possible one to act at a long distance, all the others (chemical bonds, Van der Waals and electrostatic forces) are in fact either intrinsically acting at very short distances, or are screened by the freely moving small ions in the cytosol. Actually, this is a longstanding theoretical scenario 4-6 which, for several reasons, has been discarded. However, the upgrade of Fröhlich's theoretical proposition in 1,2 and the experimental outcomes reported in 1 , represent a first crucial leap forward to ascertain whether the above mentioned hypotheses can be given experimental confirmation or refutation that can be attempted with the nowadays available technology 7,8 .
The Wu-Austin Hamiltonian from which Fröhlich rate equations can be derived by resorting to time dependent perturbation theory 9  www.nature.com/scientificreports/ where a ω i ,â ω i are the quantum creation/annihilation operators for the vibrational normal modes of a biomolecule with frequency ω i . A thermal bath at temperature T B toward which the normal modes of the biomolecule dissipate energy is represented by a collection of harmonic oscillators with characteristic frequencies j whose annihilation/creation operators are b j and b j . In order to put the biomolecule out of thermal equilibrium, the external energy pumping is modeled by another thermal bath at a temperature T S ≫ T B represented by a collection of harmonic oscillators with frequencies ′ k , the quantum annihilation/creation operators of which are ĉ ′ k and c ′ k . Then, besides linear interactions among the thermal baths modes and the biomolecule modes, mode-mode interactions among the biomolecule normal modes are considered to be mediated by the modes of the former thermal bath.
The aim of the present paper is to understand, qualitatively and quantitatively, how the model of external energy feeding of proteins -through a high temperature heat bath -can be improved to better represent the experimental conditions realised in 1 where a laser light shines on dye-labeled proteins. This means that the external energy supply is done by an electron excitation which has to be converted into vibrational energy of the chain of subunits (amino acids) composing a macromolecule (protein), a process lacking description in Hamiltonian (1). This is a very fast process with respect to the phonon condensation phenomenon, therefore it is meaningful to study it separately from the latter which has already been satisfactorily understood. In what follows, we will tackle a simplified model -with respect to that described by Eq. (1) -suitably defined to separately describe the mentioned fast process of energy transfer from the light-excited electrons of the fluorophores to the phonons of a chain of particles representing a chain of amino-acids.
As we shall see, it is found that only a fraction of the initially available electron energy is released to the phonons of a biomolecule. Even an approximate estimate of this energy transfer process is very important for a better assessment of the physical conditions that are necessary to activate the intramolecular collective vibrations.

Definition of the model
In 1 the Wu and Austin 9 quantum model leading to the original Fröhlich rate equations of 4 was reformulated (by one of us among the others) in a classical framework. However, in order to refine the description of the external energy supply, by replacing a heat bath with an electron excitation, one has to get back to a quantum description. In fact, in both cases of photo-excitation and, presumably, of ionic collisions, the excitation mechanism is supposed to be mediated by the molecular electron cloud. Therefore, the model describing the phenomenon that we want to investigate is borrowed from the standard Davydov and Holstein-Fröhlich models 10-12 to account for electron-phonon interaction. Hence, the following energy operator is assumed where the first term Ĥ el is the electron energy operator with B n and B † n the annihilation and creation operators for the electron at any site n (n = 1, 2, , . . . , N) which labels the amino acid along the protein. The term E 0B † nB n accounts for the initial "bare" electron energy distributed on several lattice sites according to initial shape of the electron wavefunction. The constant J is the nearest neighbour coupling energy of the electron tunnelling across two neighbouring amino acids, and ǫ is the energy scale of the nonlinear electron-electron coupling. In this model we have considered only a longitudinal chain of amino acids. The moving electron -yielded by the excited fluorophore-interacts on its way with almost free electrons in each amino acid, and it may just propagate along the chain of amino-acids or make a disturbance which will allow a next electron to continue on the trip. Anyway, due to electron indistinguishability, the net effect is a traveling electron along the chain of amino-acids rather than an excitonic transfer because in this latter case there is no moving mass along the chain. Then the term ǫ�B † nB n �B † nB n has been introduced to take into account non-linear effects due to the interaction between the electron in motion along the chain and the electrons of the substrate of amino acids. In particular, the term takes into account effects related to the Coulombic repulsion between the traveling electron and the charges localized on the amino acids. The averaging is intended as the expectation value of B † nB n on the dynamically evolving state of the system. The second term Ĥ ph in (2) is the phonon energy operator where p n and û n are momentum and position operators for longitudinal displacements of amino acids at site n, respectively. Furthermore, M and are average values of the mass of the amino acids of a protein and of the spring constants of two neighbouring amino acids, respectively. The quartic term is a correction stemming from the power series which gives the harmonic term at the lowest order expansion around the minimum of (1) (2) H =Ĥ el +Ĥ ph +Ĥ int , www.nature.com/scientificreports/ interparticle interaction potential (typically nonlinear, as is the case, for example, of the Van der Waals potential). This term is responsible for phonon-phonon interaction, absent in the harmonic approximation; the parameter µ sets the strength of the phonon-phonon coupling. Finally, the third term Ĥ int in (2) is the electron-phonon interaction operator where χ is the energy coupling parameter. The above defined model is considerably different with respect to the original Davydov's model because of the presence of nonlinear terms in both the electron and phonon sectors of the Hamiltonian. Moreover, the assumption-discussed in the following-of the phonon sector initially at equilibrium at room temperature makes the electron wave-function spread everywhere throughout the molecular chain. This is very far from the localized electrosoliton solutions of the original Davydov's model. And this is also coherent with what is experimentally observed 1 , thus making our model reasonable.

Derivation of the dynamical equations with TDVP
In order to derive from the model Hamiltonian (2) the corresponding dynamical equations, we make a simplifying ansatz about the state vectors by assuming the following factorization in which | � describes an electron given a single quantum excitation and supposed to be free to propagate along the chain of N amino acids composing a protein where |0� el is the vacuum state of the Amide-I oscillators, and We then set where β n (t) and π n (t) are the average values of the longitudinal displacement and momentum of an amino acid, respectively.
To derive dynamical equation we now resort to the time-dependent variational principle (TDVP) in quantum mechanics. TDVP is a formulation of the time-dependent Schrödinger equation through variation of an action functional 13,14 . The Schrödinger equation is obtained by requiring that the action functional be stationary under free variation of the time-dependent state. According to this principle, we define a new wave function |φ� in terms of |ψ� in Eq. (6) as where S(t) is a time-dependent phase factor (S(t) ∈ R) , which will be determined in a self-consistent manner and the normalization condition is �φ|φ� = 1 . The wave function |φ� satisfies the Schrödinger equation which according to Eq. (10) becomes

Integrating, we obtain
We can now derive the equations of motion by requiring that the action with the Lagrangian to be stationary From Eqs. (6), (7), and (8) we write and then arrive at It is worth noting that the TDVP, being a variational approach, applies generically to any quantum system and its effectiveness depends on a reasonable choice of the initial ansatz for the state vector. Moreover, the remarkable fact is that the dynamical equations worked out by means of the TDVP are formally classical but give the time evolution of actual quantum expectation values. Now, by inspection of the equations above, the electron/phonon energy transfer mechanism appears to be due the spatial inhomogeneities of the electron wavefunction, in fact the term χ(|C n | 2 − |C n−1 | 2 ) enters the equation for β n and, reciprocally, the local displacements χ(β n+1 − β n ) enter the equations for the Ċ n . In the original Davydov's model, this feedback mechanism produces electrosolitons which is not the case for our model where the phonons are initially thermalized at room temperature and nonlinear coupling terms enter both the electron and phonon sectors of the Hamiltonian.

Definition of the physical parameters for numerical simulations
Let us see how to make a physically reasonable choice of the coupling parameters entering the Hamiltonian. We borrow from 15-17 the estimates of the interaction energy between an electron and each of all the 20 amino acids (reported in Table 1). The average value of these interaction energies is � E� = 0.74 eV with a dispersion σ E = 0.47 eV. As a first rough picture of an electron tunnelling across the sequence of amino acids constituting a protein we can consider the electron of energy E 0 moving in a periodic sequence of square potential barriers of height V 0 = 0.74 eV and of width a = 4.5Å , the average distance between two nearest neighboring amino acids 10 . We can then weigh the electron displacement operators between neighbouring sites with the probability P(n → n ± 1) of tunnelling from one potential well to the nearest ones. This is achieved by computing the transmission coefficient Moreover, the coefficient of the electron displacement term in the Hamiltonian has to be a characteristic energy scale of the process, thus a natural choice is to set J ∝ � E�T , then, www.nature.com/scientificreports/ assuming that an electron is initially excited at any given point of the chain of amino acids and that it has the same probability of moving to the left or to the right, we add a factor 1/2 so that finally we have J = 1 2 � E�T . Now, assuming E 0 = 0.19 eV as initial value of the electron energy, we find J = 0.0585 eV, whereas assuming that only a fraction δ ∈ [0, 1] of the maximum available energy is kept by the electron, for example for δ = 0.5 , we find J = 0.031 eV. For what concerns the electron-phonon coupling constant χ , we make a rough estimate of its value as In what follows, in dimensionless units, we have χ ′ = 0.81 , and J ′ = 5 with δ = 0.5 , while J ′ = 9 with δ = 1.   (29) and (30) and the above system have been numerically integrated by combining a finite differences scheme and a leap-frog scheme as follows  18 . Therefore energy is well conserved without any drift, just zero-mean fluctuations around a given energy value fixed by the initial conditions. By using sufficiently small time steps t any desired precision of energy conservation can be attained. About the initial conditions, we aim at simulating a physical situation where each photon absorbed by a fluorophore attached to a protein releases -in the average -0.19 eV of energy to the surrounding electron cloud. This energy is the difference between the energies of the absorbed photon of 4880Å and that of the re-emitted one as fluorescent radiation of 5150Å . We assume, as already stated above, that the effect of a single photon excitation is to make one electron move across the protein by tunnelling through a sequence of potential barriers.
In the experiments to which we are referring 1 each protein is labelled with 5-6 fluorochromes, and a laser light is continuously shined on the labelled proteins, therefore what we are after is modelling an elementary process and assuming, in a first approximation, a property of additivity of the same elementary process. In other words, if more than one electron is activated we assume that the resulting physical effect is the sum of a single electron effect. As a consequence, the electron initial condition is assumed to be described by a wavefunction C n (t = 0) centered at the site n = n 0 at time t = 0 10 : Then, coming to the initial conditions of the phonon component of the system, we assume a thermalized macromolecule at room temperature, that is at T = 310K . At equilibrium, the energy equipartition theorem for the Hamiltonian (4) reads where k B is the Boltzmann constant. This assumption means that just before photo-excitation the molecular vibrations are at thermal equilibrium. As a matter of fact, the phononic sector of our system is described by the Hamiltonian of the celebrated Fermi-Pasta-Ulam β-model for which energy equipartition among the normal modes is always attained 18 thus supporting this assumption. At thermal equilibrium, energy is equally shared among all the degrees of freedom and, in particular, between kinetic and potential energies, therefore at t = 0 the velocities and the displacements have been initialized with random values of zero mean and fulfilling the conditions expressed in dimensionless form. Let us remark that the physical state so modelled consists of a large molecule which is initially at thermal equilibrium, thus the amino-acids constituting the large molecule have random configurations and movements, and then, at some initial time, "hot points" are created on the molecule, bringing it (transiently) out of equilibrium.
In Table 2 the values chosen for the physical parameters are reported. These are: the initial excitation energy E 0 , an average value of the mass M of the amino acids, the electron displacement parameter J, the elasticity constant used in the numerical studies of 10 , and the electron-phonon coupling χ . In Table 2 also the corresponding dimensionless values of the same physical quantities are reported, these are obtained by using (27) and the frequency ω = 10 13 s −1 .

Results
All the numerical computations have been performed using an integration time step t = 5 × 10 −5 entailing a very good energy conservation, with typical relative error �E/E ≃ 10 −5 . The length of the chain is N = 500 rounding the number of amino acids of the protein in 1 . Figures 1 and 2 show the spatial distribution of the probability |ψ(n, t)| 2 of finding the moving electron at any site n versus time for the electron-phonon coupling χ = 100 pN and χ = 366 pN, respectively. The electron is initially centered around the site n = 250 . Figure 1 shows that the electron wavefunction quickly spreads over the whole substate of amino acids, a phenomenon somewhat less pronounced in Fig. 2 and to some extent counterintuitive since the latter corresponds to a stronger electron-phonon coupling. Figure 3 shows As is seen from the plots in Figs. 5, the value of the phonon-phonon coupling parameter µ ′ does not seem crucial to control the release of the electron energy to the phonons, the process appears to be mainly driven by the electron-phonon coupling constant. In fact, for χ = 488 pN the relaxation time to the oscillatory state is quick and practically independent of the value of µ ′ . Neither at the lower value χ = 61 pN significative differences in the relaxation rate are observed by varying µ ′ , and even for µ ′ = 0 the energy transfer takes place in both cases of χ = 61 pN and χ = 488 pN. There are two distinct phenomena: first, the release of the electron energy to the whole ensemble of phonons, and, second, the sharing (thermalization) among the normal modes of the energy received in the first process. The rate of the first process is controlled by the coupling constant χ , and the rate of the second process is controlled by the value of µ ′ , as is clearly shown by Fig. 9.
Then we have checked how the phenomenology changes as a consequence of the introduction of the nonlinear coupling in the electron Hamiltonian. In Figs. 6 and 7 the effects of different values of the parameter ǫ are reported, again for χ = 61 pN and χ = 488 pN respectively. At χ = 61 pN the electron energy relaxation is much slower than for χ = 488 pN. In particular, for ǫ = 0 meV the electron energy appears to oscillate in a rather narrow range of values with no evident tendency of a relaxation toward a value smaller than the initial one. Whereas, for non-vanishing values of ǫ ′ the electron energy is clearly decaying in time with some oscillations. Figure 7 confirms what has already been reported in Fig. 5, that is, for χ = 488 pN the electron energy fastly decreases in time, even in the case of ǫ = 0 meV.
Let us remark that a non-vanishing value of ǫ , that is, the presence of the nonlinear coupling term in the electron Hamiltonian, plays a relevant role to ensure a more efficient transfer of part of the electron energy to the phonons of the chain of amino acids. But, in any case, it is the electron-phonon coupling constant which mainly rules the energy transfer process.
For any chosen set of physical parameters, except possibly for ǫ = 0 , the electron always transfers part of its energy to the phonons, and eventually this energy is equally shared among the phonons. In order to work out the typical time scales of this thermalization process we have computed the spectral entropy of the normal modes of the chain of amino acids, that is, of the phonons. For the harmonic term H h of the dimensionless Hamiltonian (26) we have and then, by following 18 , the coordinate transformations Q m = S mn b n and P m = S mnḃn , with where Of course, these oscillators are the normal modes (phonons) of the system. Then a spectral entropy S(t) is defined as where E T (t) = N m=1 E m (t) and E m (t) = (P 2 m + � ′ ω 2 m Q 2 m )/2 , so that the weights p m (t) are normalized. The maximum value of S(t) is attained when all the p m (t) are equal to 1/N. Thus, at equipartition, when the energy content of each normal mode is the same, entropy attains its maximum. In principle, the complete Hamiltonian for the phonon part is and the equipartition theorem now would read with f(E) a function independent of the mode. Since we have N = 500 , for each normal mode at each time step we should compute 125 × 10 6 terms in the sum for a total of 62.5 × 10 9 terms for each value of the spectral entropy at any time, and this is computationally prohibitive. On the other hand it has been shown 18,19 that even considering only the harmonic energies the thermalization process can be actually detected, even if at energy equipartition S(t) is not exactly equal to S max .
A normalized entropy is then defined as so that when the phonon oscillators are "frozen" it is S(t) = S(0) and consequently η = 1 ; but at equipartition, when S(t) = S max (t) , it is η = 0 . By following the time decay of η , it is thus possible to find out if and on which time scale the energy released by the electron is definitely transferred to the phonons. In Fig. 8 η(t) is plotted    Fig. 1. It is evident that equipartition of energy is always attained, and the time needed for this to happen is rather weakly dependent on the electron-phonon coupling constant. In fact, the equipartition rate is controlled by µ ′ as can be seen in Fig. 9. The case µ ′ = 0 is special, in the sense that the phonon-phonon coupling is indirectly made by the nonlinear electron-phonon interaction, and, in fact, the comparison between the left and right panels of Fig. 9 shows that for χ ′ = 4 the decay pattern of η(t) is suggestive of some tendency to thermalization, which is apparently absent (possibly very slow) for χ ′ = 0.5 . In fact, the decay time is approximately varying between 0.5 ns and 1 ns (the unit time scale being 10 −13 s). Importantly, the two time scales, of the electron energy release to the amino acids and of equipartition of this energy among all the normal modes of the lattice, are not equal and need not to be equal. The chaotic behavior of the particles representing the amino-acids generically prevents the reversibility of the electron energy transfer to the phonons, but then the equipartition among the phonons of the energy received from the electron depends on the phonon-phonon coupling strength and on the degree of chaoticity of the phonon dynamics 18,19 .

Discussion and conclusions
Summarizing the physical meaning of the results reported in the present work, we have to keep in mind that the parameter space of the system investigated here is very large, thus we have limited our investigation to a basic choice of physically meaningful parameters to tackle the fast process that we aimed at modelling, as stated in the Introduction. Then we have checked the robustness of the phenomenology so observed by changing some parameters, as is the case of the nonlinear coupling constants ǫ and µ , or the electron-phonon coupling constant χ . The numerical simulations of the evolution equations have shown that after having given 0.19 eV of initial excitation energy to an electron, the electron wavefunction spreads through the chain by releasing to the phonons only a small fraction of the electron energy, approximately 0.02eV. This is a somewhat unexpected result but physically interesting because it helps in understanding why exciting a collective intramolecular oscillation of the BSA protein required a very long time 1 . Of course, the contributions of several fluorophores add up, and the continuous illumination of the labelled proteins with an intense laser light allows to accumulate energy in the www.nature.com/scientificreports/ protein until the activation threshold of the coherent oscillation of all its atoms is reached and exceeded. The phonon part of the model tackled here has been simplified with respect to the model derived by the de-quantisation of the original Fröhlich's model 4 because it is focussed only on the fast mechanism of down-conversion of the energy of the photons, harvested by the protein through its fluorophore-receptors, to the internal vibrations of the chain of amino acids. Although no more than 10% of energy is dissipated by electron to phonons, in the studied regime no coherent transport of information can occur on the amino acids (as sometimes one could expect in a spin chain model 20 ) due to the fact that the electron wave function spreads over all sites. The model studied here can be easily adapted to estimate the efficiency of other excitation mechanisms of biomolecular collective oscillations, and giving an outlook at future developments, we are faced with the problem of understanding what might replace the laser action in living cells. There are several possible candidates to play the role of external energy suppliers, for instance, the hydrolysis of Adenosine Triphosphate (ATP) releases a highly energetic phosphate group that could operate a momentum transfer on some target electron via Coulomb collisions. Redox reactions and mitochondria produce weak UV photons that might excite Tryptophan and Tyrosine amino acids 21,22 in proteins, as well nucleotides of DNA and RNA. Also an anisotropic momentum transfer operated by water molecules or ions could make the job 23 . In either cases of metabolically generated photons or of ion collisions (phosphate stemming from ATP hydrolysis or other) we can assume that the external energy input for a biomolecule occurs through the generation of "hot points", as in the case of light activated fluorophores, and mediated by either radiative or collisional electronic excitation. This implicit assumption of the cumulative character of the fluorophore-contributions comes from the experimental facts observed while working out the results reported in 1 . We observed that trying to excite the collective vibrations of BSA molecules with different numbers of attached fluorophores, only with at least five attached fluorophores -in the average -the collective vibrations of the BSA molecules were activated. Let us conclude by mentioning that, for a broad class of Hamiltonian systems, long-living Quasi Stationary States (QSS) can be dynamically generated which keep a system out of thermodynamic equilibrium. Among many other systems where QSS are produced 24 , let us mention a beam of fast particles interacting with the set of waves describing a physical system 25,26 , a situation which is reminiscent, for example, of the above mentioned fast phosphate groups -produced by ATP hydrolysis.