Topology, Landscapes, and Biomolecular Energy Transport

TL;DR

This paper explores how the topology and nonlinear dynamics of biomolecular energy landscapes influence transient thermal transport, revealing insights into energy migration and biomolecular function through advanced simulations and analysis.

Contribution

It introduces a novel approach linking energy landscape topology to thermal transport in biomolecules, enabling detailed analysis of energy migration mechanisms.

Findings

Transport pathways are determined by molecular contact topology.

Nonlinear processes dominate over coherent transport at short scales.

Vibrational energy transport probes biomolecular dynamics.

Abstract

While ubiquitous, energy redistribution remains a poorly understood facet of the nonequilibrium thermodynamics of biomolecules. At the molecular level, finite-size effects, pronounced nonlinearities, and ballistic processes produce behavior that diverges from the macroscale. Here, we show that transient thermal transport reflects macromolecular energy landscape architecture through the topological characteristics of molecular contacts and the nonlinear processes that mediate dynamics. While the former determines transport pathways via pairwise interactions, the latter reflects frustration within the landscape for local conformational rearrangements. Unlike transport through small-molecule systems, such as alkanes, nonlinearity dominates over coherent processes at even quite short time- and length-scales. Our exhaustive all-atom simulations and novel local-in-time and space analysis, applicable to both theory and experiment, permit dissection of energy migration in biomolecules. The approach demonstrates that vibrational energy transport can probe otherwise inaccessible aspects of macromolecular dynamics and interactions that underly biological function.

Figures

Equations

$$\begin{gathered} {\dot{E}}_j(t)=\sum _i[k_{ij}(t)E_i(t)-k_{ji}(t)E_j(t)] \\ -k_{\text{s},j}(t)[E_{\text{s},j}(t)-E_j(t)]. \end{gathered}$$

$$\mathbf{k}(t_n)=\underset{\mathbf{k}\ge 0}{\min }\frac{1}{2}||\mathbf{G}(t_n)\cdot \mathbf{k}-\mathbf{d}(t_n)|{|}^2$$

Full text

Topology, Landscapes, and Biomolecular Energy Transport

Justin E. Elenewski

Biophysics Group, Microsystems and Nanotechnology Division, Physical Measurement Laboratory, National Institute of Standards and Technology, Gaithersburg, MD 20899, USA

Maryland Nanocenter, University of Maryland, College Park, MD 20742, USA

  

Kirill A. Velizhanin

Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545 USA

  

Michael Zwolak

mpz@nist.gov

Biophysics Group, Microsystems and Nanotechnology Division, Physical Measurement Laboratory, National Institute of Standards and Technology, Gaithersburg, MD 20899, USA

Abstract

While ubiquitous, energy redistribution remains a poorly understood facet of the nonequilibrium thermodynamics of biomolecules. At the molecular level, finite–size effects, pronounced nonlinearities, and ballistic processes produce behavior that diverges from the macroscale. Here, we show that transient thermal transport reflects macromolecular energy landscape architecture through the topological characteristics of molecular contacts and the nonlinear processes that mediate dynamics. While the former determines transport pathways via pairwise interactions, the latter reflects frustration within the landscape for local conformational rearrangements. Unlike transport through small–molecule systems, such as alkanes, nonlinearity dominates over coherent processes at even quite short time– and length–scales. Our exhaustive all–atom simulations and novel local–in–time and space analysis, applicable to both theory and experiment, permit dissection of energy migration in biomolecules. The approach demonstrates that vibrational energy transport can probe otherwise inaccessible aspects of macromolecular dynamics and interactions that underly biological function.

Introduction

Biological systems are characterized by a persistent nonequilibrium state, maintained by the open metabolic reactions that drive self–replication. Directed redistribution of energy is an intrinsic feature, serving to generate mechanical motion Andrieux2006 ; Hwang2017 , mediate allosteric communication Tu2008 ; Wang2017 ; Buchenberg2017 , and drive bioenergetic processes Ansari1985 ; Nedergaard2005 ; Reidel2015 . The physical scales of these processes can be surprising: Common enzymatic reactions liberate up to 2 eV of heat repeatedly over micro– to milli–second catalytic cycles Reidel2015 . This energy is redistributed throughout the surrounding protein scaffold within picoseconds and is either dissipated to mitigate thermally–induced stress, leveraged to induce mechanical motion, or employed to promote further catalytic activity. Irrespective of the endpoint, efficient and directed energy transport is critical to the function of these nanoscale machines.

At the macroscale, Fourier’s law, $J=-\kappa\nabla T$ and its time–dependent version capture diffusive heat flow, given by the flux $J$, in response to a temperature gradient $\nabla T$. Those two quantities are related by the thermal conductivity $\kappa$ (or the diffusivity $D$), which can be anisotropic. This situation is more complicated at the nanoscale, where competing ballistic and diffusive transport pathways impede a universal description Cahill2003 ; Cahill2013 . In this context, ballistic wavepackets propagate at the speed of sound in a given vibrational band, up the vibrational mean free path, even without the local thermal gradients required for diffusive transport.

Despite the ubiquity of energy redistribution and flow in biomolecular systems, experiments are difficult Ansari1985 ; Srajer1996 ; Botan2007 ; Helbing2012 ; Barends2015 ; Levantino2015 . In a pioneering work, Botan et al. Botan2007 developed an approach to monitor real–time heat migration in a polypeptide of 2–aminoisobutyric acid (Aib). The approach employs a photoexcitable azobenzene tag as a heater and backbone carbonyl modes as local vibrational thermometers. The results are complex, suggesting a ‘dynamical transition’ temperature above which transport is enhanced Backus2008 ; Backus2008b ; Schade2009 ; Backus2009 . Quantum and nonequilibrium molecular dynamics (NEMD) simulations support the presence of a transition in transport properties, and also suggest that a classical description is realistic Nguyen2010 ; Kobus2010 ; Kobus2011 ; Goj2011 (unlike for small molecules Wang2007 ; Rubtsova2015 ; Qasim2019 ; Rubtsov2019 ; Liu2019 ). However, both the nature of the transition and mechanism of transport remain unclear, with theory giving conflicting accounts Botan2007 ; Backus2008 ; Nguyen2010 ; Kobus2010 ; Schade2012 .

In this work, we utilize molecular dynamics simulations and a new space– and time–local analysis method to explore energy propagation in a paradigmatic polypeptide. We find that Fourier behavior captures the bulk of transient energy flow, provided that one accounts for the fact that fluxes and diffusivities are temperature dependent. Departures from a simple realization of Fourier’s law happen at large temperature gradients, beyond about 15 K/residue, even though transport is still diffusive. The identification of these regimes is not possible through all-atom molecular dynamics alone Nguyen2009 ; Nguyen2010 ; Kobus2010 ; Kobus2011 ; Goj2011 ; Brinkman2016 ; Buchenberg2016 or normal-mode analysis (even when treating anharmonicity as a correction) Leitner2001 ; Yu2003 ; Yu2003b ; Yu2005 ; Leitner2009 ; Gnanasekaran2011 ; Leitner2015 . The former does not unravel the atomic–scale mechanisms of transport and the latter reflects dynamics only at potential energy minima Yu2005 . Within this context, we further demonstrate how the graph–theoretic topology of molecular contacts can define directed pathways for molecular energy redistribution.

Results

Topology and energy propagation pathways. We initiate our investigations using a series of replica–exchange molecular dynamics (REMD) simulations, as the lack of symmetries, granularity, and high-dimensional free energy landscapes of biomolecules necessitate an exhaustive exploration of conformational space Wales1998 ; Wales2006 ; Wales2015 . Our simulation system is a ten–residue Aib helix (Aib10) solvated by chloroform, similar to experimental efforts Botan2007 ; Backus2008 ; Backus2008b ; Schade2009 ; Backus2009 . We previously generated temperature–dependent free energy landscapes for Aib10 at high resolution with replica–exchange simulations Elenewski2018 . From the resulting conformational ensemble, we extract 4000 conformers for each environmental (bath) temperature $T_{\text{B}}$ according to a Boltzmann distribution. This includes structures from both left– and right–handed folding funnels, ensuring a uniform distribution of configurations (Fig. 1a,b). We initiate NEMD simulations in a manner that mimics photoexcitation, distributing $\approx 1.6$ eV of energy between designated vibrational degrees of freedom in each conformer. This is achieved by thermostatting the C–terminal residue to a temperature $T^{\prime}=T_{\text{B}}+\Delta T$, with $\Delta T=670$ K, while holding the remainder of the system at $T_{\text{B}}$. The simultaneous heating of all vibrational degrees of freedom in the heater residue is well–founded, as it yields thermal transport profiles that are indistinguishable from mode–selective heating Botan2007 ; Nguyen2010 . This excess energy then propagates freely within the microcanonical ensemble (i.e., without thermostatting).

The conformational ensemble of Aib10 comprises three general structural motifs (Fig. 1b) corresponding to (i) $3_{10}$–/$\alpha$–helical conformers ($\approx 45$ % of ensemble) with hydrogen bonding between residue $j$ and residue $j+3$ or $j+4$, respectively; (ii) hairpin–like configurations, with hydrogen bonds between the first and last residues of Aib10 ($\approx 15$%); and (iii) unstructured or extended conformers that have no consistent hydrogen bonding ($\approx 40$%) Elenewski2018 . We index these subensembles with $\lambda$. This partition is defined by the underlying free energy landscape, and is thus independent of our thermal transport simulations Elenewski2018 .

In Fig. 1c,d, we present transport profiles for Aib10 versus the ensemble–averaged temperature elevation $\Delta T_{\text{B},j}(t)=\langle T_{j}(t)\rangle-T_{\text{B}}$ of the $j^{\mathrm{th}}$ residue, or $\Delta T_{\text{B},j}^{\lambda}-\Delta T_{\text{B},j}$ for subensemble $\lambda$. The full-ensemble profile exhibits a weak thermal front that traverses the peptide within $2$ ps, which is also apparent in the helical ensemble (Fig. 1d). This corresponds to backbone propagation at $v=1.7$ nm ps*-1*, approaching ballistic transport velocities in biomolecular materials and alkyl chains Botan2007 ; Yue2015 ; Rubtsova2015 ; Qasim2019 ; Rubtsov2019 . While this channel is weak, additional ballistic pathways may exist at lower group velocities in different vibrational bands Yue2015 ; Qasim2016 , though these will inevitably be obscured by more prominent diffusive features. There is also rapid transport with both ballistic and diffusive characteristics across hydrogen-bonded regions, which can be seen in the helical and hairpin conformers (see discussion below).

While a ballistic pathway exists, the majority of energy transport is nonetheless diffusive — yielding a broad profile that is sensitive to both temperature and molecular conformation. We separate diffusive and ballistic behavior by coarse–graining in time (into 100 fs bins), averaging away signatures of very fast dynamics, but retain spatial coarse–graining into individual amino acid residues. We will develop time–dependent quantitative methods to extract diffusivities, free energies, and other characteristics from temperature–based data. However, to facilitate comparison with prior theory and experiment, we initially calculate diffusivities via the time to reach the maximal temperature for each residue. Considering just the helical subensemble for fitting, the temperature-dependent thermal diffusivity $D(T_{\text{B}})$ has distinct low– and high–temperature regimes (Fig. 2a), which are also reflected in the net heat transfer (Fig. 2b). This qualitative behavior agrees with experimental Botan2007 ; Backus2008b ; Backus2009 and theoretical Botan2007 ; Nguyen2010 ; Kobus2010 efforts. These, though, report diffusivities of 0.02 nm2 ps*-1* and 0.1 nm2 ps*-1*, respectively. Theoretical $D(T_{\text{B}})$ from this type of estimate consistently exceed experimental values for Aib10 but are comparable to bulk materials Rubtsov2019 and other proteins Leitner2015 . Force-field parameterization likely contributes to this discrepancy in part. We will see, through an alternate analysis, that residual ballistic components also play a role. The crossover near 270 K is consistent with prior efforts, which ascribe this behavior to a glass–like dynamical transition Botan2007 ; Backus2008b ; Backus2009 ; Nguyen2010 . We will return to this point.

Given this diverse ensemble, it is natural to ask how transport behaves in different conformers. This question was not addressed by prior computational efforts, as they remained below the timescale for structural interconversion in forming their ensemble, sampling only helical configurations and thus a fixed secondary connectivity Botan2007 ; Nguyen2010 . Figure 1d shows the transport profile of the full Aib10 ensemble compared to ensembles that contain only helical motifs, hairpin motifs, or randomly oriented conformers without fixed secondary structure. On a residue–by-residue basis, helical conformers propagate heat more readily than the full ensemble. This is evidenced by less energy retention at the heater site for $t\leq 25$ ps, commensurate with enhanced transfer to its hydrogen–bonded contacts at early times (mostly site 4 for the helix). The randomly oriented conformers transport heat less efficiently, underscored by enhanced energy localization at the first three residues for short times and, later, a rate of energy migration that lies slightly below the full ensemble. We expect a dominant backbone contribution in this case, as longer range contacts are sporadic. Hairpin configurations are intermediate with enhanced transport to certain hydrogen bond contacts (site 10), in turn reducing the amount of heat transport through others (to the fourth site). It should be noted that, while hydrogen bonding can lead to more efficient heat transport for certain conformers, backbone channels always carry the majority of heat. Changes in energy migration are not due to local solvent heating, as the mean temperature of the first two solvation shells increases by at most 5 K over the entire simulation. While the overall cooling rate involves an interplay between heat diffusivity and surface area–dependent solvent coupling, these effects are minor for the systems considered herein (see the Supplementary Discussion).

These observations indicate that topologically nontrivial configurations yield efficient pathways for vibrational energy migration. The importance of secondary and tertiary contacts has been previously invoked when describing transport within a single conformer of HP36 Leitner2015 ; Buchenberg2016 . We extend this observation, demonstrating that representative heat transport characteristics can be obtained only when the conformational landscape is comprehensively sampled. This is particularly important for metrologies, where insufficient sampling can lead to erroneous diffusivities and the misidentification of transport pathways. Moreover, changing conditions (temperature, pH, presence of denaturants, etc.) can shift the conformational ensemble, particularly near structural transitions. This will be detected by the energy transport, including the capture of additional information about underlying interactions Velizhanin2011 ; Chien2013 ; velizhanin_crossover_2015 .

Heat fluxes and energy landscape topography. While molecular connectivity clearly determines transport pathways, NEMD simulations and existing analysis frameworks afford no immediate means to reconcile temperature–dependent features with microscopic processes and the underlying free energy landscape. To directly address this, we analyze the intermediate–timescale dynamics of NEMD trajectories – restricting to helical Aib10 conformers for both structural heterogeneity and consistency with prior work – using a master equation for the kinetic energy $E_{j}$ of the $j^{\mathrm{th}}$ residue in the peptide:

[TABLE]

In this case, $k_{ij}(t)$ is a rate constant for energy transfer from residue $i$ to residue $j$ and $k_{ji}(t)$ is a distinct rate for the reverse process (see Methods), $k_{\text{s},j}(t)$ is the rate of heat transfer to the solvent bath, and $E_{\text{s},j}(t)$ is the kinetic energy density of the solvent surrounding the $j^{\mathrm{th}}$ residue (scaled to match the residue degrees of freedom). We diverge from earlier work by treating the $k_{ij}(t)$ as parameters that depend on both position and time — thereby implying a temperature dependence. This accommodation is key to our subsequent analysis. Given this arrangement, one can identify two distinct intra–peptide couplings: (i) direct transfer between nearest–neighbors in the peptide backbone $(k_{j,j+1}$ and $k_{j,j-1})$ and (ii) a long distance coupling between hydrogen bonding partners ($k_{j,j+3}$, $k_{j,j+4}$ for ideal $3_{10}$– and $\alpha$–helices, respectively). With additional approximations, the system in Equation (1) becomes well–posed and solvable at all times (see Methods). This diverges from existing master equation analyses which assume rate constants that are time– and space–independent, and thus independent of the local temperatures and gradients Buchenberg2016 . These prior works nonetheless treat a broad network of nonlocal contacts, which combined with the analysis here would constitute a logical extension of our methods.

Our remaining discussion is driven by the pairwise heat fluxes $J_{i,j}(t_{n})=-k_{i,j}(t_{n})[E_{i}(t_{n})-(f_{j}/f_{i})E_{j}(t_{n})]$ and rate constants between coupled residues. Here $f_{j}$ is the number of degrees of freedom for residue $j$ and $t_{n}$ indexes the time domain coarse–graining of the simulation trajectory into $n\leq N$ bins via block averaging. This approach is a finite difference decomposition of the diffusion equation $\dot{E}(x,t)=D\,\nabla^{2}E(x,t)$ at the timescale $\Delta t=t_{n+1}-t_{n}$ and a length-scale $\Delta x$ defined by the distance between adjacent residues. The fluxes come from the finite difference decomposition of $J(x,t)=-D\,\nabla E(x,t)$.

The rate constants $k_{i,j}(t_{n})=D(t_{n})/(\Delta x)^{2}$, in particular, capture biomolecular heat diffusivity $D(t_{n})$ while giving a metric for energy landscape features. We are interested in the distribution of barriers between low–lying conformational minima, specifically those connected by the energy–transmitting structural displacements that are associated with vibrational energy propagation. This latter property is reflected by the local, activated conformational changes underlying transport $k_{i,j}=\Omega_{i,j}\,\exp[-\Delta G_{i,j}/k_{\text{B}}T]$, where $\Delta G_{i,j}$ is the free energy barrier between heat–accepting microstates and $(\Omega_{i,j})^{-1}$ is an effective timescale for free diffusion, influenced by both the protein and its environmental coupling. While each pair of microstates is characterized by a distinct $\Delta G_{i,j}$, these values evolve during heat transport — commensurate with changes in the free energy landscape.

We employ this kinetic approach with an intermediate timescale ($\Delta t=100$ fs), long enough to average over most coherent motion but short enough not to obscure the evolution of energy in time. The distribution of backbone fluxes $J_{\text{BB}}$ is parameterized by an effective temperature gradient $\Delta_{ij}T_{\text{eff}}=2[E_{i}-(f_{j}/f_{i})E_{j}]/3Nk_{\text{B}}$ between residues $i$ and $j$, where the flux is incident on a residue containing $N$ atoms. While transport is explicitly quantified through $J_{\text{BB}}$ for simplicity, the effect of hydrogen bonding is present when fitting the backbone flux distribution at hydrogen bonding sites. The results for $J_{\text{BB}}$ are presented in Fig. 3a. A complimentary analysis for $J_{\text{HB}}$ and a validation of fitting methods are presented in Supplementary Figures 10–15.

Region A. The forward flux $J_{\text{BB}}$ has a linear region for small $\Delta_{ij}T_{\text{eff}}$ (less than about 15 K), although it does not go to zero at $\Delta_{ij}T_{\text{eff}}=0$. Purely diffusive transport will not afford a heat flux in the absence of a local temperature gradient. Thus, a finite $J_{\text{BB}}$ at $\Delta_{ij}T_{\text{eff}}=0$ is a signature of ballistic/coherent behavior. Supporting this interpretation, we find that the zero–gradient flux to decrease with increasing $\Delta t$ during coarse–graining, while only exhibiting small error bands at all scales (thus it is not due to short–timescale fluctuations). A linear fit to this regime gives an effective diffusivity of $D_{\text{eff},\text{A}}=2.3\times 10^{-2}$ nm2 ps*-1* (or conductivity $\kappa_{\text{eff},\text{A}}=3.9\times 10^{-3}$ eV K*-1* ps*-1*). Fitting for small $\Delta_{ij}T_{\text{eff}}$, while ignoring the residual ballistic contribution right around $\Delta_{ij}T_{\text{eff}}=0$, removes high rate constant artifacts. Encouragingly, the magnitude of the resulting diffusivity is consistent with experimental values Botan2007 ; Backus2008 . Employing the time to reach the maximum temperature, as done in prior theoretical work (see discussion above), affords much higher diffusivities. This linear regime has the same slope regardless of whether the lattice is in the low– or high–temperature regime (Fig. 3b).

The lack of a dependence on temperature indicates that this regime of transport occurs in a lightly corrugated landscape — that is, with low–lying barriers separating the minima associated with thermal transport. In this case, the characteristic barrier scale is below 15 meV, and thus the mean energy at the lowest background temperature ($T_{\text{B}}=230$ K) is above the landscape corrugation. Lower temperature observations are necessary to identify the precise scale, requiring an accurate treatment of quantum effects and different experimental protocols. Stated more succinctly, the equality of the low– and high–temperature diffusivity indicates that the characteristic time $\Omega^{-1}$ is the same and no free energy barrier exists at this level of landscape hierarchy.

Region B. As $\Delta_{ij}T_{\text{eff}}$ goes above 15 K, the flux decreases with the increasing temperature gradient. This suggests the appearance of a vibrational mismatch between adjacent residues due to nonlinearity. That is, adjacent residues separated by a sufficiently large temperature gradient will see different tiers of the energy landscape hierarchy and thus access different vibrational mode structures. As a consequence, the molecular conformation is pushed into an activated region of the free energy landscape where the energy barrier is larger than the available kinetic energy and increases with $\Delta_{ij}T_{\text{eff}}$. Moreover, the average temperature elevation does not substantially change for $\Delta_{ij}T_{\text{eff}}$ in region B where the flux dips (Fig. 3c). Thus, barrier crossing is not aided by energy remaining from the initial deposition. This is further supported by the separation of low– and high–temperature curves, indicating that transport increases with temperature — a signature of a free energy barrier. The characteristic barriers can be estimated from the ratio of high– and low–temperature fluxes (or rates), $J_{\text{H}}/J_{\text{L}}\approx 1.2\approx\exp(-\Delta F/k_{\text{B}}T_{\text{H}}+\Delta F/k_{\text{B}}T_{\text{L}})$, giving values of $\Delta F$ that span from $28$ meV to $167$ meV when we use the average temperature in each regime (i.e., $T_{\text{L}}=250$ K and $T_{\text{H}}=310$ K). These effective barriers are precisely the energy scale leading to conformational changes that restore efficacious vibrational coupling.

Region C. As $\Delta_{ij}T_{\text{eff}}$ increases beyond 30 K, there is a substantial increase in flux for both low– and high–temperature structures. In this case, a large $\Delta_{ij}T_{\text{eff}}$ implies a larger average temperature elevation for a given residue pair (Fig. 3c), as large gradients are primarily found at early times (and near the heater site) when a substantial fraction of initially deposited energy is present (Fig. 3d). If we assume $\Omega$ remains the same, the temperature elevation $\Delta T_{\text{B},j}$ is enough to once again put transport in a stable regime of the landscape at this level of hierarchy, with a typical barrier energy of $67$ meV. This yields an approximately linear region for $J_{\text{BB}}$ with a diffusivity $D_{\text{eff},\text{C}}=1.9\times 10^{-2}$ nm2 ps*-1* ($\kappa_{\text{eff},\text{C}}=3.2\times 10^{-3}$ eV K*-1* ps*-1*).

Region D. Increasing $\Delta_{ij}T_{\text{eff}}$ even further, beyond 50 K, leads to a transport region with a larger diffusivity $D_{\text{eff},\text{D}}=8.0\times 10^{-2}$ nm2 ps*-1* ($\kappa_{\text{eff},\text{D}}=1.3\times 10^{-2}$ eV K*-1* ps*-1*), corresponding to over–the–barrier diffusion. In this case, a new level of the energy landscape hierarchy becomes accessible, which would otherwise require strong activation at lower energies.

Figure 4a shows the effective free energy barriers in the different regimes, which are also reflected in the backbone rate constants (Fig. 4b,c). The $k_{\text{BB}}$ initially decrease with $\Delta_{ij}T_{\text{eff}}$ (from 0 K to 4 K) due to a diminishing residual ballistic component when averaging at $\Delta t=100$ fs. Overestimation of this signature (e.g., through an improper coarse–graining scale), can lead to the discrepancies with experiment found in earlier theoretical analyses Botan2007 ; Nguyen2010 . This is followed by a plateau in $k_{\text{BB}}$ at about 1.5 ps*-1* between 4 K to 15 K, followed by a drop as the landscape is pushed into a new, barrier–dominated region. After this, though, the larger $\Delta_{ij}T_{\text{eff}}$ correspond to a larger temperature elevation, bringing the events above the features in the energy landscape and raising $k_{BB}$ further. {Our methods extract the dependence on the local temperature gradients and, by spatiotemporal correlation, the temperature elevation. Beyond $\Delta_{ij}T_{\text{eff}}=77$ K, the rate constants and fluxes decline sharply, reflecting very early dynamics where strong dynamical localization processes dominate. These barriers collectively define the energy scales, and thus the rate of diffusion in conformational space Best2006 , that is associated with the mechanical dynamics of heat propagation at different temperatures.

Discussion

While our NEMD simulations support that a transition Botan2007 ; Backus2008b ; Backus2009 ; Nguyen2010 in diffusivity is present, they do not support that the transition happens solely due to the existence of energy barriers, as stated in Ref. Botan2007 ; Backus2008 ; Kobus2010 , or glassy dynamics (which is certainly the case but does not pinpoint the particular processes that occur here). Rather, the transition is due to the development of region C physics: Energy flow, which largely happens from 0 to 10 ps, is in the presence of large $\Delta_{ij}T_{\text{eff}}$ (see initial time, high gradient line in Fig. 3d) on top of equilibrium fluctuations ($\Delta_{ij}T_{\text{eff}}\approx\pm 10$ K). We interpret this to indicate that large gradients give a vibrational mismatch via nonlinear energy localization, introducing a barrier to energy transport. In this context, localization would then mediate the transition into a higher diffusivity regime — thereby suggesting an origin of the sharpness of the transition. The increase of the base temperature reduces the vibrational mismatch by pushing the dynamics onto a different level of the landscape hierarchy. Simultaneous Arrhenius activation and barrier reduction conspire to give a sharp transition. More extensive simulations are necessary to make this precise.

These findings demonstrate that energy transport gives quantitative information regarding the biomolecular free energy landscape, its nonlinearity, and overall connectivity. Going beyond what we present here, the experimental analogues of our simulations offer potential probes of structural transitions, where a temperature–dependent change in the transport profile is a manifestation of the graph–theoretic topology associated with molecular contacts and nonlinear interactions of the dominant conformer(s). In other words, thermal transport can be employed to devise ‘tomographies’ that provide a complementary mapping of biomolecular structure, conformational dynamics, and folding pathways. While dominated by local contacts and secondary structure within the simple Aib10 peptide, we expect higher aspects of fold (tertiary, quaternary) to define these dynamics in increasingly complex biomolecules. Furthermore, such probes might excel for highly fluctuating systems such as intrinsically disordered proteins (IDPs), where efficacious thermal transport may still persist (addressed in the Supplementary Discussion), or as a means to dissect local shifts in vibrational mode structure during molecular signaling or allostery. These dynamics have been impervious to other spectroscopies. Our approach provides the conceptual foundations and analysis tools that are directly applicable to experimental data, permitting the immediate interpretation of measurements that leverage local vibrational thermometry. In addition to the functional implications, the approach will also enable the development of a better understanding of what interactions look like at the atomic scale, and therefore better force-fields, and facilitate the design of nanodevices with directed, environmentally responsive heat transport mechanisms.

Methods

Molecular dynamics simulations. Our simulations consist of a modified Aib10 peptide (AcOHN-Aib10-COOCH3), embedded in a box of 922 chloroform molecules. Equilibration and ensemble generation are described in Ref. Elenewski2018 . Prior to NEMD runs, structures are further equilibrated for 100 ps at each base ($T_{\text{B}}$) temperature (NPT; time step $\delta t=1.0$ fs) followed by a 50 ps run with shorter time step (NPT; $\delta t=0.1$ fs). Using the final configurations, NEMD (NVT; $\delta t=0.1$ fs) is initiated by heating the first residue of Aib10 to $T_{\text{B}}+\Delta T$ ($\Delta T=670$ K) for 1 ps, while holding the remaining atoms at $T_{\text{B}}$. Thermostatting is then disabled and heat propagation monitored in the microcanonical ensemble. Similar thermostatting protocols have been established as surrogates for explicit photoexcitation Nguyen2010 ; Buchenberg2016 . NVT simulations employ a velocity Verlet integrator and modified Nosé–Hoover thermostat (damping = 100 fs), while NPT runs add a Martyna–Tobias–Klein barostat (damping = 1000 fs, eight member chain) Parrinello1981 ; Martyna1994 ; Shinoda2004 . Isotropic cell fluctuations are allowed for NPT runs and initial velocities are assigned according to a Gaussian distribution. Simulations employ CHARMM36 force field parameters MacKerell2004 ; Best2012 , CHARMM pair potentials (without CMAP parameters, as rationalized in Ref. Elenewski2018 ), transferrable parameters for CHCl3 Norbeg1998 , PPPM electrostatics (force cutoff $6.95\times 10^{-3}$ pN; pair coupling rescaled at 1.0 nm, terminated at 1.35 nm) and the LAMMPS codebase Plimpton1995 . We have adopted a thermostat timescale that is faster than backbone amide relaxation and azobenzene isomerization in order to preserve transport–relevant dynamics. While a slight overpopulation of long–range modes remains possible, it would only serve to underestimate the impact of nonlinear localization while overestimating ballistic signatures — thus leaving our conclusions unaffected.

Kinetic fitting. While physically descriptive, the master equation, Equation (1), is underdetermined when fitting the simulated transport profiles $E_{j}(t)=3/2N_{j}k_{\text{B}}\langle T(t)\rangle$ for the $N_{j}$ atoms of the $j^{\mathrm{th}}$ residue. As a simplifying approximation, we relate forward and reverse rate constants $k_{ij}=(f_{i}/f_{j})k_{ji}$ through the degrees of freedom of each residue $f_{j}$, as required for detailed balance to hold at equilibrium. We also restrict analysis to structurally homogeneous (helical) conformers, where the rate constants for hydrogen bond energy transfer $k_{j,j+3}\approx k_{j,j+4}\approx k_{\text{HB}}$ and solvent coupling $k_{\text{s},j}\approx R_{j}\,k_{\text{s}}$ as can be approximated as uniform (up to a fixed geometric factor $R_{j}$ for the surface area of terminal residues). Under these conditions, we may fit the time dependence of the solvent $k_{\text{s}}\rightarrow k_{\text{s}}(t)$ and peptide rate constants, $k_{ij}\rightarrow k_{ij}(t)$ and $k_{\text{HB}}\rightarrow k_{\text{HB}}(t)$, to account for the local temperature (which changes in time). This is in contrast to prior efforts that assume a uniform and time–independent backbone rate constant $k_{j,j+1}=k_{\text{BB}}$ Buchenberg2016 .

Rate constants $\mathbf{k}_{j}=(k_{1,2},\dots,k_{N-1,N},k_{\text{H}})$ at the $n^{\mathrm{th}}$ simulation time step are estimated for the linear system of Equation (1) though a constrained optimization

[TABLE]

where $\mathbf{d}_{j}(t_{n})=[E_{j}(t_{n})-E_{j}(t_{n-1})]+k_{\text{s}}(t_{n})[E_{j}(t_{n})-E_{\text{s}}]$ captures energy redistribution among residues of the peptide. The matrix $\mathbf{G}(t)$ is similarly defined so that $\mathbf{G}_{i,j}(t)=-\mathbf{G}_{i+1,j}(t)=-[E_{i}(t)-E_{j}(t)]$ accommodates backbone energy transport and $\mathbf{G}_{i,N}(t)=\sum_{\ell}[E_{i}-E_{\ell}]$ describes its hydrogen bonding counterpart to the $i^{\mathrm{th}}$ residue. The solvent coupling rate $k_{\text{s}}(t_{n})=\sum_{j}[E_{j}(t_{n})-E_{j}(t_{n-1})]/[E_{j}(t_{n})-E_{\text{s}}(t_{n})]$ is then given by the energy exchanged between the peptide and the solvent at each time step (the solvent bath energy $E_{\text{s}}(t_{n})=3N_{j}k_{\text{B}}T_{\text{B}}/2$ is treated a constant).

Data availability

The authors declare that all data supporting the findings in this manuscript are available within the paper and its supplementary information.

Acknowledgements

The authors would like to thank Thomas LeBrun for his insightful comments. J. E. acknowledges support under the Cooperative Research Agreement between the University of Maryland and the National Institute for Standards and Technology Physical Measurement Laboratory, Award 70NANB14H209, through the University of Maryland. K. V. was supported by the U.S. Department of Energy through the LANL/LDRD Program. Computing resources were made available through the Los Alamos National Laboratory Institutional Computing Program, which is supported by the U.S. DOE National Nuclear Security Administration under contract no. DE-AC52-06NA25396, as well as the Maryland Advanced Research Computing Center (MARCC).

Author contributions

J. E. performed the simulations and analysis. J. E. and M. Z. formulated the theoretical concepts. J. E., K. V., and M. Z. all contributed to the development of the ideas and preparation of the manuscript.

Competing interests

The authors declare no competing interests.