Topology, landscapes, and biomolecular energy transport

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 [1,2], mediate allosteric communication [3][4][5], and drive bioenergetic processes [6][7][8]. 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 [8]. 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 = −κ∇T and its time-dependent version capture diffusive heat flow, given by the flux J, in response to a temperature gradient ∇T . Those two quantities are related by the thermal conductivity κ (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 [9,10]. 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.
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 allatom molecular dynamics alone [21][22][23][24][31][32][33] or normalmode analysis (even when treating anharmonicity as a correction) [34][35][36][37][38][39][40]. The former does not unravel the atomic-scale mechanisms of transport and the latter reflects dynamics only at potential energy minima [37]. Within this context, we further demonstrate how the graph-theoretic topology of molecular contacts can define directed pathways for molecular energy redistribution.  States for a right-handed helix are colored from blue (more chiral) to green (less chiral), while those of a left-handed helix are uniformly grey. Helicity parameters and ensemble determinations follow Ref. [11]; (b) Major conformers in the Aib10 structural ensemble. The C-terminal heater residue is denoted by a red asterisk (*), and hydrogen bonds are colored green; (c) Thermal transport profile from NEMD simulations, characterized as a per-residue kinetic temperature elevation ∆TB,j(t) = Tj(t) − TB with respect to the solvent bath. The dashed, white line demarcates the ballistic front; (d) Differential heat transport between a full structural ensemble and those (∆T λ B,j ) containing only λ = helical, hairpin, or unstructured populations. Upper and lower limits on the temperature elevation (e.g., on ∆TB,j) provide a cutoff for all values lying above or below the bound, respectively.

Results
Topology and energy propagation pathways. We initiate our investigations using a series of replicaexchange 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 [41][42][43]. Our simulation system is a ten-residue Aib helix (Aib 10 ) solvated by chloroform, similar to experimental efforts [13,[17][18][19][20]. We previously generated temperaturedependent free energy landscapes for Aib 10 at high resolution with replica-exchange simulations [11]. From the resulting conformational ensemble, we extract 4000 conformers for each environmental (bath) temperature T 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 ≈ 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 = T B + ∆T , with ∆T = 670 K, while holding the remainder of the system at T B . The simultaneous heating of all vibrational degrees of freedom in the heater residue is wellfounded, as it yields thermal transport profiles that are indistinguishable from mode-selective heating [13,21]. This excess energy then propagates freely within the mi-crocanonical ensemble (i.e., without thermostatting).
The conformational ensemble of Aib 10 comprises three general structural motifs (Fig. 1b) corresponding to (i) 3 10 -/α-helical conformers (≈ 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 Aib 10 (≈ 15%); and (iii) unstructured or extended conformers that have no consistent hydrogen bonding (≈ 40%) [11]. We index these subensembles with λ. This partition is defined by the underlying free energy landscape, and is thus independent of our thermal transport simulations [11].
In Fig. 1c,d, we present transport profiles for Aib 10 versus the ensemble-averaged temperature elevation ∆T B,j (t) = T j (t) − T B of the j th residue, or ∆T λ B,j − ∆T B,j for subensemble λ. 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 [13,[26][27][28]44]. While this channel is weak, additional ballistic pathways may exist at lower group velocities in different vibrational bands [44,45], 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 temperaturedependent thermal diffusivity D(T B ) has distinct lowand high-temperature regimes (Fig. 2a), which are also reflected in the net heat transfer (Fig. 2b). This qualitative behavior agrees with experimental [13,18,20] and theoretical [13,21,22] efforts. These, though, report diffusivities of 0.02 nm 2 ps −1 and 0.1 nm 2 ps −1 , respectively. Theoretical D(T B ) from this type of estimate consistently exceed experimental values for Aib 10 but are comparable to bulk materials [28] and other proteins [40]. 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 [13,18,20,21]. 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 [13,21]. Figure 1d shows the transport profile of the full Aib 10 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 ≤ 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 [33,40]. 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 [46][47][48].
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 intermediatetimescale dynamics of NEMD trajectories -restricting to helical Aib 10 conformers for both structural heterogeneity and consistency with prior work -using a master equation for the kinetic energy E j of the j th residue in the peptide: 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 s,j (t) is the rate of heat transfer to the solvent bath, and E s,j (t) is the kinetic energy density of the solvent surrounding the j 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 timethereby implying a temperature dependence. This accommodation is key to our subsequent analysis. Given this arrangement, one can identify two distinct intrapeptide couplings: (i) direct transfer between nearestneighbors 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 α-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 [33]. 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 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 ≤ N bins via block averaging. This approach is a finite difference decomposition of the diffusion equationĖ(x, t) = D ∇ 2 E(x, t) at the timescale ∆t = t n+1 − t n and a length-scale ∆x defined by the distance between adjacent residues. The fluxes come from the finite difference decomposition of The rate constants k i,j (t n ) = D(t n )/(∆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 prop- agation. This latter property is reflected by the local, activated conformational changes underlying transport where ∆G i,j is the free energy barrier between heat-accepting microstates and (Ω 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 ∆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 (∆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 BB is parameterized by an effective tempera- 3N k B between residues i and j, where the flux is incident on a residue containing N atoms. While transport is explicitly quantified through J BB for simplicity, the effect of hydrogen bonding is present when fitting the backbone flux distribution at hydrogen bonding sites. The results for J BB are presented in Fig. 3a. A complimentary analysis for J HB and a validation of fitting methods are presented in Supplementary Figures 10-15.
Region A. The forward flux J BB has a linear region for small ∆ ij T eff (less than about 15 K), although it does not go to zero at ∆ ij T eff = 0. Purely diffusive transport will not afford a heat flux in the absence of a local temperature gradient. Thus, a finite J BB at ∆ ij T eff = 0 is a signature of ballistic/coherent behavior. Supporting this interpretation, we find that the zero-gradient flux to decrease with increasing ∆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 eff,A = 2.3 × 10 −2 nm 2 ps −1 (or conductivity κ eff,A = 3.9×10 −3 eV K −1 ps −1 ). Fitting for small ∆ ij T eff , while ignoring the residual ballistic contribution right around ∆ ij T eff = 0, removes high rate constant artifacts. Encouragingly, the magnitude of the resulting diffusivity is consistent with experimental values [13,17]. 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 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 hightemperature diffusivity indicates that the characteristic time Ω −1 is the same and no free energy barrier exists at this level of landscape hierarchy.
Region B. As ∆ ij T 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 ∆ ij T eff . Moreover, the average temperature elevation does not substantially change for ∆ ij T 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 Effective free energy barriers ∆F corresponding to different regions of the JBB flux profile. Region A has nearly no barrier, but as the gradient becomes large, a barrier starts to form and increases in region B. In C, this barrier decreases until in D it is zero to within statistical errors (albeit, the uncertainty is large in this last region due to the limited number of samples for large temperature gradients, which occur only at short times). (b) Backbone rate distributions (kBB) for helical Aib10 conformers. Rate constants are parameterized by the temperature gradient ∆ijT eff between adjacent residues and (c) partitioned into low-temperature (blue; 230 K to 270 K) and high-temperature (red; 290 K to 330 K) regimes. The error bars are plus/minus one standard error.
further supported by the separation of low-and hightemperature 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), , giving values of ∆F that span from 28 meV to 167 meV when we use the average temperature in each regime (i.e., T L = 250 K and T H = 310 K). These effective barriers are precisely the energy scale leading to conformational changes that restore efficacious vibrational coupling.
Region C. As ∆ ij T eff increases beyond 30 K, there is a substantial increase in flux for both low-and hightemperature structures. In this case, a large ∆ ij T 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 Ω remains the same, the temperature elevation ∆T 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 BB with a diffusivity D eff,C = 1.9 × 10 −2 nm 2 ps −1 (κ eff,C = 3.2 × 10 −3 eV K −1 ps −1 ).
Region D. Increasing ∆ ij T eff even further, beyond 50 K, leads to a transport region with a larger diffusivity D eff,D = 8.0 × 10 −2 nm 2 ps −1 (κ eff,D = 1.3 × 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 BB initially decrease with ∆ ij T eff (from 0 K to 4 K) due to a diminishing residual ballistic component when averaging at ∆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 [13,21]. This is followed by a plateau in k 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 ∆ ij T 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 ∆ ij T 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 [49], that is associated with the mechanical dynamics of heat propagation at different temperatures.

Discussion
While our NEMD simulations support that a transition [13,18,20,21] in diffusivity is present, they do not support that the transition happens solely due to the existence of energy barriers, as stated in Ref. [13,17,22], 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 ∆ ij T eff (see initial time, high gradient line in Fig. 3d) on top of equilibrium fluctuations (∆ ij T eff ≈ ±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. Simultane-ous 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 Aib 10 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.

Molecular dynamics simulations.
Our simulations consist of a modified Aib 10 peptide (AcOHN-Aib 10 -COOCH 3 ), embedded in a box of 922 chloroform molecules. Equilibration and ensemble generation are described in Ref. [11]. Prior to NEMD runs, structures are further equilibrated for 100 ps at each base (T B ) temperature (NPT; time step δt = 1.0 fs) followed by a 50 ps run with shorter time step (NPT; δt = 0.1 fs). Using the final configurations, NEMD (NVT; δt = 0.1 fs) is initiated by heating the first residue of Aib 10 to T B + ∆T (∆T = 670 K) for 1 ps, while holding the remaining atoms at T B . Thermostatting is then disabled and heat propagation monitored in the microcanonical ensemble. Similar thermostatting protocols have been established as surrogates for explicit photoexcitation [21,33]. NVT simulations employ a velocity Verlet integrator and modified Nosé-Hoover thermostat (damping = 100 fs), while NPT runs add a Martyna-Tobias-Klein barostat (damp-ing = 1000 fs, eight member chain) [50][51][52]. Isotropic cell fluctuations are allowed for NPT runs and initial velocities are assigned according to a Gaussian distribution. Simulations employ CHARMM36 force field parameters [53,54], CHARMM pair potentials (without CMAP parameters, as rationalized in Ref. [11]), transferrable parameters for CHCl 3 [55], PPPM electrostatics (force cutoff 6.95 × 10 −3 pN; pair coupling rescaled at 1.0 nm, terminated at 1.35 nm) and the LAMMPS codebase [56]. 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 B T (t) for the N j atoms of the j 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 ≈ k j,j+4 ≈ k HB and solvent coupling k s,j ≈ R j k 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 s → k s (t) and peptide rate constants, k ij → k ij (t) and k HB → k HB (t), to account for the local temperature (which changes in time). This is in contrast to prior efforts that assume a uniform and timeindependent backbone rate constant k j,j+1 = k BB [33].
Rate constants k j = (k 1,2 , . . . , k N −1,N , k H ) at the n th simulation time step are estimated for the linear system of Equation (1) though a constrained optimization captures energy redistribution among residues of the peptide. The matrix G(t) is similarly defined so that is then given by the energy exchanged between the peptide and the solvent at each time step (the solvent bath energy E s (t n ) = 3N j k B T B /2 is treated a constant).

Data availability
The authors declare that all data supporting the find-ings in this manuscript are available within the paper and its supplementary information.  There is a small negative Q(t). This is within the error, but may be due to the high flexibility of residue 10 leading to intermittent contact with the remainder of the peptide even with in the helical structure. This could give a thermal transport pathway from residue 10 to 9, or some residual effect from the solvent. The error bands are plus/minus one standard error. n P E y q 1 D m r / P t z Y 2 y I 2 X O + d Y m 5 T 7 r o U T 6 O G w U 3 4 J m I r M I Y p R a W F v j L M V 6 L g w q q a F T C D I L q Z C X 4 h x q j s Y i A l v P e x 9 0 6 L J T m r S V G F c x 0 p 7 6 f S I X k b X t K H S d k c A L + 9 v r i n 9 5 t Q x b W / V c x W m G E M v P R a 1 M U 0 x o N w 7 a V A Y k 6 r Y j Q h r l b q X y Q h g h 0 Y V W c C F 8 f U r / J y d r P l / 3 2 d F G a W e l H 8 c Y W S R L Z I V w U i E 7 Z J 8 c k i q R 5 J b c k 0 f y 5 N 1 5 D 9 6 z 9 / L Z O u D 1 Z + b J D 3 i v H y Y h l V w = < / l a t e x i t >   25 50 75 0 ΔF 0 T B  270 K < l a t e x i t s h a 1 _ b a s e 6 4 = " P t n e I T k f c s c F H a p n x l d k V U t s e k k = " > A A A C A X i c d V D L S g N B E J z 1 b X x F v Q h e B o P g Q Z Y Z H 4 n e R C + C F w W j Q j a E 2 U l H B 2 c f z v S K Y Y k X f 8 W L B 0 W 8 + h f e / B s n M Y K K F j Q U V d 1 0 d 4 W p V h Y Z e / c G B o e G R 0 b H x g s T k 1 P T M 8 X Z u R O b Z E Z C V S Y 6 M W e h s K B V D F V U q O E s N S C i U M N p e L n X 9 U + v w V i V x M f Y T q E e i f N Y t Z Q U 6 K R G c e G 4 s U s D D V d 0 r c J o s E o D h B v M D z q N Y o n 5 n P E y q 1 D m r / P t z Y 2 y I 2 X O + d Y m 5 T 7 r o U T 6 O G w U 3 4 J m I r M I Y p R a W F v j L M V 6 L g w q q a F T C D I L q Z C X 4 h x q j s Y i A l v P e x 9 0 6 L J T m r S V G F c x 0 p 7 6 f S I X k b X t K H S d k c A L + 9 v r i n 9 5 t Q x b W / V c x W m G E M v P R a 1 M U 0 x o N w 7 a V A Y k 6 r Y j Q h r l b q X y Q h g h 0 Y V W c C F 8 f U r / J y d r P l / 3 2 d F G a W e l H 8 c Y W S R L Z I V w U i E 7 Z J 8 c k i q R 5 J b c k 0 f y 5 N 1 5 D 9 6 z 9 / L Z O u D 1 Z + b J D 3 i v H y Y h l V w = < / l a t e x i t > Δ ij T eff (K)

II. SUPPLEMENTARY DISCUSSION
As a general rule, different folds of the Aib 10 peptide (or other biomolecules) have a conformationallydependent solvent-accessible surface area (SASA). Since thermal conduction between the peptide and the solvent bath occurs at this interface, it is expected that the overall thermalization profile will also depend on molecular conformation. To quantify this, we exploit the largely diffusive nature of heat conduction in Aib 10 to construct a simple model for thermal relaxation. In this case, we assume two compartments -consisting of the peptide and the solvent -which undergo strictly conductive heat transfer according to Fourier's law (i.e., no convective contribution). Assuming that the solvent bath is much larger than the peptide, with minimal local solvent heating, we can write the total kinetic energy content E(t) of Aib 10 in the time-dependent form where E B is the net kinetic energy content of Aib 10 when in thermal equilibrium with the bath, E(0) is the net kinetic energy content of peptide immediately following heating, and k c is a characteristic Our simple cooling model is generally robust when applied to thermal relaxation in Aib 10 .
Modest deviations between the resulting fits and simulation data are observed at early times (t ≤ 2.5 ps) and high temperatures, where ballistic processes likely shunt heat to the solvent more rapidly than allowed by a diffusive mechanism. Outside of this region, the thermal transport dynamics are relatively similar for the helix, hairpin, and extended coil, with the most prominent (relative) conformational variation observed for the cooling rate constant k c . In our simulations, the extended coil dissipates heat most rapidly to solvent, followed by the structured α-helix fold and comparatively globular hairpin conformations. While statistically significant, this effect is also small -largely due to the fact that Aib 10 dynamics are surface-dominated in any conformer due to finite-size effects.
The overall scaling trend for k c is mirrored when comparing the mean SASA between conformers (Supplementary Table 1). Quantitative differences nonetheless exist between SASA and k c ratios [k c (coil)/k c (hairpin) = 1.10, while we find SASA(coil)/SASA(hairpin) = 1.26], with the surface 18 area contribution underestimated in the rate constants. While a variety of factors may collude in this effect, the conformer-dependent heat transport rates within the peptide likely make the largest contribution. That is, the helical and hairpin conformations possess high thermal diffusivities and auxiliary conduction pathways (due to molecular topology), which facilitate the transfer of heat to away from the heater and throughout the peptide (Supplementary Figure 17). Since the local cooling rateĖ j (t) at the j-th residue is proportional to the local temperature gradient between the peptide and the solventĖ j (t) ∝ T j (t) − T s,j (t), this accelerates solvent relaxation. The peptide thus acts as a radiator, relinquishing heat to the bath while circumventing local solvent heating -and thus partially mitigating the insulating effect of molecular conformation. We expect the contribution of this effect to be less in larger biomolecules -in a manner that depends on the surface-to-volume ratio -where redistribution pathways may also lead deeper into the molecular 'bulk' and thus away from the solvent interface.