Out-of-equilibrium collective oscillation as phonon condensation in a model protein

TL;DR

This paper models out-of-equilibrium collective oscillations in proteins as phonon condensation and confirms experimentally that proteins can exhibit such coherent oscillations when driven out of thermal equilibrium, linking theory and spectroscopy.

Contribution

It introduces a classical phonon condensation model for protein oscillations and experimentally verifies this phenomenon using THz spectroscopy on BSA proteins.

Findings

Observation of a 0.314 THz absorption feature in BSA under optical pumping

Qualitative agreement between experimental data and the phonon condensation model

Identification of collective molecular oscillation as the spectral feature

Abstract

In the first part of the present paper (theoretical), the activation of out-of-equilibrium collective oscillations of a macromolecule is described as a classical phonon condensation phenomenon. If a macromolecule is modeled as an open system, that is, it is subjected to an external energy supply and is in contact with a thermal bath to dissipate the excess energy, the internal nonlinear couplings among the normal modes make the system undergo a non-equilibrium phase transition when the energy input rate exceeds a threshold value. This transition takes place between a state where the energy is incoherently distributed among the normal modes, to a state where the input energy is channeled into the lowest frequency mode entailing a coherent oscillation of the entire molecule. The model put forward in the present work is derived as the classical counterpart of a quantum model proposed long time ago by H. Fr\"ohlich in the attempt to explain the huge speed of enzymatic reactions. In the second part of the present paper (experimental), we show that such a phenomenon is actually possible. Two different and complementary THz near-field spectroscopic techniques, a plasmonic rectenna, and a micro-wire near-field probe, have been used in two different labs to get rid of artefacts. By considering a aqueous solution of a model protein, the BSA (Bovine Serum Albumin), we found that this protein displays a remarkable absorption feature around 0.314 THz, when driven in a stationary out-of-thermal equilibrium state by means of optical pumping. The experimental outcomes are in very good qualitative agreement with the theory developed in the first part, and in excellent quantitative agreement with a theoretical result allowing to identify the observed spectral feature with a collective oscillation of the entire molecule.

Equations

$${\hat{H}}_{Tot}={\hat{H}}_0+{\hat{H}}_{Int}$$

$${\hat{H}}_0=\sum _{{\omega }_i\in {\mathcal{I}}_{sys}}\hslash {\omega }_i{\hat{a}}_{{\omega }_i}^{†}{\hat{a}}_{{\omega }_i}+\sum _{{\Omega }_j\in {\mathcal{I}}_{bath}}\hslash {\Omega }_j{\hat{b}}_{{\Omega }_j}^{†}{\hat{b}}_{{\Omega }_j}+\sum _{{\Omega }_k^{{}^{\prime }}\in {\mathcal{I}}_{src}}\hslash {\Omega }_k^{{}^{\prime }}{\hat{c}}_{{\Omega }_k^{{}^{\prime }}}^{†}{\hat{c}}_{{\Omega }_k^{{}^{\prime }}}.$$

$$\begin{array}{rl}{\hat{H}}_{IntWA}=& {\hat{H}}_{sys-bth}+{\hat{H}}_{src-sys}+{\hat{H}}_{sys-bath-sys}=\sum _{{\omega }_i\in {\mathcal{I}}_{sys},{\Omega }_j\in {\mathcal{I}}_{bth}}{\eta }_{{\omega }_i{\Omega }_j}{\hat{a}}_{{\omega }_i}^{†}{\hat{b}}_{{\Omega }_B}+\\ & +\sum _{{\omega }_i\in {\mathcal{I}}_{sys},{\Omega }_k^{{}^{\prime }}\in {\mathcal{I}}_{src}}{\xi }_{{\omega }_i{\Omega }_k^{{}^{\prime }}}{\hat{a}}_{{\omega }_i}^{†}{\hat{c}}_{{\Omega }_k^{{}^{\prime }}}+\sum _{{\omega }_{A_i},{\omega }_{A_j}\in {\mathcal{I}}_{sys},{\Omega }_k\in {\mathcal{I}}_{bth}}{\chi }_{{\omega }_i{\omega }_j{\Omega }_k}{\hat{a}}_{{\omega }_i}^{†}{\hat{a}}_{{\omega }_j}{\hat{b}}_{{\Omega }_k}^{†}+h.c.\end{array}$$

$$\begin{split}\widehat{H}_{intQ}=&\sum\limits_{\omega_{i},\omega_{j},\omega_{k},\omega_{l}\in\mathcal{I}_{sys}}\,\,\Bigr{[}\kappa_{(1)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\,\,\widehat{a}^{{\dagger}}_{\omega_{i}}\widehat{a}^{{\dagger}}_{\omega_{j}}\widehat{a}_{\omega_{k}}\widehat{a}_{\omega_{l}}+\\ &+\kappa_{(2)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\,\,\widehat{a}^{{\dagger}}_{\omega_{i}}\widehat{a}^{{\dagger}}_{\omega_{j}}\widehat{a}^{{\dagger}}_{\omega_{k}}\widehat{a}_{\omega_{l}}+\kappa_{(3)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\widehat{a}^{{\dagger}}_{\omega_{i}}\widehat{a}^{{\dagger}}_{\omega_{j}}\widehat{a}^{{\dagger}}_{\omega_{k}}\widehat{a}^{{\dagger}}_{\omega_{l}}\Bigr{]}+\mathrm{h.c.}\end{split}$$

$$|\psi ⟩=|\psi (x_1,...,x_N)⟩.$$

$$\delta S=0\text{with}S={\int }_0^{t}L(\psi ,\bar{\psi })\mathrm{d}{t}^{\prime }$$

$$L(\psi,\bar{\psi})=\dfrac{\imath\hslash}{2}\dfrac{\left\langle\psi\Big{|}\dot{\psi}\right\rangle-\left\langle\dot{\psi}\Big{|}\psi\right\rangle}{\left\langle\psi\Big{|}\psi\right\rangle}-\dfrac{\left\langle\psi\right|\widehat{H}_{Tot}\left|\psi\right\rangle}{\left\langle\psi\Big{|}\psi\right\rangle}\,\,\,.$$

$$H_{Tot}=⟨\psi (x_1,...,x_N)|{\hat{H}}_{Tot}|\psi (x_1,...,x_N)⟩.$$

$$w_{i}=\imath\hslash\left\langle\psi\Big{|}\dfrac{\partial\psi}{\partial x_{i}}\right\rangle=-\imath\hslash\left\langle\dfrac{\partial\psi}{\partial x_{i}}\Big{|}\psi\right\rangle$$

$$W_{ij}=\frac{\partial w_j}{\partial x_i}-\frac{\partial w_i}{\partial x_j}$$

$$\sum _{j=1}^N{W}_{ij}{\dot{x}}_j=\frac{\partial H_{Tot}}{\partial x_i}.$$

$$\left\{f,g\right\}=\sum _{i,j}^N\frac{\partial f}{\partial x_i}{\mathcal{W}}_{ij}\frac{\partial g}{\partial x_j}.$$

$$|\Psi (t)⟩=\prod _{{\omega }_i\in {\mathcal{I}}_{sys},{\Omega }_j\in {\mathcal{I}}_{bth},{\Omega }_k^{{}^{\prime }}\in {\mathcal{I}}_{src}}{|z_{{\omega }_i}(t)⟩}_{sys}\otimes {|z_{{\Omega }_j}(t)⟩}_{bth}\otimes {|z_{{\Omega }_k^{{}^{\prime }}}(t)⟩}_{src}$$

$$|z⟩=\exp \left[-\frac{|z{|}^2}{2}\right]\sum _{k=0}^{+\infty }\frac{z^k}{\sqrt{k!}}|k⟩=\exp \left[-\frac{|z{|}^2}{2}\right]\sum _{k=0}^{+\infty }\frac{{\left(z{\hat{a}}^{†}\right)}^k}{k!}|0⟩\text{where}z=z(t)\in \mathbb{C}.$$

$$⟨n(t)⟩=⟨z(t)|{\hat{a}}^{†}\hat{a}|(t)⟩=|z(t){|}^2$$

$$z_i=n_i^{1/2}\exp \left[-ı{\theta }_i\right]⟹z_i^{\ast }=n_i^{1/2}\exp \left[ı{\theta }_i\right]n_i=|z_i{|}^{2}i\in {\mathcal{I}}_{\mathcal{S}}.$$

$$\begin{split}&w_{n_{i}}=\imath\hslash\left\langle\left\{(n_{j},\theta_{j})\right\}_{j\in\mathcal{I}_{\mathcal{S}}}\Big{|}\dfrac{\partial}{\partial n_{i}}\left\{(n_{j},\theta_{j})\right\}_{j\in\mathcal{I}_{\mathcal{S}}}\right\rangle=\\ &=\dfrac{\imath\hslash}{\partial n_{i}/\partial|z_{i}|}\left\langle\left\{(n_{j},\theta_{j})\right\}_{j\in\mathcal{I}_{\mathcal{S}}}\right|\dfrac{\partial}{\partial|z_{i}|}\bigotimes_{j\in\mathcal{I}_{\mathcal{S}}}\left(\exp\left[-\dfrac{|z_{j}|^{2}}{2}\right]\sum_{k=0}^{+\infty}\dfrac{|z_{j}|^{k}\exp\left[-\imath k\theta_{j}\right]\left(\widehat{a}^{{\dagger}}_{j}\right)^{k}}{k!}\left|0\right\rangle\right)=\\ &=\dfrac{\imath\hslash}{2|z_{i}|}\Bigr{(}-|z_{i}|\left\langle\left\{(n_{j},\theta_{j}\right\}_{j\in\mathcal{I}_{\mathcal{S}}}\Big{|}\left\{(n_{j},\theta_{j})\right\}_{j\in\mathcal{I}_{\mathcal{S}}}\right\rangle+\exp[-\imath\theta_{i}]\left\langle\left\{(n_{j},\theta_{j})\right\}_{j\in\mathcal{I}_{\mathcal{S}}}\right|\widehat{a}^{{\dagger}}_{j}\left|\left\{(n_{j},\theta_{j})\right\}_{j\in\mathcal{I}_{\mathcal{S}}}\right\rangle\Bigr{)}=\\ &=\dfrac{\imath\hslash}{2|z_{i}|}(-|z_{i}|+|z_{i}|\exp[\imath\theta_{i}]\exp[-\imath\theta_{i}])=0\\ \end{split}$$

$$\begin{array}{rl}& w_{{\theta }_i}=ı\hslash ⟨{\left\{(n_j,{\theta }_j)\right\}}_{j\in {\mathcal{I}}_{\mathcal{S}}}|\frac{\partial }{\partial {\theta }_i}|{\left\{(n_j,{\theta }_j)\right\}}_{j\in {\mathcal{I}}_{\mathcal{S}}}⟩=\\ & =ı\hslash ⟨{\left\{(n_j,{\theta }_j)\right\}}_{j\in {\mathcal{I}}_{\mathcal{S}}}|\frac{\partial }{\partial {\theta }_i}\underset{j\in {\mathcal{I}}_{\mathcal{S}}}{⨂}\left(\exp \left[-\frac{|z_j{|}^2}{2}\right]\sum _{k=0}^{+\infty }\frac{|z_j{|}^k\exp \left[-ık{\theta }_j\right]{\left({\hat{a}}_j^{†}\right)}^k}{k!}|0⟩\right)=\\ & =ı(-ı)\hslash |z_i|\exp [-ı{\theta }_i]⟨{\left\{(n_j,{\theta }_j)\right\}}_{j\in {\mathcal{I}}_{\mathcal{S}}}|{\hat{a}}_i^{†}|{\left\{(n_j,{\theta }_j)\right\}}_{j\in {\mathcal{I}}_{\mathcal{S}}}⟩=\hslash |z_i{|}^2=\hslash n_i.\end{array}$$

$$W_{{\theta }_i{\theta }_k}=W_{n_i{n}_k}=0$$

$$W_{n_i{\theta }_k}=-W_{{\theta }_k{n}_i}=\frac{\partial w_{n_k}}{\partial {\theta }_i}-\frac{\partial w_{{\theta }_i}}{\partial n_k}=\hslash {\delta }_{i,k}.$$

$${\mathcal{W}}_{{\theta }_i{\theta }_k}={\mathcal{W}}_{n_i{n}_k}=0$$

$${\mathcal{W}}_{{\theta }_i{n}_k}=-{\mathcal{W}}_{n_k{\theta }_i}=\frac{{\delta }_{i,k}}{\hslash }.$$

$$\begin{array}{rl}{H}_0=& ⟨\Psi ({\theta }_{\omega },J_{\omega })|{\hat{H}}_0|\Psi ({\theta }_{\omega },J_{\omega })⟩=\sum _{{\omega }_i\in {\mathcal{I}}_{sys}}{\omega }_i{J}_{{\omega }_i}+\sum _{{\Omega }_k\in {\mathcal{I}}_{bth}}{\Omega }_k{J}_{{\Omega }_k}+\\ & +\sum _{{\Omega }_p^{{}^{\prime }}\in \in {\mathcal{I}}_{src}}{\Omega }_p^{\prime }{J}_{{\Omega }_p^{{}^{\prime }}}\end{array}$$

$$\begin{array}{rl}{H}_{IntWA}& =⟨\Psi ({\theta }_{\omega },J_{\omega })|{\hat{H}}_{IntWA}|\Psi ({\theta }_{\omega },J_{\omega })⟩=\\ & =\sum _{{\omega }_i\in {\mathcal{I}}_{sys}}\sum _{{\Omega }_k\in {\mathcal{I}}_{bth}}\frac{\left|{\eta }_{{\omega }_i{\Omega }_k}\right|}{\hslash }{J}_{{\omega }_i}^{1/2}{J}_{{\Omega }_k}^{1/2}\cos \left({\theta }_{{\omega }_i}-{\theta }_{{\Omega }_k}+{\theta }_{{\eta }_{{\omega }_i{\Omega }_k}}\right)+\\ & +\sum _{{\omega }_i\in {\mathcal{I}}_{sys}}\sum _{{\Omega }_p^{{}^{\prime }}\in {\mathcal{I}}_{src}}\frac{\left|{\xi }_{{\omega }_i{\Omega }_p^{{}^{\prime }}}\right|}{\hslash }{J}_{{\omega }_i}^{1/2}{J}_{{\Omega }_k^{{}^{\prime }}}^{1/2}\cos \left({\theta }_{{\omega }_i}-{\theta }_{{\Omega }_p^{{}^{\prime }}}+{\theta }_{{\xi }_{{\omega }_i{\Omega }_p^{{}^{\prime }}}}\right)+\\ & +\sum _{{\omega }_i,{\omega }_j\in {\mathcal{I}}_{sys}}\sum _{{\Omega }_k\in {\mathcal{I}}_{bth}}\frac{\left|{\chi }_{{\omega }_i{\omega }_j{\Omega }_k}\right|}{{\hslash }^{3/2}}{J}_{{\omega }_i}^{1/2}{J}_{{\omega }_j}^{1/2}{J}_{{\Omega }_k}^{1/2}\cos \left({\theta }_{{\omega }_i}-{\theta }_{{\omega }_j}+{\theta }_{{\Omega }_k}+{\theta }_{{\chi }_{{\omega }_i{\omega }_j{\Omega }_k}}\right)\end{array}$$

$$\begin{split}H_{IntQ}&=\left\langle\Psi(\theta_{\omega},J_{\omega})\right|\widehat{H}_{IntQ}\left|\Psi(\theta_{\omega},J_{\omega})\right\rangle=\sum_{\omega_{i},\omega_{j},\omega_{k},\omega_{l}\in\mathcal{I}_{sys}}J_{\omega_{i}}^{1/2}J_{\omega_{j}}^{1/2}J_{\omega_{k}}^{1/2}J_{\omega_{l}}^{1/2}\\ &\Biggr{[}\dfrac{\left|\kappa_{(1)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\right|}{\hslash^{2}}\cos\left(\theta_{\omega_{i}}+\theta_{\omega_{j}}-\theta_{\omega_{k}}-\theta_{\omega_{l}}+\theta_{\kappa_{(1)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}}\right)+\\ &\dfrac{\left|\kappa_{(2)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\right|}{\hslash^{2}}\cos\left(\theta_{\omega_{i}}+\theta_{\omega_{j}}+\theta_{\omega_{k}}-\theta_{\omega_{l}}+\theta_{\kappa_{(2)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}}\right)+\\ &+\dfrac{\left|\kappa_{(3)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\right|}{\hslash^{2}}\cos\left(\theta_{\omega_{i}}+\theta_{\omega_{j}}+\theta_{\omega_{k}}+\theta_{\omega_{l}}+\theta_{\kappa_{(3)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}}\right)\Biggr{]}\end{split}$$

$$\begin{array}{rl}& {\theta }_{{\eta }_{{\omega }_i{\Omega }_k}}={\theta }_{{\xi }_{{\omega }_i{\Omega }_p^{{}^{\prime }}}}={\theta }_{{\chi }_{{\omega }_i{\omega }_j{\Omega }_k}}={\theta }_{{\kappa }_{(1,2,3){\omega }_i{\omega }_j{\omega }_k{\omega }_k}}=0\\ & \frac{\left|{\eta }_{{\omega }_i{\Omega }_k}\right|}{\hslash }\to {\eta }_{{\omega }_i{\Omega }_K},\frac{|{\xi }_{{\omega }_i{\Omega }_p^{{}^{\prime }}}|}{\hslash }\to {\xi }_{{\omega }_i{\Omega }_p^{{}^{\prime }}},\frac{\left|{\chi }_{{\omega }_i{\omega }_j{\Omega }_k}\right|}{{\hslash }^{3/2}}\to {\chi }_{{\omega }_i{\omega }_j{\Omega }_k},\frac{\left|{\kappa }_{(1,2,3){\omega }_i{\omega }_j{\omega }_k{\omega }_l}\right|}{{\hslash }^2}\to {\kappa }_{(1,2,3){\omega }_i{\omega }_j{\omega }_k{\omega }_l}\end{array}$$

$$\begin{split}&H_{Tot}(\left\{(J_{\omega_{i}},\theta_{\omega_{i}})\right\}_{\omega_{i}\in\mathcal{I}_{\mathcal{S}}})=\sum_{\omega_{i}\in\mathcal{I}_{sys}}\,\omega_{i}J_{i}+\sum_{\Omega_{j}\in\mathcal{I}_{bth}}\,\Omega_{j}J_{\Omega_{j}}+\sum_{\Omega_{k}^{{}^{\prime}}\in\mathcal{I}_{src}}\,\Omega_{k}^{{}^{\prime}}J_{{\Omega_{k}}^{{}^{\prime}}}+\\ &+\sum_{\omega_{i}\in\mathcal{I}_{sys}}\sum_{\Omega_{j}\in\mathcal{I}_{bth}}\,\eta_{\omega_{i}\Omega_{j}}J_{\omega_{i}}^{1/2}J_{\Omega_{j}}^{1/2}\,\cos\left(\theta_{\omega_{i}}-\theta_{\Omega_{j}}\right)+\sum_{\omega_{i}\in\mathcal{I}_{sys}}\sum_{\Omega_{k}^{{}^{\prime}}\in\mathcal{I}_{src}}\,\xi_{\omega_{i}\Omega_{k}^{{}^{\prime}}}J_{\omega_{i}}^{1/2}J_{\Omega_{k}^{{}^{\prime}}}^{1/2}\cos\left(\theta_{\omega_{i}}-\theta_{\Omega_{k}^{{}^{\prime}}}\right)+\\ &\sum\limits_{\omega_{i},\omega_{j}\in\mathcal{I}_{sys}}\sum_{\Omega_{k}\in\mathcal{I}_{bth}}\chi_{\omega_{i}\omega_{j}\Omega_{k}}\,\,J_{\omega_{i}}^{1/2}J_{\omega_{j}}^{1/2}J_{\Omega_{k}}^{1/2}\cos\left(\theta_{\omega_{i}}-\theta_{\omega_{j}}+\theta_{\Omega_{k}}\right)+\\ &+\sum_{\omega_{i},\omega_{j},\omega_{k},\omega_{l}\in\mathcal{I}_{sys}}J_{\omega_{i}}^{1/2}J_{\omega_{j}}^{1/2}J_{\omega_{k}}^{1/2}J_{\omega_{l}}^{1/2}\Biggr{[}\kappa_{(1)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\cos\left(\theta_{\omega_{i}}+\theta_{\omega_{j}}-\theta_{\omega_{k}}-\theta_{\omega_{l}}\right)+\\ &\kappa_{(2)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\cos\left(\theta_{\omega_{i}}+\theta_{\omega_{j}}+\theta_{\omega_{k}}-\theta_{\omega_{l}}\right)+\kappa_{(3)\omega_{i}\omega_{j}\omega_{k}\omega_{l}}\cos\left(\theta_{\omega_{i}}+\theta_{\omega_{j}}+\theta_{\omega_{k}}+\theta_{\omega_{l}}\right)\Biggr{]}\,\,.\end{split}$$

$$⟨J_{{\omega }_i}(t)⟩=\int J_{{\omega }_i}\rho ({\left\{\left(J_{\omega },{\theta }_{\omega }\right)\right\}}_{\omega \in {\mathcal{I}}_{\mathcal{S}}},t)\prod _{\omega \in {\mathcal{I}}_{\mathcal{S}}}\mathrm{d}{J}_{{\omega }_i}\mathrm{d}{\theta }_{{\omega }_i}.$$

$$\frac{\partial \rho }{\partial t}=\left\{H,\rho \right\}=-ı{\mathcal{L}}_H(\rho )$$

$$\begin{array}{rl}\left\{f,g\right\}=& \sum _{{\omega }_i}\left(\frac{\partial f}{\partial J_{{\omega }_i}}\frac{\partial g}{\partial {\theta }_{{\omega }_i}}-\frac{\partial g}{\partial J_{{\omega }_i}}\frac{\partial f}{\partial {\theta }_{{\omega }_i}}\right)+\sum _{{\Omega }_j}\left(\frac{\partial f}{\partial J_{{\Omega }_j}}\frac{\partial g}{\partial {\theta }_{{\Omega }_j}}-\frac{\partial g}{\partial J_{{\Omega }_j}}\frac{\partial f}{\partial {\theta }_{{\Omega }_j}}\right)+\\ & +\sum _{{\Omega }_k^{{}^{\prime }}}\left(\frac{\partial f}{\partial J_{{\Omega }_k^{{}^{\prime }}}}\frac{\partial g}{\partial {\theta }_{{\Omega }_k^{{}^{\prime }}}}-\frac{\partial g}{\partial J_{{\Omega }_k^{{}^{\prime }}}}\frac{\partial f}{\partial {\theta }_{{\Omega }_k^{{}^{\prime }}}}\right).\end{array}$$

$$\left\langle f\Big{|}g\right\rangle=\int_{\lambda_{\{(J,\alpha)\}_{\omega\in\mathcal{I}_{\mathcal{S}}}}}\,f^{*}g\,\,\prod\limits_{\omega\in\mathcal{I}_{\mathcal{S}}}\mathrm{d}J_{\omega}\,\mathrm{d}\theta_{\omega}=\prod\limits_{\omega\in\mathcal{I}_{\mathcal{S}}}\int_{0}^{2\pi}\mathrm{d}\theta_{\omega}\int_{0}^{+\infty}\mathrm{d}J_{\omega}\,\,f^{*}g$$

Full text

I Introduction

Recent progress in terahertz technology has enabled to look at biological systems with terahertz radiation, that is, in an energy domain (a few meV) which is of the order of the activation energy of many biological processes. These include the excitation of collective modes of vibration of biomolecules for which experimental evidence has been originally provided at thermal equilibrium by means of Raman spectroscopy painter and, recently, by several terahertz spectroscopic studies mainly carried on using dry or low-hydrated powders because of the very strong absorption of water turton ; acbas ; falconer ; markelz . More recent studies also addressed solvated proteins xu ; ebbinghaus . Collective oscillations of biomolecules, possibly driven by metabolic activity in living matter, by bringing about large oscillating dipole moments, could activate resonant (thus selective), intermolecular attractive electrodynamic forces acting at a long distance. But we have shown Preto:2015 that this is not possible at thermal equilibrium, where all the mentioned experiments have been hitherto performed. Thus, understanding whether long-range intermolecular electrodynamic forces can be really activated requires to figure out if out-of-thermal equilibrium collective oscillations of biomolecules can be realised. The possibility of activating these forces could help explaining the fast encounters of the cognate partners of biochemical reactions in living matter, even in the context of low molecular concentration of one of the reactants, as it might be the case of some ligand-receptor reaction or of the transcription machinery of the DNA. From a theoretical viewpoint, collective oscillations of macromolecules have been described by H. Fröhlich Frohlich:1968 ; Froeh:1970 ; Froeh:1972 ; Frohlich:1977 as a Bose-like condensation of their normal vibrational modes.

However, Fröhlich’s theory has been criticized and marginalized because of several reasons. Among the others the original theory is certainly oversimplified and lacks many important explanations to qualify as really relevant to understand basic mechanisms of life. Recently, a remarkable progress was done in the direction of assessing the possible relevance of Fröhlich condensates in a revisitation of the theory reimers which has led to identify three different kinds of Fröhlich condensates: weak, strong, and coherent. The authors explain that weak condensates may have profound effects on chemical and enzyme kinetics, while coherent condensates are shown to involve very large amounts of energy, to be very fragile and not to be produced by the Wu-Austin Hamiltonian from which Fröhlich’s rate equations can be derived. The work of Ref.reimers is still theoretical and developed in a quantum framework.

From the experimental side, very interesting observations of phonon condensates of quantum nature at thermal equilibrium have been recently reported in the literature. Remarkably, optical phonon condensates have been observed in heterostructures at room temperature altfeder , arising from phase and frequency dynamical synchronisation of phonons due to suitable scattering processes around defects. Another intriguing result concerns polariton Bose condensation at thermal equilibrium in a GaAs high-Q micro cavity in a temperature range of $\displaystyle 10-25K$ polaritons , still a quite high temperature with respect to Bose-Einstein condensation of alkali atoms. These nice results concern optical phonons or polaritons weakly interacting with phonons, respectively, in condensed matter systems. But what about low frequency vibrations of biomolecules at room temperature? We can wonder whether a further and complementary step forward can be done going beyond the quantum treatment of Fröhlich condensates, also after the revisitation of Ref.reimers .

In fact, the frequency of collective oscillations of bio-macromolecules is expected in the sub-THz domain, around $\displaystyle 10^{11}$ Hz, hence at room temperature it is $\displaystyle k_{B}T/\hbar\omega\gg 1$ (where $\displaystyle k_{B}$ and $\displaystyle\hbar$ are the Boltzmann and Planck’s constants respectively) and the average number of phonons estimated at this frequency with the Bose-Einstein statistics ($\displaystyle\langle n\rangle\simeq 62.01$) is the same as that given by the Boltzmann statistics ($\displaystyle\langle n\rangle\simeq 62.51$) within a $\displaystyle 1\%$ accuracy, as a consequence a classical description seems reasonable. Moreover, when the system is brought out of thermal equilibrium by external energy injection the effective temperature of each normal mode is higher, thus increasing the occupation number of each mode and making the classical approximation even more reasonable. Furthermore, for example, the vibrational properties of proteins are very well described by Molecular Dynamics simulations performed in a classical context. Therefore, the first important question addressed in the present work is whether a phonon condensation phenomenology can be retrieved also in a classical framework and out of thermal equilibrium. To this aim, by resorting to a de-quantization method, we have worked out a classical version of the original Fröhlich model, finding that - remarkably - in a classical context too a Bose-like phonon condensation is possible classicalcondens . This possibility requires considering a biomolecule as an open system - that is, far from thermal equilibrium with its environment - through which energy flows under the simultaneous actions of an external energy supply and of dissipation due to radiative, dielectric, and viscous energy losses. We have found that the classical Bose-like condensation in the lowest vibrational mode occurs when the energy input rate exceeds some threshold value. Then, this a-priori non obvious result has motivated an experimental effort to find out if the theoretically predicted phenomenon is actually possible in the physical world. Two independent experiments, in geographically distinct laboratories, have confirmed the existence of out-of-equilibrium collective oscillations for a model protein. This is a proof of concept of which the most significant implication is that, in compliance with a theoretical prediction Preto:2015 , a crucial pre-requisite is fulfilled to excite intermolecular long-range electrodynamic interactions. In their turn, as already said above, these interactions could affect the biomolecular dynamics by contributing to drive the high efficiency and rapidity of mutual encounters of the partners of biochemical reactions in living matter, encounters which sometimes could hardly appear to be the result of Brownian diffusion only.

The derivation of the mentioned classical model - together with the numerical evidence of phonon condensation - is reported in Section II. The results of two complementary experiments are reported and discussed in Section III. Materials and methods are described in Section IV. In Section V we draw some conclusion.

II Classical out-of-equilibrium phonon condensation

Some decades ago, the study of open systems far from thermodynamic equilibrium showed, under suitable conditions, the emergence of self-organization. Striking similarities were observed among very different physical systems having in common the fact of being composed of many non-linearly interacting subsystems. When a control parameter, typically the energy input rate, exceeds a critical value (that is, when the energy gain exceeds the energy losses) then the subsystems act cooperatively to self-organize in what is commonly referred to as a non-equilibrium phase transition haken1 ; haken2 . This phenomenon is a generic consequence of the presence of the mentioned basic ingredients, that is, nonlinearly interacting subsystems, dissipation to a thermal bath, and external energy supply. Therefore, even idealised models of real systems are capable of correctly catching the occurrence of these collective behaviors, at least qualitatively. This fascinating topic was pioneered in the late ’60s by H. Fröhlich with the model mentioned in the Introduction. However, Fröhlich rate equations for the vibrational mode amplitudes, as occupation numbers in Fock’s representation, of a generic biomolecule were originally put forward heuristically: the original formulation was lacking a microscopic model. A microscopic Hamiltonian - from which the Fröhlich’s rate equations can be derived - was later given by Wu and Austin wu1 ; wu2 ; wu3 ; wu4 to model the dynamics of the normal modes of a macromolecule, of the thermal bath surrounding it, and of the external energy pump. The Wu-Austin Hamiltonian is formulated in second quantization formalism.

In what follows we indicate with $\displaystyle\widehat{a}^{{\dagger}}_{\omega_{i}},\widehat{a}_{\omega_{i}}$ the quantum creation/annihilation operators for the vibrational normal modes of the main system (i.e. a biomolecule) with frequency $\displaystyle\omega_{i}\in\mathcal{I}_{sys}$.

Such a system is put in contact with a thermal bath which represents the degrees of freedom of the environment surrounding the protein and, possibly, other normal modes of the protein which can be considered at thermal equilibrium with the surrounding environment. The thermal bath is characterized by a temperature $\displaystyle T_{B}$ and it is represented by a collection of harmonic oscillators with characteristic frequencies $\displaystyle\Omega_{j}\in\mathcal{I}_{bth}$ whose annihilation/creation operators are $\displaystyle\widehat{b}_{\Omega_{j}}$ and $\displaystyle\widehat{b}^{{\dagger}}_{\Omega_{j}}$, respectively.

In order to put the system representing normal modes of a biomolecule out of thermal equilibrium, an external source of energy is necessary: such an external source is represented as another thermal bath at a temperature $\displaystyle T_{S}\gg T_{B}$. Also in this case the corresponding thermal bath is described by a collection of harmonic oscillators with frequencies $\displaystyle\Omega_{k}^{{}^{\prime}}\in\mathcal{I}_{src}$, the quantum annihilation/creation operators of which are $\displaystyle\widehat{c}_{\Omega_{k}^{{}^{\prime}}}$ and $\displaystyle\widehat{c}^{{\dagger}}_{\Omega_{k}^{{}^{\prime}}}$. These three sets of harmonic oscillators can be regarded as three subsystems of a larger isolated system $\displaystyle\mathcal{S}$ (we coherently indicate with $\displaystyle\mathcal{I}_{\mathcal{S}}$ the set of all the normal modes of the system) whose quantum dynamics is described by the Hamiltonian

[TABLE]

where $\displaystyle\widehat{H}_{0}$ is the free Hamiltonian of the three sets of harmonic oscillators representing the molecular normal modes and the two heat baths

[TABLE]

The interactions among normal modes are described by $\displaystyle\widehat{H}_{Int}$; in the original formulation by Wu and Austin the interaction term has the form:

[TABLE]

where $\displaystyle\eta_{\omega_{i}\Omega_{j}},\xi_{\omega_{i}\Omega_{k}^{{}^{\prime}}},\chi_{\omega_{i}\omega_{j}\Omega_{k}}\in\mathbb{C}$ are the coupling constants describing the linear interactions among the thermal bath modes and the biomolecule modes, the linear interactions between the external source and the biomolecule, and the mode-mode interactions among the biomolecule normal modes mediated by the thermal bath, respectively.

From these terms it is possible to derive the Fröhlich rate equations, by resorting to time dependent perturbation theory; details are given in Refs. wu1 ; wu2 and in the reference book in Ref. pokorny .

However, the mode coupling term corresponds to a potential energy unbounded from below and, consequently, this would give rise to dynamical instability of the system and, in the quantum context, to the absence of a finite energy ground state. This problem led to strong criticism against the Wu-Austin Hamiltonian model and also against the ensemble of Fröhlich condensation theory bolterauer . This problem can be easily fixed by adding a term with a quartic dependence on the creation/annihilation operators of the form

[TABLE]

so that the lower bound to the ground energy level does not go to $\displaystyle-\infty$ for large values of $\displaystyle N_{\omega_{i}}$.

This quartic interaction stands for a anharmonic interaction among the normal modes of the biomolecule, a broadly studied topic of relevance to energy transport in biomolecules like proteins.

II.1 Dequantization of Wu-Austin Hamiltonian by Time Dependent Variational Principle (TDVP)

In this section, a classical Hamiltonian system with its canonical equations is associated to the quantum system described by the quantum Hamiltonian of eqs. (1), (3) and (4). This result is obtained by applying the Time Dependent Variational Principle (TDVP) technique kramer1980 ; kramer2008 . The same technique, more deepened from both the mathematical and conceptual points of view, has been proposed as a “dequantization” technique in jauslin as a kind inverse procedure with respect to the geometrical quantization. This consists in the evaluation of the time-dependent operator action on coherent states of quantum harmonic oscillators. The scalar parameters describing the coherent states become generalized coordinates of a classical dynamical system whose equations of motion can be derived from a variational principle.

In more details, one begins with the ansatz that the wavefunction depends on $\displaystyle N$ parameters $\displaystyle\{x_{i}\}_{i=1,...,N}$

[TABLE]

where the parameters $\displaystyle x_{i}=x_{i}(t)$ are in general functions of time. For a quantum system with Hamiltonian $\displaystyle\widehat{H}_{Tot}$ the equations of motion of the $\displaystyle x_{i}$ can be derived using the following variational principle (equivalent to the Least Action principle)

[TABLE]

where $\displaystyle L(\psi,)$ is the Lagrangian associated to the system

[TABLE]

The equations of motions derived from Eq.(6) can be worked out in the framework of classical Hamiltonian dynamics.

The classical Hamiltonian is associated with the quantum one by simply taking the expectation value of the Hamiltonian operator $\displaystyle\widehat{H}_{tot}$ over the state $\displaystyle\left|\psi(x_{1},...,x_{N})\right\rangle$, that is

[TABLE]

The Poisson brackets $\displaystyle\{\cdot,\cdot\}$ depend only on the chosen parametrization for the wavefunction. Starting from the variables

[TABLE]

one defines the antisymmetric tensor

[TABLE]

so that the equations of motion are implicitly given by

[TABLE]

If the condition $\displaystyle\mathrm{Det}W_{ij}\neq 0$ holds, then the matrix $\displaystyle\mathcal{W}_{ij}=\left(W^{-1}\right)_{ij}$ defines the Poisson brackets for the classical Hamiltonian system

[TABLE]

This formalism can be applied to the quantum system described by the quantum Hamiltonian of Eq.(1) to associate to it a classical Hamiltonian system. The choice of the parametrization for the wavefunction is quite arbitrary and the TDVP, as any other variational principle, restricts the dynamics to a certain region of the Hilbert space. Since the Hamiltonian is expressed in terms of creation/annihilation operators of the quantum harmonic oscillators describing the system, the wavefunction is been chosen as a product of the corresponding coherent states. In particular

[TABLE]

where $\displaystyle\left|z_{A_{i}}(t)\right\rangle_{sys},\left|z_{B_{k}}(t)\right\rangle_{bth},\left|z_{C_{p}}(t)\right\rangle_{src}$ are normalized coherent states for the normal modes of the main system, of the thermal bath and of the external source, respectively: their general form is given by

[TABLE]

From the definition of coherent states in Eq.(14) it follows that the expectation value for the occupation number $\displaystyle n_{i}$ is given by the squared norm of $\displaystyle z$

[TABLE]

so, as we are interested in writing rate equations for these quantities, we parametrize the wavefunction $\displaystyle\psi(t)$ with the set of real parameters $\displaystyle\left\{(n_{i},\theta_{i})\right\}_{i\in\mathcal{I}_{\mathcal{S}}}$ such that

[TABLE]

Using Eq.(9) it is possible to derive the Poisson brackets associated with the variables $\displaystyle\left\{(n_{i},\theta_{i})\right\}_{i\in\mathcal{I}_{\mathcal{S}}}$:

[TABLE]

[TABLE]

and, consequently, using the definition (10) the entries of the matrix $\displaystyle W$ are

[TABLE]

and its inverse has the form

[TABLE]

Consequently it follows that the variables $\displaystyle J_{\omega}=\hslash n_{\omega}$ and $\displaystyle\theta_{\omega}$ are canonically conjugated variables. The classical Hamiltonian $\displaystyle H=H_{0}+H_{IntWA}+H_{IntQuad}$ for the variables $\displaystyle\left\{(\theta_{\omega},J_{\omega})\right\}_{\omega\in\mathcal{I}_{\mathcal{S}}}$ is given by a free classical part

[TABLE]

by a semi-classical Wu and Austin interaction part

[TABLE]

and by the quartic term

[TABLE]

where each complex coupling constant is given in polar representation. In what follows, the coupling constants are considered real and rescaled s.t.

[TABLE]

with these choices, the total Hamiltonian of the system reads

[TABLE]

In order to derive Fröhlich-like rate equations, the dynamics of the action variables $\displaystyle J_{\omega_{i}}$ of the system has to be studied. We could choose to investigate the dynamics of the system by letting observable quantities evolve in time (according to Hamilton’s equations of motion) and then performing time averaging, and averaging on different initial conditions compatible with the assumption that the two subsystems $\displaystyle\mathcal{I}_{bth}$ and $\displaystyle\mathcal{I}_{src}$ are two thermal baths with different temperatures. However, this method has some disadvantages: the integration should be done numerically because of the presence of non-linear interactions terms in the Hamiltonian, and a very large number of degrees of freedom would be necessary to adequately simulate the dynamics of the thermal baths. Moreover, long integration times would be necessary to attain the convergence of time averages of the observables, and, finally, this computational effort would not be worth as it would provide a redundant information on the dynamics. Another way to derive rate equations for the $\displaystyle J_{\omega_{i}}$ consists in a statistical approach: the relevant dynamical variable to consider is the phase space distribution function $\displaystyle\rho(\left\{(J_{\omega},\theta_{\omega})\right\}_{\omega\in\mathcal{I}_{\mathcal{S}}},t)$ so that the rate equations are written for the statistical averages of actions variables

[TABLE]

In the following Section classical Fröhlich-like rate equations are derived by resorting to the time evolution of the distribution function $\displaystyle\rho(\left\{(J_{\omega},\theta_{\omega})\right\}_{\omega\in\mathcal{I}_{\mathcal{S}}},t)$ satisfying the Liouville equation.

II.2 Derivation of classical rate equations using the Koopman-Von Neumann formalism

Let $\displaystyle\rho(\left\{(J_{\omega},\theta_{\omega})\right\}_{\omega\in\mathcal{I}_{\mathcal{S}}};t)$ be a probability density function for the whole system described by the Hamiltonian in Eq.(27); according to the Liouville Theorem the evolution of $\displaystyle\rho$ associated with this Hamiltonian is given by

[TABLE]

where $\displaystyle\imath^{2}=-1$ and $\displaystyle\left\{\cdot,\cdot\right\}$ are the canonical Poisson brackets

[TABLE]

and $\displaystyle\mathcal{L}_{H}(\cdot)=\imath\left\{H,\cdot\right\}$ is the Liouville operator acting on functions defined on the phase space of the system. An interesting method to study and solve Liouville equations relies on the Koopman-Von Neumann (KvN) formalism developed in the 30’s: a formal analogy among Liouville and Schrödinger equation is established so that also classical mechanics can be formulated in the framework of a Hilbert space of square integrable functions.

In our case the Hilbert space of complex square integrable functions in phase space is $\displaystyle L^{2}(\Lambda_{\{(J_{\omega},\theta_{\omega})\}_{\omega\in\mathcal{I}_{\mathcal{S}}}})$ with the inner product defined by

[TABLE]

with $\displaystyle f,g\in L^{2}(\Lambda_{\{(J_{\omega},\theta_{\omega})\}_{\omega\in\mathcal{I}_{\mathcal{S}}}})$. On this space we can define the action of the Liouville operator

[TABLE]

and consider the domain $\displaystyle\mathcal{D}_{\widehat{\mathcal{L}}_{H}}\subseteq L^{2}(\Lambda_{\{(J_{\omega},\theta_{\omega})\}_{\omega\in\mathcal{I}_{\mathcal{S}}}})$ where the Liouville operator is self-adjoint, namely, $\displaystyle\widehat{\mathcal{L}}^{{\dagger}}_{H}=\widehat{\mathcal{L}_{H}}$ and $\displaystyle\mathcal{D}_{\widehat{\mathcal{L}}_{H}}=\mathcal{D}_{\widehat{\mathcal{L}}^{{\dagger}}_{H}}$.

Let $\displaystyle\psi(\{(J,\theta)\}_{\mathcal{I}_{\mathcal{S}}};t)\in\mathcal{D}_{\widehat{\mathcal{L}}_{H}}$ be a normalized time-dependent function nota such that

[TABLE]

then it can be proved that $\displaystyle\rho=\|\psi\|_{L^{2}(\Lambda_{(J,\theta)})}=\psi^{*}\psi$ is a normalized function for which (29) holds. Moreover as $\displaystyle\widehat{\mathcal{L}}_{H}$ is a self-adjoint operator, it represents the unitary time evolution of the “wave function” as

[TABLE]

in analogy with quantum mechanics. With this formalism the rate equations for the average values of the actions (which are the counterpart of quantum occupation numbers) associated with the normal modes $\displaystyle\omega_{i}$ of the main system are given by

[TABLE]

where $\displaystyle\widehat{\mathcal{M}}_{J_{\omega_{i}}}$ is a multiplicative operator acting on $\displaystyle L^{2}(\Lambda_{\{(J_{\omega},\theta_{\omega})\}_{\omega\in\mathcal{I}_{\mathcal{S}}}})$ as follows

[TABLE]

The Liouville operator can be decomposed as

[TABLE]

and since the eigenfunctions of the operator $\displaystyle\widehat{\mathcal{L}}_{H_{0}}$ are known, the action of the operator $\displaystyle\widehat{\mathcal{L}}_{H_{Int}}$ can be treated as a time dependent perturbation, which is adiabatically turned on and off from $\displaystyle t_{0}=0$ to $\displaystyle t=+\infty$. Then a classical analogous of the interaction representation formalism is used. Thus, if $\displaystyle\left|\psi(t)\right\rangle_{S}$ is the “wave function” in the Schrödinger representation, then in interaction representation $\displaystyle\left|\psi(t)\right\rangle_{I}$ reads

[TABLE]

and, given a generic operator $\displaystyle\widehat{A}_{S}$ in the Schrödinger picture, its expression in the interaction picture $\displaystyle\widehat{A}_{I}$ ris

[TABLE]

With this formalism the time evolution of $\displaystyle\left|\psi(t)\right\rangle_{I}$ can be written through the unitary evolution operator $\displaystyle\widehat{U}(t;t_{0})$ satisfying

[TABLE]

and the formal solution of (41) is given by

[TABLE]

At first order in $\displaystyle\widehat{\mathcal{L}}_{H_{int}}(t^{\prime})$, the unitary evolution operator $\displaystyle\widehat{U}(t;t_{0})$ in the right hand side of eq. is substituted by the identity operator meaning that the state $\displaystyle\left|\psi(t_{0})\right\rangle_{I}$, if the perturbation is turned on at $\displaystyle t_{0}$, at the lowest order can be approximated by $\displaystyle\left|\psi(t_{0})\right\rangle_{I}\simeq\left|\psi(t_{0})\right\rangle_{S}$ and assumed to be coincident with Schrödinger picture (i.e. $\displaystyle\left|\psi(0)\right\rangle=\left|\psi_{0}\right\rangle$), so that

[TABLE]

and, as $\displaystyle\widehat{\mathcal{L}}_{H_{int}}(t^{\prime})$ is self-adjoint, the time evolution for the ”bra” has the form

[TABLE]

The time derivative of the multiplicative operator $\displaystyle\widehat{\mathcal{M}}_{J_{\omega_{i}}}$ in interaction picture is derived according to Eq.(39)

[TABLE]

This means that the average (35) can be entirely rewritten using the interaction picture as

[TABLE]

After lengthy computations (see SM for the details) one finally arrives at the following rate equations for the expectation values of the actions (that is, the amplitudes of vibrational modes) which are worked out as follows

[TABLE]

II.3 Properties of the classical rate equations

The rate equations (47) that have been derived by solving the Liouville equation with the Hamiltonian (27) display a functionally similar structure to the original Fröhlich’s equations Frohlich:1968 , apart from the additional quartic terms. As in the original formulation given by Fröhlich, the condensation phenomenon is found by considering the stationary solutions of the rate equations. It is convenient to rewrite Eqs.(47) in non-dimensional form after introducing the following variables

[TABLE]

so that Eqs.(47) read

[TABLE]

By inspection, one can notice the following properties of the equations above

• •

if the normal modes thermalize at the heat bath temperature $\displaystyle T_{B}$, so that $\displaystyle\left\langle J_{\omega_{i}}\right\rangle=k_{\text{B}}T_{B}/\omega_{i}$, it follows that the variables $\displaystyle y_{\omega_{i}}$ are equal and take the value $\displaystyle y_{\omega_{i}}=1$;

• •

by switching-off the external source of energy, i.e. putting $\displaystyle S_{\omega_{i}}=0$, the thermal solution above, that is $\displaystyle y_{\omega_{i}}=1$ for all the $\displaystyle\omega_{i}$, is a stationary solution of the system, namely $\displaystyle\dot{y}_{\omega_{i}}=0$.

In order to understand whether the phonon condensation phenomenon exists also in the above defined classical framework, one has to work out non trivial out of equilibrium stationary states of the model equations. Unfortunately, doing this analytically for a system with a large number of nonlinear equations is hopeless, therefore one has to resort to numerical integration of the dynamical equations (49). A difficulty of the numerical approach is to provide a priori estimates of the coupling constants for a real biomolecule. To overcome this problem we borrowed from Ref.piazza the values of physically realistic nonlinear coupling parameters. Moreover, the non-dimensional form of the equations partially simplifies the problem as only the ratios among the coefficients has to be known. On the other hand, since we are interested only in a qualitative study of purely theoretical kind, some simplifying assumptions are in order. Thus, the coupling constants appearing in Eqs.(49) are supposed to be independent of the frequency, in particular, also the condition $\displaystyle S_{\omega_{i}}=S$ is assumed (the energy injection rate is the same for all the normal modes). Then we have performed the numerical integration of the dynamical equations by changing the total number $\displaystyle N+1$ of oscillators (normal modes) in a fixed range ($\displaystyle 1\leq\omega_{i}\leq 2$) of equally spaced frequencies given (in arbitrary units) by $\displaystyle\omega_{n}=1+n/N$ with $\displaystyle n=0,1,2,..,N$. This choice models an increasing density of modes in a given fixed frequency interval. The classical energy condensation phenomenon, numerically found and described below, is characterised by a strong deviation from energy equipartition among the normal modes in a stationary state of the dynamical system (47) (and, equivalently, of the system (49) ): as the energy injection rate increases, the system undergoes a major change in the energy distribution among its normal modes resulting in a more “organized” phase (the energy is mainly channeled into the lowest frequency mode). Such a condensation phenomenon can be viewed as associated with the breaking of the permutation symmetry for the set of variables $\displaystyle\left\{y_{\omega_{0}},...,y_{\omega_{N}}\right\}$ in the stationary state, where $\displaystyle y_{\omega_{i}}$ is the energy content in the mode of frequency $\displaystyle\omega_{i}$ in $\displaystyle k_{\text{B}}T_{B}$ units. It is clear that such a symmetry holds at energy equipartition. In analogy with the characterization of a symmetry breaking in the classical theory of equilibrium phase transitions, we have to define some (possibly scalar) indicators to detect the permutation symmetry breaking. In order to detect energy equipartition, and its violation in presence of the condensation phenomenon, we borrow and adapt from Ref.reimers a condensation index $\displaystyle\mathcal{E}_{y}$ which reads

[TABLE]

where $\displaystyle y_{\omega_{0}}$ is the energy content of the fundamental mode and $\displaystyle{\widetilde{y}_{\omega_{0}}}=1+S/B$ is the energy content of the fundamental mode given by the exact solution of Eqs.(49) in the absence of all the nonlinear couplings (mode-mode and system-bath). Here we replace the occupation numbers of the modes entering the condensation index of Ref.reimers with the mode energy contents. $\displaystyle\mathcal{E}_{y}=0$ corresponds to energy equipartition among the normal modes, and $\displaystyle\mathcal{E}_{y}=1$ corresponds to full concentration of the energy into the lowest frequency mode. In order to easily compare the results obtained at different number $\displaystyle N$ of modes, the energy fractions $\displaystyle y_{\omega_{i}}$ in each normal mode are normalised to give the variables $\displaystyle p_{i}=y_{\omega_{i}}/\sum_{\omega_{i}\in\mathcal{I}_{sys}}y_{\omega_{i}}$ .

Another possible choice for the indicator detecting the permutation symmetry breaking is any of the quantities $\displaystyle y_{\omega_{i}}/y_{\omega_{0}}=p_{i}/p_{0}$ with $\displaystyle i\neq 0$: these are expected to attain their maximum value equal to one at equipartition and to take values between zero and one for a Fröhlich-like condensation phenomenon. Numerical simulations have been performed in order to characterize this major change of the stationary states of the dynamical system (49), that is, to find the asymptotic stationary solutions of the variables $\displaystyle y_{\omega_{i}}(t)$ which correspond to the fixed points $\displaystyle\dot{y}_{\omega_{i}}=0$ of the system. As a reasonable initial condition, at $\displaystyle t=0$ the system is taken at thermal equilibrium with the heath bath, that is, all the $\displaystyle y_{\omega_{i}}(0)$ are set equal to $\displaystyle 1$.

II.4 Results of numerical simulations

In Figs.1 and 2 the results of numerical integration of the rate equations (49) are reported for the set of parameters $\displaystyle B_{\omega_{i}}=B=1$, $\displaystyle C_{\omega_{i}\omega_{j}}=C=0.1$ with the initial conditions specified in the previous subsection. The coefficient $\displaystyle B$ has been chosen as reference parameter because it is directly related with the characteristic thermalization time scale $\displaystyle\tau_{\text{therm}}\approx B^{-1}$. The quartic coupling constant has been chosen $\displaystyle\tilde{\Upsilon}_{(1)}=\tilde{\Upsilon}_{(2)}=10^{-4}$, which is sufficient to guarantee a finite bound from below of the energy. The choice $\displaystyle C/B=0.1$ is coherent with respect to what has been reported in the literature pokorny concerning the already investigated qualitative aspects of quantum Fröhlich condensation. The integration has been performed using a fourth-order Runge-Kutta algorithm, in each case for a time interval $\displaystyle\tau_{\text{int}}>0$ sufficiently long so as to guarantee that $\displaystyle|\epsilon_{\omega_{i}}(\tau)-\epsilon_{\omega_{i}}(\tau_{\text{stat}})|<10^{-5}$ for $\displaystyle\tau_{\text{stat}},t\in[0,t_{\text{int}}]$. In Figure 1 the reported results display the fraction of energy content in each mode for a system with $\displaystyle N=20$, that is $\displaystyle 21$ normal modes and different values of the energy injection rate. It appears quite clearly that the larger the energy injection rate $\displaystyle S$ the stronger the deviation from energy equipartition with increasing energy concentration in the lowest frequency mode $\displaystyle\omega_{0}$. This energy concentration phenomenon, that we can call phonon condensation, stems from classical Liouville equation and displays several analogies with the phenomenology of the pseudo-Bose condensation predicted by Fröhlich in a quantum framework, and subsequently largely studied and debated in the same quantum framework.

A relevant feature of the Fröhlich condensation is the existence of a threshold in the energy injection rate which corresponds to the breaking of the above discussed permutation symmetry. In order to find out if a similar threshold effect is displayed also by our classical model, we studied the behavior of the above defined indicators of symmetry breaking (the condensation index $\displaystyle\mathcal{E}_{y}$ and the ratio $\displaystyle y_{\omega_{i}}/y_{\omega_{0}}$) as a function of $\displaystyle S$ and for different values of $\displaystyle N$. Figure 1 shows that, for the chosen values of $\displaystyle S$, it is $\displaystyle p_{i}>p_{i+i}$ hence the parameter $\displaystyle p_{1}/p_{0}$ appears as a natural choice to detect a possible threshold effect. The results of numerical simulations of the order parameters have been reported in Figure 2. Both indicators pass from their respective equipartition values ($\displaystyle\mathcal{E}_{y}\simeq 0$, $\displaystyle p_{1}/p_{0}\simeq 1$) to their condensation values ($\displaystyle\mathcal{E}_{y}\simeq 1$, $\displaystyle p_{1}/p_{0}\simeq 0$) by displaying a clear steepening at increasing $\displaystyle N$. In particular, at increasing $\displaystyle N$, the curves of the condensation index $\displaystyle\mathcal{E}_{y}(S,N)$ display a tendency to superpose at large $\displaystyle S$-values, and a sharpening of the transition knee at lower $\displaystyle S$-values, suggesting an asymptotic convergence to a sharp bifurcation pattern. In fact, this sharpening of the transition curves with $\displaystyle N$ is reminiscent of the finite-size patterns of the order parameter as a function of the control parameter of equilibrium second order phase transitions. An analogous phenomenon is displayed by the $\displaystyle S$-pattern $\displaystyle(p_{1}/p_{0})(S)$ of the indicator $\displaystyle p_{1}/p_{0}$ at increasing $\displaystyle N$. Actually, the right panel of Figure 2 shows that, around the value $\displaystyle S\simeq 10$, the different curves $\displaystyle(p_{1}/p_{0})(S,N)$ cross each other and get steeper at increasing $\displaystyle N$, again suggesting that the larger the numbers of modes the sharper the transition between almost-equipartition among the modes and phonon condensation in the lowest frequency mode. It is worth noting that this phenomenology is qualitatively confirmed by the outcomes reported in Figure 5 (aimed at highlighting a saturation effect) where again a sharpening - at increasing $\displaystyle N$ - of a threshold effect can be observed. The lower panel of Figure 2, for the case of $\displaystyle N_{sys}=201$, provides some interesting and more detailed information about the dynamics underlying the condensation phenomenon reported in the upper panels of the same Figure. In fact, the lower panel shows how the relative energy contents $\displaystyle p_{i}$ of a randomly chosen subset of modes ($\displaystyle p_{0},p_{1},p_{5},p_{10},p_{100}$) deviate from energy equipartition at increasing $\displaystyle S$. We can observe that the highest frequency mode selected, $\displaystyle p_{100}$, deviates from equipartition at smaller values of $\displaystyle S$ than the lower frequency modes chosen; then, somewhat below $\displaystyle S=10$, the energy content of all the selected modes starts decreasing except the energy content $\displaystyle p_{0}$ of the lowest frequency mode which gathers more and more energy approaching the value $\displaystyle p_{0}=1$, that is, full condensation at high $\displaystyle S$.

The main result of this theoretical part of the paper is an original derivation of a classical counterpart of the well-known Wu-Austin quantum Hamiltonian from which Fröhlich’s rate equations can be derived. This classical counterpart consists of a system of nonlinear dynamical equations (49) for the (adimensionalized) energy contents $\displaystyle y_{\omega_{i}}$ of each normal mode $\displaystyle\omega_{i}$ of an ensemble of oscillators coupled with two thermal baths at different temperatures. The ensemble of oscillators - ideally modelling the vibrations of the atoms or of groups of atoms in a macromolecule - is both linearly and nonlinearly coupled to the set of oscillators representing a thermal bath at temperature $\displaystyle T_{B}$; this thermal bath accounts for the environment (mainly water molecules) and for other fast internal oscillations of a macromolecule. The second thermal bath at a temperature $\displaystyle T_{S}\gg T_{B}$ models an external energy supply, in close analogy with the Wu-Austin model. The numerical study of Eqs. (49) has shown that the condensation phenomenon is not peculiar of a quantum model, which a-priori was not an obvious fact, hence the occurrence of a phonon condensation phenomenon in a macromolecule at room temperature is more plausible on the basis of our classical model than on the basis of Fröhlich’s quantum model. Therefore, even though Equations (49) still represent a biomolecule in a very idealized way, the genericity of qualitative features of out-of-equilibrium phase transitions - mentioned at the beginning of this Section - suggests that the condensation phenomenon should actually occur in real biomolecules, even though our model is not expected to predict the quantitative details of the phenomenon. This is enough to make the predicted phenomenology worth an experimental effort to detect it.

III Experimental detection of non-equilibrium collective mode

Even though experimental evidence of the existence of collective modes of vibration of biomolecules has been provided at thermal equilibrium by means of Raman spectroscopy painter already many years ago, and is still being the object of many investigations turton ; acbas ; falconer ; markelz ; xu ; ebbinghaus ; lundholm , no experimental evidence was hitherto available of the possibility of exciting out-of-thermal-equilibrium collective oscillations of a biomolecule. Unveiling whether these can be activated amounts to understanding whether a necessary condition to activate long-range intermolecular electrodynamic forces Preto:2015 can be fulfilled. This is what strongly motivated the joint theoretical and experimental work reported in the present paper. The experimental counterpart of the theoretical work relies on two complementary experiments.

For both experiments a model protein has been chosen: the BSA (Bovine Serum Albumine) protein. This is mainly made out of $\displaystyle\alpha$-helices, and is a ”model” since it is largely studied in the biophysical chemistry literature. Our strategy to create a stationary out-of-thermal-equilibrium state of this molecule is to induce it by means of optical pumping, without involving any optical transition of the protein, through the excitation of some fluorochromes bound to each protein molecule. The optical excitation of these fluorochromes creates on each protein some “hot points” acting as the epicentres of a so-called “proteinquake” ansari ; levantino - better discussed in the following - and resulting in an energy transfer to the vibrational part of the protein. We used the Alexa488 fluorochrome which is covalently bonded at the lysine residues of the BSA and which is excited by means of an Argon laser (wavelength $\displaystyle 488$ nm). Some $\displaystyle 0.19$ eV per fluorochrome and per incident photon (the average energy difference between the absorbed and re-emitted photons) is thus available for an energy transfer to the protein and, partly, also to its environment. By attaching an average number of 5 fluorochromes per protein a considerable amount of energy ( $\displaystyle\gg k_{B}T$) can be continuously pumped into each protein. Two THz-near-field absorption spectroscopy setups of aqueous solutions of the protein (at 1 mg/mL concentration) operating into two distinct laboratories, have been used at room temperature (Fig. 3 (a), (c)). In both experiments THz radiation is produced by tunable, highly-spectrally-resolved ($\displaystyle<300$ Hz) and continuous-wave sources with an average power of 1 mW allowing an accurate detection of possible resonances. A typical experiment consists in three phases, during all of them the aqueous solution of proteins is illuminated with THz radiation performing a sweep in frequency, and thus allowing to measure the frequency dependence of the absorbed electromagnetic power (detected by the near-field probes) in the solution. During the first phase no extra illumination with the Argon laser is done; during the second phase the Argon laser is switched on to excite the fluorochromes bound to the proteins; finally, during the third phase the laser is switched off to check whether some memory and irreversible photochemistry effect (photobleaching) or sample heating are present. The use of near-field coupling of metal probes to the sample eliminates the Fabry-Perot interferences often seen in optical spectra taken in fluidic cells globus .

The first setup (Fig. 3 (a), (b)) used a micro-coaxial near-field probe put inside a metallic rectangular waveguide enabling a modal transition from $\displaystyle TM_{01}$ Sommerfeld’s to $\displaystyle TE_{01}$ waveguide mode. The sub-wavelength diameter of the wire (12 $\displaystyle\mu$m) allows an extremely localized detection of the longitudinal component of the electric field at its end and on a volume of about 4 pL. The spectra were subtracted from the spectrum of pure water in order to remove artifacts coming either from the water absorption or from the geometry of the experimental setup. The second setup giliberti (Fig. 3 (c), (d)) used a near-field probe rectenna composed of a planar metal bow-tie antenna with dimensions close to half-a-wavelength (at 0.3 THz) that enhanced the THz field in the feed gap region (volume of about 0.2 pL) and a plasma-wave field-effect transistor (FET) integrated in the feed gap of the antenna. When illuminated by THz radiation, the antenna-coupled FET device provides a DC-voltage – between Source and Drain contacts – proportional to the THz-field intensity dyakonov ; knap ; nouvel . A hemispherical silicon lens pressed on the back of the semiconductor substrate focused the THz radiation on the antenna simultaneously eliminating Fabry-Perot interference in the substrate digaspare . The spectrum of the protein solution obtained in the absence of blue light illumination was subtracted from the spectrum obtained with blue light illumination.

Figure 4 (a) presents the spectra obtained using the micro-coaxial probe in the absence of blue light illumination (black squares) and in the presence of blue light illuminations of different durations (from 3 to 9 min). In the former case (no illumination), there was no specific spectral feature in the studied frequency range while in the latter case (with illumination) we observed spectral resonances which become more evident for increasing duration of illumination. In particular, the strongest resonance appeared at 0.314 THz accompanied by two other minor resonances situated at 0.278 and 0.285 THz; these values did not depend on the time of illumination and the strength of the resonances saturated after 9 min of illumination. These results have full reproducibility. Figure 4 (b) presents the spectra obtained using the rectenna probe for two durations of blue light illumination. The spectra with illumination-on versus illumination-off were taken several times (two of them reported to show reproducibility). Also in this case we observed the appearance of evident resonances whose strength saturated at increasing durations of illumination (10 and 15 minutes in the two experimental runs reported). Long excitation times are needed because, under our experimental conditions, the energy dissipation rate and the energy supply rate are almost equal, so that a long time is needed to accumulate enough energy into each protein in order to make intramolecular nonlinear interaction terms sufficiently strong to activate the condensation phenomenon.

The spectra obtained using the two previously described methods and for the longer durations of illumination are compared in Fig. 4 (c). The principal resonance at 0.314 THz is perfectly reproduced using two completely different and complementary setups. Since THz extinction in water is huge ($\displaystyle 2000$ dB/cm), the emergence of this spectral feature of the protein out of the water absorption background must be associated with the activation of a giant dipole moment. And this can happen only as a consequence of the activation of a coherent oscillation of the whole molecule, possibly together with a dipole moment strengthening due to an electret-like structuring of water dipoles of the first hydration layers of the proteins as hydration shell might contribute to the magnitude of the protein dipole moment susko .

Since the BSA is a heart-shaped globular protein, a possible lowest frequency of a global oscillatory mode is roughly estimated by schematising the molecule as composed of two masses $\displaystyle m$, equal to half the total protein mass, joined by a spring of elastic constant $\displaystyle k$ given by $\displaystyle k=EA_{0}/l_{0}$, where $\displaystyle E$ is the Young modulus of the protein, $\displaystyle A_{0}$ and $\displaystyle l_{0}$ are its transverse section and length at rest, respectively. Using $\displaystyle m=33$ kD, $\displaystyle A_{0}\simeq 1.2\times 10^{-13}$ cm2, $\displaystyle l_{0}\simeq 1.2\times 10^{-7}$ cm, and $\displaystyle E=6.75\ \rm{GPa}$, we find $\displaystyle\nu=(1/2\pi)\sqrt{k/m}\simeq 0.300$ THz which is close to the main resonance observed at $\displaystyle 0.314$ THz. Since the BSA molecule can be modelled to first order as a three-dimensional elastic nanoparticle xu , a more refined approximation is obtained by modelling the protein with an elastic sphere and then considering its vibrational frequencies.

The fundamental frequency of a spheroidal deformation mode of an elastic sphere is given by the formula bastrukov

[TABLE]

which holds for $\displaystyle l\geq 2$. Using the following data for the BSA protein: Young modulus $\displaystyle E=6.75\ \rm{GPa}$ obtained at room temperature using Brillouin light scattering of hydrated BSA proteins perticaroli , hydrodynamic (Stokes) radius $\displaystyle R_{H}=35\textrm{\AA}$, and specific volume $\displaystyle 1/\rho=0.74$ derived from small-angle X-ray scattering (SAXS) experiments mylonas , we find for the lowest mode ($\displaystyle l=2$) the frequency $\displaystyle\nu_{0}=0.308\ THz$ which agrees within an error of about $\displaystyle 2\%$ with the observed peak value at $\displaystyle\nu=0.314\ THz$. Though such a modeling is unrealistic in what it does not take into account the details of the protein structure and the associated normal modes suhre , it nonetheless catches a relevant aspect of the global deformation dynamics of the BSA molecule, namely the activation of a collective oscillation, also suggesting that the physical parameters adopted correspond quite well to the situation investigated bosonpeak . Secondary resonances are also present in both spectra. A possible explanation could be tentatively given considering torsional modes. These could be activated at the frequencies given by the relation bastrukov

[TABLE]

where $\displaystyle\nu_{0}$ is given by equation (51), whence, for $\displaystyle l=2$ and $\displaystyle l=3$, one finds $\displaystyle\nu=0.257\ THz$ and $\displaystyle\nu=0.246\ THz$ respectively. These could be associated with the two weaker absorption lines observed at $\displaystyle\nu=0.278\ THz$ and $\displaystyle\nu=0.285\ THz$. Here the larger discrepancy can be attributed to the non-spherical shape of the BSA, what entails the existence of different moments of inertia according to the rotation axis, whereas the breathing mode is insensitive to this fact. Minor peaks are observed at higher frequencies with the rectenna (falling outside the accessible frequency range of the near-field probe) when the protein solution is illuminated with blue laser light. However, the blue light illumination could produce spurious signals from the 2D electron gas of the FET junction as a consequence of electron-hole pairs excitation causing a change of the transistor channel conductivity. This effect is well known and studied in the literature bubanga so that minor peaks could be instrumental artifacts due to this electron-hole pair creation effect.

Let us stress an important point: computational normal mode analysis for proteins has shown nearly continuous vibrational density of states suhre which have also been proved nearly uniformly optically active. Moreover, the coupling of these vibrational modes with water results in broad absorption features heyden ; sun . But this is true at thermal equilibrium, whereas under out-of-equilibrium phonon condensation the energy content of all the normal modes is strongly depleted with the exception of the collective mode. Whence a narrow absorption feature. And, in fact, the computation of the function $\displaystyle L(\omega)$ in Section IV shows how a dipole actively oscillating at a given frequency entails an absorption feature of shape similar to the experimentally observed one.

According to our classical version of the Fröhlich model, it is also expected that the appearance of a collective oscillation should exhibit a threshold-like behaviour when increasing the energy flowing through the protein. Actually, Figure 5 (a) presents a clear threshold in the intensity of the resonance peak at $\displaystyle 0.314$ THz when the optical input power exceeds 10 $\displaystyle\mu$W. By using a classical formalism for the analysis of the out-of-equilibrium phonon condensation we have calculated the intensity of the normal vibrational modes of the BSA-protein as a function of the source power injected through the protein. Figure 5 (b) highlights a threshold-like behaviour of the intensity of the fundamental mode that accumulates the energy at the expenses of the excited modes, in qualitative agreement with the experimental outcome. By increasing the number of modes included in the calculation this threshold becomes more and more evident. The experimental and theoretical results reported in Figure 5 agree also in displaying a saturation effect occurring at large values of the energy input rate.

The observed spectra are certainly due to the light-excited protein because the spectral feature at $\displaystyle 0.314$ THz was not observed by illuminating: i) water alone; ii) a aqueous solution of the fluorochrome in the absence of the protein; iii) a aqueous solution of the BSA protein without the bound fluorochrome (see Section IV). On the other hand, the observed spectral line at $\displaystyle 0.314$ THz immediately disappears by switching off the laser. Remarkably, the spectra obtained with two independent and different experimental setups, based on two different methods of detection of the THz radiation, operated in two different laboratories, show a strikingly good overlap of the respective absorption line profiles at $\displaystyle 0.314$ THz. A result in excellent agreement with the frequency of $\displaystyle 0.308$ THz predicted for the spheroidal (collective) vibrational mode computed on the basis of recent experimental measurements of the relevant parameters of the BSA protein. This triple concordance among the results so far obtained should be enough to exclude experimental artifacts. This notwithstanding various sources of artifact for the observed phenomenology were considered. A first objection suggests itself, namely that the observed phenomenology is just a trivial heating effect due to the laser light. This would be true in the absence of a frequency dependent response to the injection of energy into the proteins. Heating indeed means increasing of the kinetic energy of the atoms and group of atoms of the protein entailing energy equipartition among the vibrational modes of the protein. Hence neither a frequency dependent effect nor a threshold effect for the energy input rate would have been measured, but just a new thermal equilibrium state would have been achieved. To the contrary, each protein - submitted to continuous energy feeding and energy dissipation - behaves as an open system undergoing a non-equilibrium phase transition: when the ratio between energy gain and losses exceeds a critical value a collective behavior sets in producing the phonon condensation, as it has been already discussed in Section II. Then another question arises about the conversion mechanism of the visible light energy absorbed by the electrons of the complex protein/dyes into the vibrational modes of the proteins. This mechanism of rapid intramolecular dissipation of energy through quake-like structural motions as a consequence of a perturbation (such as the breaking of a chemical bond or the absorption of photons through electronic transitions) is being given increasing experimental attention ansari ; levantino ; arnlund ; brinkmann and is referred to as “protein quake”. Similarly to an earthquake, this effect describes how a protein strain is released at a focus or “hot-point” (in our case the fluorochromes) and then rapidly spreads as a structural deformation through waves, thus exciting protein vibrational modes. Another source of artifact could be the apparition of standing waves and related interferences, but these would have been easily identified. Furthermore, there would be no reason for such interferences - if any - to manifest themselves as a consequence of the excitation of fluorochromes through the blue light [see Figure (6b)].

IV Materials and methods

IV.1 Sample preparation

Bovine Serum Albumine (BSA) was purchased from Sigma-Aldrich (A7030) (St. Louis, MO) in lyophilised powder, protease free, fatty acid free ( $\displaystyle\leq 0.02\%$), essentially globulin free; being the degree of purity of the BSA used higher than $\displaystyle 98\%$, it was used without further purification and was dissolved directly into bidistilled water. The BSA molecules have been labeled with the fluorochrome AF488 5-TFP (Alexa Fluor $\displaystyle 488$ Carboxylic Acid, $\displaystyle 2,3,5,6-$Tetrafluorophenyl Ester), $\displaystyle 5$-isomer (A$\displaystyle 30005$) which was purchased from Molecular Probes Invitrogen.

The dye has excitation/emission of 495/515 nm and a molar extinction coefficient of $\displaystyle\varepsilon_{495}=71000$ M cm*-1*. Both the protein and dye concentrations have been determined measuring their absorbance with a Nanodrop 1000 Spectrophotometer (ThermoScientific), at $\displaystyle 280$ nm with a molar extinction coefficient of $\displaystyle 36000$ M cm*-1* and at 495 nm with a molar extinction coefficient of $\displaystyle 71000$ M cm*-1*, respectively. The chemically labeled BSA molecules have been obtained by homemade labelling using 2 mg/mL aqueous solution concentration of proteins with an initial ratio of concentrations [A488 5-TFP]/[BSA] equal to 8 in sodium bicarbonate at a pH $\displaystyle\approx$ 8.5, during one hour. Unconjugated dye was removed using a PD-mini-trap G25 (GE Health Care) according to the instructions of the manufacturer, using gravity protocol, and the degree of labeling was determined spectroscopically. After purification each BSA molecule has between 5 and 6 fluorochromes attached.

IV.2 Experimental detection of the collective mode

Two separated and different experiments, performed in Montpellier (France) and in Rome (Italy), respectively, have been performed to get the terahertz non-equilibrium spectra of the model protein chosen. The former used a microwire as local probe whereas the latter used a nanorectenna.

IV.2.1 Microwire-based THz spectroscopy

A constitutive element of the experimental setup is a tunable and continuous-wave primary source from Virginia Diodes Inc.(https://www.vadiodes.com/en/) thus emitting in the 0.22-0.33 THz frequency range with an average power of 1 mW. The high spectral resolution ($\displaystyle<300$ Hz) of the continuous-wave source allows an accurate detection of resonances. The emitted radiation beam is focused on the samples of protein solutions on which, at option, a 488 nm light beam – produced by an Argon laser – can also be focused. The blue light provides the proteins with the necessary energy to activate a collective vibrational state. The latter being an out-of-equilibrium state because it is kept by a non-thermal energy supply.

The THz near-field scanning spectroscopy technique in aqueous medium is performed by resorting to a homemade micro-coaxial (i.e. subwavelength) near-field probe put inside a metallic rectangular waveguide connected to an heterodyne head and an electrical spectrum analyser. The subwavelength diameter of the wire allows an extremely focused enhancement of the longitudinal component of the electric field at its end over a volume of 4 pL. In 1995, F. Keilmann keilmann highlighted the advantages of introducing a metal wire in a circular metal waveguide to produce probes for near-field microscopy in the far–infrared and microwaves frequency domains. The advantage of this method is to avoid the frequency cut–off when the diameter of the guide is a subwavelength one. The circular waveguide is transformed into a coaxial waveguide which does not have low frequency cut-off and which makes superfocalisation and high–resolution imaging possible.

However, for the experiments reported in the present work, a rectangular waveguide was used, and this entails some different phenomena since the micro-wire is soldered along the long axis of the waveguide and bent to exit it. This has two main consequences: on the one hand, the micro-wire enables a modal transition and, on the other hand, it serves as a waveguide. The bent portion of the micro-wire allows the conversion of the fundamental mode $\displaystyle TE_{01}$ inside the rectangular waveguide into a $\displaystyle TM_{01}$ mode along the wire. To optimize the coupling efficiency between the near-field and the probe, a special care has been first paid to the positioning of the micro-wire inside the guide. More precisely, the maximum of the near field signal is theoretically attained for a wire positioned at $\displaystyle L_{inside}=p\frac{\lambda_{TE_{01}}}{4}$ from the open side of the waveguide, where $\displaystyle\lambda_{TE_{01}}$ is the wavelength of the $\displaystyle TE_{01}$ mode at the considered frequency. At 0.3 THz, the micro-wire must be fixed at an entire multiple $\displaystyle p$ of 250 $\displaystyle\mu$m. The best compromise between efficiency and technical possibilities was found for $\displaystyle p=4$ that gives $\displaystyle L_{inside}=1$ mm. The second parameter to take into account is the angular positioning of the bent portion of the micro-wire inside the rectangular waveguide. The coupling is maximum when the wire is parallel to the orientation of the electric field in the guide, that is, along the long-axis of the rectangular waveguide. Finally we also paid attention to two relevant parameters that are the total length and the diameter of the wire. The intensity of the electric field is roughly sinusoidal and its maximum is attained when the total length $\displaystyle L$ is a multiple of the half-wavelength. The best compromise between theory and technological possibilities has given a total length of $\displaystyle L=2$ mm. Since the electric field intensity at the micro-wire extremity increases when the diameter decreases we used a wire of 12 $\displaystyle\mu$m diameter. All the previously mentioned parameters have also been simulated using CSTmicrowaves studio*®* (https://www.cst.com/products/cstmws) to ensure the better coupling efficiency as possible.

All the measurements have been performed at room temperature. A drop of the protein solution sample is placed under the near-field probe which is directly immersed inside the solution. A typical experiment consists in measuring the near-field electric field intensity through a reference medium (water) and the protein solution sample. A sweep through the frequency domain accessible to the THz source is performed alternatively when the blue light is switched-on and off to yield a difference spectrum showing the absorption features that are attributed to collective vibrations of each protein molecule.

Figure 7 displays the spectra, stabilised in time, obtained by varying the power of the laser illuminating the solution of labelled proteins at the constant concentration of 1 mg/mL. Notice that the secondary peaks appear and develop together with the main peak at 0.314 THz therefore also these secondary peaks are vibrational modes of the entire protein, giving some support to the hypothesis attributing them to torsion modes. However, these kind of modes cannot be accounted for by our theoretical model because it contains no explicit information about the geometry of the molecule.

IV.2.2 Rectenna-based THz spectroscopy

In order to minimize the optical depth of water, in the second experiment (Rome) the probe domain was reduced to a volume of $\displaystyle 10\times 10$ microns in $\displaystyle xy$ (horizontal plane), and to about $\displaystyle 2$ microns in $\displaystyle z$ (vertical axis). To confine the THz radiation (wavelength $\displaystyle\lambda$ around 1 mm) to such a deeply sub-wavelength region, a plasmonic antenna is used giliberti . This device is based on two main components: a planar metal antenna with length close to $\displaystyle\lambda/2$ (bow-tie, broadband type) that produces a high THz field region with antenna feed gap of $\displaystyle 10\times 10$ microns; and a plasma wave FET transistor, which is a THz nonlinear electronic device integrated in the feed gap of the antenna, that provides an electric signal proportional to the THz field strength in the feed gap only. The plasma wave transistor was first introduced by Dyakonov and Shur in 1993 dyakonov and further developed by many authors knap , the main advantage being that standard microwave transistor technology can be employed for fabricating a device which is sensitive to THz radiation. The device is mounted in a package with a silicon lens pressed on the back of the semiconductor substrate and illuminated from below with a tunable THz oscillator (by Virginia Diodes Inc., 0.18 to 0.4 THz) through a set of off-axis parabolic mirrors. The resolution of the free-running oscillator is 2 GHz. After acquiring the empty-channel response spectrum of the device, a drop of protein solution was cast on the top (air) side of the device with a micropipette (a volume of 1 microliter was drop-casted). The drop extends over the entire antenna, but the radiation comes from below (i.e. from inside the substrate) and it is not attenuated by the whole drop. In this experiment absorption from molecules outside the antenna feedgap in the $\displaystyle xy$ directions can be disregarded. Also, in the $\displaystyle z$ direction the field extends for less than 2 microns due to the plasma wave properties (the extension of plasmonic field in the $\displaystyle z$ direction is calculated by Finite Element Modeling simulations). Therefore, the number of molecules probed in this setup is that present in a $\displaystyle 10\times 10\times 2$ micron volume, $\displaystyle 10^{-6}$ times less than the number of molecules present in a $\displaystyle 1\ mm^{3}$ diffraction limited focus.

IV.3 Spectroscopic detection of the collective mode

In both experiments the collective oscillation of the BSA protein is seen as a spectroscopic absorption feature. At variance with standard absorption spectroscopy, where the radiation entering the absorbing medium is responsible for the creation of atomic or molecular excited states, in the experiments reported here the THz radiation is used to detect an already excited state of the molecules. Actually, the collective oscillation of the proteins makes them behave as very small antennas (Hertzian dipoles) with the characteristic property of antennas of simultaneously absorbing and emitting electromagnetic radiation. However, the absorption along the THz optical path cannot be compensated by the radiation emitted by the oscillating dipoles because it spreads over all the directions in space. The net result is an absorption feature. If we denote with $\displaystyle\vec{\mu}(t)$ the dipole moment of a protein and with $\displaystyle\vec{E}(t)=\vec{E}_{0}\cos(\omega t)$ the electric field of the THz radiation, the attenuation of the electric energy density within the drop of protein solution is proportional to the work done by the electric field, that is, $\displaystyle L=-\vec{\mu}(t)\cdot\vec{E}(t)$. The oscillation of the dipole moment is necessarily damped, predominantly because of bremsstrahlung emission, so that, denoting by $\displaystyle\tau$ the lifetime of the activated collective oscillation and by $\displaystyle\omega_{c}$ its frequency, we can set $\displaystyle\vec{\mu}(t)=\vec{\mu}_{0}e^{-t/\tau}\cos(\omega_{c}t)$. Thus, after averaging over all the relative orientations and all the phase differences $\displaystyle\phi$ between $\displaystyle\vec{\mu}(t)$ and $\displaystyle\vec{E}(t)$ such that the electric field does a positive work, we obtain

[TABLE]

This is the elementary attenuation process of the THz radiation. This process is repeated in time for each molecule at a rate proportional to the intensity of the drop illumination with the blue light. Moreover, the total attenuation is proportional to the concentration of absorbing molecules in the protein solution. Equation (52) gives for $\displaystyle L(\omega)$ a Lorentzian shape centered at $\displaystyle\omega_{c}$, the resonance frequency of the collective oscillation of the BSA protein. Figure 10 shows three different shapes of the function $\displaystyle L(\omega)$ obtained for different values of $\displaystyle\tau$ (in arbitrary units). These line shapes show that an already vibrating underdamped dipole absorbes the weak terahertz-radiation-probe with the same frequency pattern of the $\displaystyle 0.314$ THz absorption line reported in Figure 3 of the main text. The latter is well fitted by the $\displaystyle L(\omega)$ function by using a $\displaystyle Q$ factor of $\displaystyle 50$ ($\displaystyle Q=\Delta\nu/\nu$, that is the ratio between the line width and the line frequency).

In Figures 6a and 6b some control spectra are reported which have been obtained with the microwire antenna. The two groups of spectra refer to the blue light emitting laser switched off and on, respectively. The observed absorption line at $\displaystyle 0.314$ THz is clearly found only when the fluorochrome AF488 is bound to the BSA molecules and in presence of $\displaystyle 488$ nm laser light illumination of the protein solution. These spectra rule out any other origin of the observed absorption feature beside the proposed one in the main text of the present work.

In Figure 8 the photoresponse spectra obtained with the rectenna probe highlight the same phenomenology: the absorption line at $\displaystyle 0.314$ THz is present only when the aqueous solution of AF488-labelled BSA molecules is illuminated with $\displaystyle 488$ nm laser light. A clear difference is again observed between the spectra when the blue light illumination is switched off and on. Artefacts possibly due to blue light illuminated water alone, or possibly due to the blue light illuminated AF488 dye in aqueous solution, are ruled out by the control spectra reported in Figure 9, and are in agreement with the analogous controls performed with the microwire antenna reported in Figures 6a and 6b.

IV.4 Experimental activation time scale of the collective mode

In both cases of Montpellier and Rome experiments, the absorption feature, which is attributed to the activation of an out-of-equilibrium collective oscillation of the protein, appeared after several minutes of blue light irradiation: about $\displaystyle 4$ minutes in the Montpellier experiment, after $\displaystyle 10$ minutes in the Rome experiment (though in this latter case no measurement was performed before $\displaystyle 10$ minutes). In what follows we provide a consistency estimate of this long activation time with the physical conditions of the experiments. Going beyond a qualitative estimate is hardly feasible, and in any case would not affect the meaning of the experimental outcomes reported in the present work. Thus we proceed by estimating the balance between the energy input rate and the energy loss rate of each molecule. As we shall discuss at the end of this Section, this long time is not an hindrance to the biological relevance of the phenomenon reported in the present work.

An elementary account of the balance between energy gain and loss for each protein can be given by the equation

[TABLE]

where $\displaystyle E$ in the l.h.s is the numerical value of the energy of the system described by the Hamiltonian (2) of the previous Section, the first term in the r.h.s. represents the radiative losses (bremsstrahlung) of the accelerated dipole of moment $\displaystyle Ze\mathbf{x}$, (where $\displaystyle Z$ is the number of elementary charges $\displaystyle e$ and $\displaystyle\mathbf{x}$ is the separation between the center of positive and negative charges) and $\displaystyle W$ is the energy input rate due to the $\displaystyle 488$ nm light irradiation. $\displaystyle\Gamma$ stands for other kind of energy losses leading to thermalization of the protein with its aqueous environment.

For what concerns the energy input rate $\displaystyle W$, the energy difference between the entering photons, of wavelength $\displaystyle 4.88\times 10^{-5}$ cm, and the fluorescently emitted photons, of wavelength $\displaystyle 5.3\times 10^{-5}$ cm, amounts to $\displaystyle 3.18\times 10^{-13}{\rm erg}=0.19\ {\rm eV}$. When the Argon laser operates at $\displaystyle 500$ $\displaystyle\mu$W the number of emitted photons per second is $\displaystyle 1.2\times 10^{15}$. The cross section of the fluorochrome Alexa488 is $\displaystyle 3.2\times 10^{-16}$ $\displaystyle{\rm cm}^{2}$ (free dye in water) lili , so that assuming a Gaussian beam section (at $\displaystyle 1/e^{2}$) of $\displaystyle 7.8\times 10^{-3}$ $\displaystyle{\rm cm}^{2}$ (corresponding to a diameter of about $\displaystyle 0.1$ $\displaystyle{\rm cm}$), we find that each Alexa488 molecule receives approximately $\displaystyle 50$ photons per second. Each protein molecule has an average number of $\displaystyle 5-6$ fluorochromes bound to it so that, considering that part of them could be partly shadowed by the protein itself, we can assume that the average number of photons received by each protein per second belongs to the interval $\displaystyle 120\div 300$. Hence the upper bound for $\displaystyle W$ is estimated between $\displaystyle 3.8\times 10^{-11}$ erg s*-1* and $\displaystyle 9.5\times 10^{-11}$ erg s*-1*. This estimate is obtained under the assumption that all the photons hitting the Alexa488 molecules are absorbed; the quantum yield of these dye molecules is high (92 $\displaystyle\%$) so that practically all the harvested photons are converted into fluorescence thus providing an upper bound of the energy available for the protein excitation. Of course part of this energy can be dissipated in the sorrounding environment in the form of heat.

For what concerns the radiation losses, we proceed to estimate the time average value of $\displaystyle{\overline{|\ddot{\mathbf{x}}|^{2}}}$. Since the $\displaystyle 63\%$ of the N-isoform of the BSA protein consists of $\displaystyle\alpha$-helices, we can assume that the so-called accordionlike vibrational modes of these helices provide the largest contribution to the protein extension and thus to the variation of the protein’s dipole moment. At thermal equilibrium we can assume that the dipole elongation $\displaystyle{\mathbf{x}}(t)$ can be the result of an incoherent superposition like $\displaystyle{\mathbf{x}}(t)=\sum_{i=1}^{30}A_{i}\cos(\omega_{i}t+\xi_{i})$, where the $\displaystyle\xi_{i}$ are random phases, and the sum runs over the $\displaystyle 30$ $\displaystyle\alpha$-helices of which the protein is composed. To compute the time average of $\displaystyle|\ddot{\mathbf{x}}|^{2}$ we have to estimate the quantity $\displaystyle\sum_{i=1}^{30}\frac{1}{2}{\overline{A_{i}^{2}\omega_{i}^{4}}}$ (since the mixed terms average to zero). Using the formula $\displaystyle\omega_{i}=2\pi(1/2L_{i})(E/\rho)^{1/2}$ for the frequencies of the accordionlike modes of the $\displaystyle\alpha$-helices alphahelix , where $\displaystyle 1/\rho=0.74$, $\displaystyle E=2.31\times 10^{11}$ dyne cm*-2* for $\displaystyle\alpha$-helices, $\displaystyle L_{i}=1.5N$ $\displaystyle\mathring{\text{A}}$, where $\displaystyle N$ is the number of amino acid residues, and assuming that $\displaystyle A_{i}$ belongs to the interval $\displaystyle 0.05L_{i}\div 0.1L_{i}$, and $\displaystyle Ze=-13\times 4.803\times 10^{-10}Fr$ for the BSA protein at neutral pH barbosa , we finally find:

[TABLE]

The range of values of the energy input rate and that of the radiative losses are thus almost overlapping. Making these estimates more precise is hardly feasible and is beyond the aims of the present work. What matters here is that since we experimentally observe the activation of the collective mode of the BSA molecules, the energy input rate $\displaystyle W$ must exceed the rate of all the losses, radiative and non, and if this happens only by a small amount, then $\displaystyle{dE}/{dt}$ can be so small as to require some minutes in order to accumulate enough energy in each molecule [this can be thought of as a steady increase of the value of the third term of the unperturbed part of the Hamiltonian in Eq.(2)]. When the energy so accumulated exceeds a critical threshold value, the system undergoes a condensation transition channelling the largest fraction of the input energy into the lowest frequency mode(s). An elementary estimate of the rate of energy losses due to the bremsstrahlung in the condensed phase is obtained by entering into the Larmor formula (54) the dipole acceleration $\displaystyle|\ddot{\mathbf{x}}|$ computed from $\displaystyle{\mathbf{x}}(t)=A\cos(\omega_{C}t)$, putting $\displaystyle\omega_{C}=0.314$ THz, which is the frequency of the collective mode, and assuming that the dipole elongation $\displaystyle A$ is about $\displaystyle 10$ $\displaystyle\mathring{\text{A}}$. The breakeven between the energy input and the radiative losses is then found for a dipole moment in the range $\displaystyle 14500\div 23000$ Debye, corresponding to an effective number of charges $\displaystyle Z$ approximately in the range $\displaystyle 290\div 460$. In other words, in the condensed phase the total dipole of the protein oscillates at a low frequency (with respect to the accordion-like modes of the $\displaystyle\alpha$-helices) so that the collective oscillation can be stable with respect to radiative losses up to the activation of large values of the protein dipole moment.

We have neglected any estimate of the collisional losses - represented by $\displaystyle\Gamma$ in equation (53) - because their estimate as viscous losses - on the basis of Stokes formula - gives an unreasonably high energy dissipation rate with respect to the energy input rate. If this were the case, no collective vibration would be observed at all. There is a vast literature about relaxation phenomena, thus thermalization, of proteins subject to different kinds of radiative excitations. The typical relaxation time scales are in the range of picoseconds to nanoseconds, with a remarkable exception - recently reported in lundholm - where the thermalization of terahertz photons in a protein crystal happens on a micro- to milli-second time scale, which is interpreted as being due to Fröhlich condensation (no matter whether quantum or classical). To the contrary, typical picoseconds to nanoseconds relaxations usually pertain side chains relaxations or small groups of atoms. Also in the case of relaxations of collective modes detected with FIR spectroscopy or with THz Time Domain Spectroscopy xie ; kneller the physical conditions are very different with respect to our present experiments. In the mentioned literature the collective modes of proteins are probed at thermal equilibrium, thus in presence of all the other vibrational modes and under the action of weak exciting fields. In our experiments the protein molecules are out-of-thermal equilibrium and are strongly excited through an internal cascade of mode-mode couplings stemming from the fluorescence decay of the attached fluorochromes, and the THz radiation is only used to read the presence of the collective oscillation (see the preceding Subsection).

On the other hand, it has been recently found that the hydration shell of the BSA protein is $\displaystyle 25$ $\displaystyle\mathring{\text{A}}$ thick susko , and this seem to be a generic property of solvated proteins ebbinghaus . The microrheology waigh of this kind of water-protein system, and in particular its high frequency viscoelasticity, is still an open research subject, making a quantitative estimate of the term $\displaystyle\Gamma$ in equation (53) hardly feasible.

A comment about the prospective biological relevance of the observed phenomenology is in order. The main energy source within living cells is provided by ATP hydrolysis. The typical intracellular concentration of ATP molecules is given around $\displaystyle 1\ \text{mM}$ implying that a protein molecule in the cell undergoes around $\displaystyle 10^{6}$ collisions with ATP molecules per second alberts . Given the standard free-energy obtained from ATP hydrolysis estimated around $\displaystyle 50\ \text{kJ.mol}^{\text{-1}}=8.306\times 10^{-13}\text{erg}$, if we assume that only $\displaystyle 1\%$ or $\displaystyle 2\%$ of the collisions with ATP will provide energy, a power supply between $\displaystyle 8.306\times 10^{-9}\text{erg s}^{\text{-1}}$ and $\displaystyle 1.6\times 10^{-8}\text{erg s}^{\text{-1}}$ is potentially available. This could be at least two orders of magnitude larger than the power supplied to each protein in our experiments, but reasonably even much more than two orders of magnitude because we have assumed a hundred percent conversion efficiency of the energy supplied by the laser into internal vibrations of the protein. But this can be hardly the case, thus the condensation mechanism in vivo can be considerably faster.

V Concluding remarks

In the present paper we have studied a classical version of a quantum model put forward many years ago by H. Fröhlich Frohlich:1968 ; Froeh:1970 ; Frohlich:1977 . The classical model displays the same phenomenon of a Bose-like phonon condensation of the normal vibrational modes of a macromolecule that was predicted by the original quantum model. Even though our classical model (just like its quantum predecessor) is very idealized, that is, even if our model cannot be quantitatively predictive, it can nevertheless qualitatively catch the presence of a robust and generic phenomenon undergone by open systems with internal nonlinear interactions: a non-equilibrium phase transition when a control parameter exceeds a critical threshold value. This kind of phase transition is found to correspond to a channelling into the lowest frequency mode of almost all the energy pumped into the system, that is to a collective global molecular vibration. This theoretical result has motivated an experimental - and successful - attempt at confirming it. So in the second part of the present paper we have reported - for a model protein - the unprecedented observation of an out-of-thermal-equilibrium collective oscillation which is in qualitative agreement with the theoretical model and in excellent quantitative agreement with the theoretically expected value of a spheroidal collective vibration of the whole protein.

The novelty is not the collective oscillation in itself, because several terahertz spectroscopic studies have reported about collective modes of proteins, but all of these studies were performed at thermal equilibrium and mainly carried on using dry or low-hydrated powders because of the very strong absorption of water turton ; acbas ; falconer ; markelz , even though more recent studies also addressed solvated proteins xu ; ebbinghaus . The novelty of both our theoretical and experimental contributions, let us stress this point again, consists of considering out of thermal equilibrium conditions. On the other hand, recent studies on solvated BSA in THz bye and sub THz frequency range have shown susko broad resonances due to an efficient coupling of low frequency modes of the protein with the surrounding water, and though all of these works were performed at thermal equilibrium, in common with these previous studies the experimental part of our present work confirms the relevance of the coupling of the protein with the surrounding water molecules. In fact, the strong absorption feature that we observed in a aqueous solution of the BSA protein put out of thermal equilibrium, reveals that the protein vibrating in its collective mode has to be dressed by ordered layers of water molecules in order to attain an effective dipole moment sufficiently large to overcome the strong absorption of bulk water.

We anticipate that the theoretical and experimental sides of the work presented in this paper could open a broad domain of systematic investigations about out-of-equilibrium activation mechanisms and properties of collective oscillations of different kinds of biomolecules. Furthermore, as already mentioned in the Introduction, the possibility of exciting out-of-thermal-equilibrium collective oscillations of macromolecules is specifically interesting as a necessary condition to activate resonant long-distance electrodynamic intermolecular interactions Preto:2015 . Thus our results explain why electrodynamic interactions between biomolecules have hitherto eluded detection, in fact no attempt has ever been done to detect them by involving biomolecules vibrating out-of-equilibrium. Consequently, our work also motivates new efforts to detect these electrodynamic intermolecular interactions nardecchia1 ; nardecchia2 .

Acknowledgments

We warmly thank Jack A. Tuszynski, Michal Cifra and Anirban Bandyopadhyay for useful discussions ; A. Penarier, T. Cohen and F. Cano for their help in the simulation and the realization of the micro-wire local probes. The project leading to this publication has received funding from the Excellence Initiative of Aix-Marseille University - A*Midex, a French “Investissements d’Avenir” programme. This work was also financially supported by the Seventh Framework Programme for Research of the European Commission under FET-Proactive grant TOPDRIM (FP7-ICT-318121), by the projects SIDERANT and NEBULA financed by the french CNRS, by the Occitanie Region and by Montpellier University through its TOP platform and NUMEV LabeX.

Author contributions I.N., M.L. and J.T. performed the experiments in Montpellier, J.T. designed and built the experimental setup in Montpellier with the support of P.N., V.G. and M.O. performed the experiments in Rome, the classical formulation of Fröhlich condensation stems from the PhD thesis of M.G., J.P. has also contributed the theoretical part, I.D. did some numerical work, J.S. contributed on the biophysical and biochemical aspects of the work, J.T. and L.V. contributed to conceive the Montpellier experiments and participated to the supervision of all the development of the project. M.P., in quality of project leader supervised and intervened in all the aspects the project. All the authors contributed to the discussion and to the analysis of the results. M.P. wrote the paper with the help of M.G., M.O., J.S., J.T., L.V., and J.T. made all the figures.