Efficient conversion of chemical energy into mechanical work by Hsp70 chaperones

To characterize the main features of chaperone-induced expansion, we performed MD simulations of the Hsp70/rhodanese complexes. We relied on a one-bead-per-residue Coarse Grained (CG) force field [15], which has been tailored to match experimental FRET data of intrinsically-disordered proteins and satisfactorily reproduce the compactness of unfolded rhodanese in native conditions without any further tuning (see SI). In particular, we focused on hydrophobic and excluded volume interactions while neglecting the electrostatic contribution which is negligible in rhodanese and plays only a minor role in Hsp70/rhodanese complexes according to FRET experiments (see SI). Hsp70 chaperones were modeled with a structure-based potential built upon the ADP-bound conformation, and restrained on binding sites on the substrate. We identified six binding sites on the rhodanese sequence using two distinct bioinformatic algorithms[17, 18]. Considering that each binding site could be either free or bound to a Hsp70 protein, we thus took into account a total of $2^{6}=64$ distinct chaperone/substrate complexes, which were exhaustively simulated. In Fig.1 we report the distributions of the substrate potential energy and of the radius of gyration ($R_{g}$) for three representative complexes with one (left), three (center) and six (right) bound chaperones. As previously noticed[16], chaperone binding leads to larger radii of gyration and higher potential energies, implying that the excluded volume interactions due to the large Hsp70s progressively expand the complex and disrupt the attractive intra-chain interactions in rhodanese.

We then calculated the conformational free energy of all the possible chaperone/rhodanese complexes to obtain a quantitative picture of the energy landscape governing the chaperone-induced expansion. To this aim, we performed extensive sets of non-equilibrium steering MD trajectories for each complex, and measured the work needed to steer it to a completely extended reference structure ($R_{gyr}>260\mbox{\AA }$), whose conformational free-energy is not affected by chaperone binding. Equilibrium free-energy differences with respect to this reference state were then estimated from non-equilibrium work distributions via the Jarzynski equality[36], thus allowing the determination of the conformational free-energy $\Delta G$ of each distinct chaperone/substrate complex.

In Fig.2 (main) we report $\Delta G$ for each complex as a function of its mean radius of gyration using different colors for different stoichiometries. The conformational free energy increased with the swelling of the substrate due to the progressive binding of the chaperones. The increase in substrate potential energy due to the loss of intra-chain interactions upon Hsp70 binding is therefore only marginally compensated by the gain in conformational entropy. Notably, the conformational free-energy is not uniquely determined by the stoichiometry, and is significantly affected by the specific binding pattern. The conformational free-energy cost $\Delta\Delta G$ of adding a single chaperone (inset in Fig.2) is positive for all complexes, but it varies from 2 kcal/mol up to 7 kcal/mol depending on the stoichiometry of the complex and on the particular choice of the binding sites. The increase of $\Delta G$ as a function of $R_{g}$ is quantitatively captured by Sanchez theory ([31] and SI) for the coil-to-globule collapse transition in polymers (see Fig.2). Remarkably, the excellent agreement is not the outcome of a fitting procedure since all the parameters were extracted from experiments (see Methods). This result further reinforces the reliability of our simulations as well as the general applicability of the present setup beyond the particular system considered in this work.

The structural and thermodynamic characterization obtained by molecular simulations can be profitably complemented by a kinetic model encompassing relevant biochemical processes in order to determine the probability of each chaperone/substrate complex as a function of the chemical conditions. Notably, a model of the Hsp70 biochemical cycle based on experimental rates was previously used to illustrate how ATP-hydrolysis may result into non-equilibrium ultra-affinity for peptide substrates [10]. Here we extend this result to the more complex case of Hsp70-induced expansion by taking into account multiple chaperone binding events and their consequences on the conformational free energy of the substrate.

In our model each state corresponds to a single configuration of the chaperone/substrate complex, which is defined by the occupation state of the six Hsp70 binding sites on rhodanese. Each site can be either free or occupied by an ADP- or ATP-bound chaperone for a total of $3^{6}=729$ different states. All the relevant molecular processes corresponding to transitions between these states are explicitly modeled, including chaperone binding/unbinding, nucleotide exchange and ATP hydrolysis (see Fig.3). We took advantage of available biochemical data for determining the rate constants associated to all the relevant reactions (see methods and SI). Importantly, kinetic rates for Hsp70 binding were modulated by the conformational free-energies determined by CG MD simulations. Indeed, the unbinding rates of Hsp70 from large-sized protein substrates were observed to be similar to the ones from small peptides, whereas the binding rates can be up to two orders of magnitude smaller[16, 9, 26]. This evidence was further corroborated by a recent NMR study[12] suggesting a conformational selection scenario where the energetic cost due to substrate expansion mostly affects the Hsp70/rhodanese binding rate. Following the experiments, we thus considered a substrate-independent unbinding rate constant $k_{off}$, while we expressed the binding rate constant as

where $\beta=1/k_{b}T$, $k^{0}_{on}$ is the binding rate measured for a peptide substrate, and $\Delta\Delta G_{ij}$ is the conformational free-energy cost of Hsp70 binding, which depends on the specific initial and final binding patterns $i$ and $j$ in the rhodanese/chaperone complex (see Fig.2, inset). The interactions with JDP cochaperones were not explicitly modeled but the cochaperones were assumed to be colocalized with the substrate, so that their effect was implicitly taken into account in the choice of the rate constants for the ATP hydrolysis [34, 14].

The analytical solution of the model provides the steady-state probability of each binding configuration and allows the exploration of their dependence on the biochemical parameters. It is particularly instructive to investigate the system behavior as a function of the ratio between the concentration of ATP and ADP, which is intimately connected to the energy released by ATP hydrolysis. At thermodynamic equilibrium, the $[ATP]/[ADP]$ ratio is greatly tilted in favor of ADP ($[ATP]_{eq}/[ADP]_{eq}\simeq 10^{-9}-10^{-8}$, [21]) whereas in the cell ATP is maintained in excess over ADP by energy-consuming chemostats ($[ATP]/[ADP]>1$, [22]). The $[ATP]/[ADP]$ ratio hence determines how far the system is from equilibrium, thus representing a natural control parameter for the non-equilibrium biochemical cycle. We thus report in Fig.4 the compound probabilities for complexes with the same stoichiometry $n$ as a function of this nucleotide ratio. In conditions close to equilibrium (very low values of $[ATP]/[ADP]$), the vast majority of the substrate proteins are free and only about 10% of them are bound to a single chaperone. The population of equimolar complexes increases for $[ATP]/[ADP]$ between $10^{-2}$ and $10^{-1}$ and gives way to larger complexes with multiple chaperones for higher values of the nucleotide ratio.

For $[ATP]/[ADP]>1$, most substrates are bound to at least $4$ chaperones, with an average stoichiometry $\left<n\right>\sim 4.9$. Further increase of the nucleotide ratio does not significantly change this scenario indicating an almost constant behaviour in large excess of ATP $([ATP]/[ADP]>10)$.

Combining the steady-state probabilities derived from the rate model with the results of the MD simulations, we can now exhaustively characterize the structural properties of the system. This provides the opportunity to directly compare our model with the results from FRET experiments both in equilibrium and non-equilibrium conditions. To this aim, we first focused on the average radius of gyration of the system at thermodynamic equilibrium ($[ATP]\ll[ADP]$) or in non-equilibrium conditions with ATP in large excess over ADP ($[ATP]/[ADP]>10$).

We also probed the robustness of this result with respect to inaccuracies in the molecular model by taking into account normally-distributed errors on the conformational free energies $\Delta G_{i}$. The results are reported as histograms in Fig.5 (top) and they suggest that at equilibrium the average radius of gyration is extremely close to what would be measured in the case of free substrate (dashed line). This is in agreement with the experimental observation that Hsp70 cannot significantly associate to rhodanese in these conditions[16]. Conversely, in large excess of ATP we observe a substantial swelling of the substrate ($75<R_{g}<95$Å) due to the ultra-affine binding of Hsp70s. This finding is fully compatible with the size of DnaK/DnaJ/rhodanese complexes determined by sm-FRET experiments in excess of ATP [16]. In this regime, the limited effects of cochaperone binding on substrate conformations, which are not explicitly included in the model, play a minor role in determining the global expansion of the complex.

A more quantitative comparison between the model and the FRET results can be achieved by back-calculating the transfer efficiencies that were experimentally measured for five distinct pairs of fluorescent dyes[16]. In equilibrium conditions, namely when $[\mbox{ATP}]/[\mbox{ADP}]\ll 1$, the calculated FRET efficiency is $\simeq$ 0.8 for all considered pairs of fluorescent dyes (blue circles) and it matches the experimental results for the compact unbound rhodanese ( 0.8). A dramatic difference is observed in excess of ATP (red circles), where the expansion of the substrate leads to a significant decrease of the calculated efficiency, in excellent agreement with the experimental values measured in similar conditions (black circles)[16]. Remarkably, the results correctly captured the non-monotonic behaviour of FRET efficiency as a function of the sequence separation between the dyes, which was not reproduced in previous calculations[16]. This agreement corroborates the prediction of the DnaK binding sites on the rhodanese sequence and the overall reliability of our model.

Molecular chaperones consume energy via ATP hydrolysis in order to expand rhodanese. It is hence important to determine how effective they are as molecular machines, as well as to assess how favourable the physiological conditions are to perform their biological task.

To this aim, we calculated the global increase in the overall conformational free energy of the substrate with respect to equilibrium conditions, $\Delta G_{Swell}$(Fig.6, top). This quantity measures the excess probability of each complex with respect to equilibrium conditions weighted by its corresponding conformational free-energy $\Delta G_{i}$.

where $p_{i}\left(\frac{[ATP]}{[ADP]}\right)$ is the probability of complex $i$ for a given value of $[ATP]/[ADP]$ and $p^{eq}_{i}$ is the same quantity computed at equilibrium conditions. In order to investigate the conversion of chemical energy into mechanical work it is instructive to focus on the ratio between $\Delta G_{Swell}$ and the free energy of hydrolysis of ATP $\Delta G_{h}$, which reports on the effectiveness of the transduction process. We plot in Fig.6 (top) this quantity as a function of the $[ATP]/[ADP]$ ratio considering the estimated inaccuracies of the model as previously done for the gyration radius. Not surprisingly, all these curves exhibit a maximum because the probabilities of the different states, and thus also $\Delta G_{swell}$, attain plateaus for $[ATP]\gg[ADP]$, whereas $\Delta G_{h}$ increases monotonically with the nucleotide ratio (see Methods). The maximal transduction regime intriguingly corresponds to values of $[ATP]/[ADP]$ that are typical of cellular conditions (grey area).

We highlight that in our model Hsp70 functioning encompasses two distinct yet intertwined processes: the ATP-dependent binding of the chaperones to the substrate, and its consequent expansion. Energy transduction occurs then through two steps, and the amount of energy available for the mechanical expansion is limited by that provided by chaperone binding. To corroborate this picture and to gain a more complete insight on the action of Hsp70s, we analyzed the energetic balance of chaperone binding to a single site. To this aim, we focused on a simplified reaction cycle, which essentially corresponds to a single triangle within the overall scheme in Fig.3 and does not take into account the conformational free-energy of the substrate. We report in Fig.6 (lower panel, black) the non-equilibrium dissociation constant, $K^{neq}_{d}$, normalized with respect to its equilibrium value $K^{eq}_{d}$, as a function of [ATP]/[ADP]. When the ratio between the concentrations of ATP and ADP approaches the physiological regime, the dissociation constant drops significantly until it settles at a value that is two order of magnitude lower than its equilibrium counterpart (this result is the core of ultra-affinity). We can convert the dissociation constant into a binding free energy excess with respect to equilibrium

that we can compare to the free-energy of ATP hydrolysis, $\Delta G_{h}$, as previously done in the case of $\Delta G_{Swell}$. Interestingly, also in this case the energy ratio is maximal in cellular conditions, suggesting that the optimality of the overall expansion process (Fig.6, top panel) does not depend on specific features of the substrate but it is a direct consequence of the intrinsic kinetic parameters of Hsp70 binding.

In our MD simulations both rhodanese and Hsp70 were coarse grained at the single-residue level. To this aim, the molecules were represented as collections of beads centered on the $C_{\alpha}$ atom of each amino acid. Rhodanese was modeled according to the force field from Ref. [15] (see SI for details). We modeled Hsp70 starting from the known crystal structure of ADP-bound Hsp70 (PDB:2KHO [37]). The influence of Hsp70 in the ATP state on substrate conformation was assumed to give similar results. Following [35], we considered both the NBD and the SBD to be rigid bodies interacting only via excluded-volume interactions, while we modeled the flexible linker in the same way as rhodanese. The residues of the binding site moved rigidly with the corresponding SBD, thus ensuring that each chaperone was irreversibly bound to the substrate. All the simulations were performed with a version of LAMMPS [38] patched with PLUMED [39]. The temperature $T=293\,\mbox{K}$ was controlled through a Langevin thermostat with damping parameter $16\,\mbox{ns}^{-1}$. The time step was set equal to 1 fs, and each residue had a mass equal to 1 Da to speed up equilibration. In equilibrium simulations, the system was equilibrated for $10^{7}$ time steps starting from a rod-like conformation, and subsequently sampled for other $10^{7}$ time steps. At least 10 independent simulations were performed for each configuration. Statistical errors on the computed quantities were estimated as standard errors of the mean. In the pulling simulations, starting from an equilibrium conformation the substrate was pulled by an external force acting on its radius of gyration $R_{g}$ at a constant pulling speed $v=10^{-5}$ Å$/$fs, until a rod-like conformation was reached, which we defined to be at $R_{g}=260$ Å. For this set of simulations, 100 independent realizations for each configuration were performed. The indetermination of the computed free energies were estimated according to the bootstrap method. For both equilibrium and pulling simulations, the statistical errors are smaller than the size of symbols reported in the figures.

For the rate model we considered all possible binding configurations for chaperones in the ATP or ADP state to the six identified binding sites, resulting in a total of $3^{6}=729$ states. The equation for the average concentration of each binding configuration $c_{i}$ has the form

where $k_{ij}$ is the transition rate from configuration $i$ to $j$. The transition rates between different states depend on the intrinsic rates of the underlying molecular processes, namely binding and unbinding of a chaperone from a given binding side, the hydrolysis and synthesis of the nucleotides, and the nucleotide exchange, and have been taken according to biochemistry experiments [16, 34, 43] (see SI for the rate values and for details about the thermodynamic constraints between the rates).

We considered the steady state solution, obtained by imposing

for every $i$. The concentration of chaperones was fixed at $[\mbox{Hsp70}]=10\,\mu M$, and the concentration of nucleotides was equal to $1\,\mbox{mM}$ as in the experiments in [16]. The free energy provided by the hydrolysis of ATP is given by the following formula:

where $[ATP]_{eq}$ and $[ADP]_{eq}$ are the concentrations of ATP and ADP in equilibrium. The mathematical and physical details of the model, as well as the values of the rates, are described in the Supplementary Information.

Molecular graphics in Figs.1 and 6 have been generated with UCSF Chimera, developed by the Resource for Biocomputing, Visualization, and Informatics at the University of California, San Francisco, with support from NIH P41-GM103311 [44].