Efficient conversion of chemical energy into mechanical work by Hsp70 chaperones

Hsp70 molecular chaperones are abundant ATP-dependent nanomachines that actively reshape non-native, misfolded proteins and assist a wide variety of essential cellular processes. Here, we combine complementary theoretical approaches to elucidate the structural and thermodynamic details of the chaperone-induced expansion of a substrate protein, with a particular emphasis on the critical role played by ATP hydrolysis. We first determine the conformational free-energy cost of the substrate expansion due to the binding of multiple chaperones using coarse-grained molecular simulations. We then exploit this result to implement a non-equilibrium rate model which estimates the degree of expansion as a function of the free energy provided by ATP hydrolysis. Our results are in quantitative agreement with recent single-molecule FRET experiments and highlight the stark non-equilibrium nature of the process, showing that Hsp70s are optimized to effectively convert chemical energy into mechanical work close to physiological conditions.


Introduction
Even though in vitro most proteins can reach their native structure spontaneously (Anfinsen, 1973), this is not always the case in cellular conditions and proteins can populate misfolded states which can form cytotoxic aggregates (Dobson, 2003). In order to counteract misfolding and aggregation, cells employ specialized proteins, called molecular chaperones, which act on non-native protein substrates by processes that stringently depend on ATP hydrolysis for most chaperone families (Hartl, 1996). Among them, the ubiquitous 70 kDa heat-shock proteins (Hsp70s) play a special role because they assist a plethora of fundamental cellular processes beyond prevention of aggregation (Clerico et al., 2019;Rosenzweig et al., 2019).
Decades of biochemical and structural studies have clarified the core elements of the Hsp70 functional cycle at the molecular level (Mayer, 2013). Hsp70s consist of two domains: the substrate binding domain (SBD) interacts with disparate substrate proteins, whereas the nucleotide binding domain (NBD) is responsible for the binding and hydrolysis of ATP. The two domains are allosterically coupled, and the nature of the nucleotide bound to the NBD affects the structure of the SBD and as a consequence the affinity for the substrate and its association/dissociation rates. More precisely, when the chaperone is in the ATP-bound state, the SBD is open and easily accessible to the substrate, whereas the SBD is closed when ADP is bound. These structural differences result in substrate binding and unbinding rates when ATP is bound that are orders of magnitude faster than when ADP is bound (Mayer et al., 2000). Furthermore, the coupling is bidirectional: the substrate, together with a co-localized J-domain protein (JDP) that serves as cochaperone (Kampinga and Craig, 2010;Kampinga et al., 2019), greatly accelerates the hydrolysis of ATP. Substrate binding thus benefits from the fast association rate of the ATP-bound state and the slow dissociation rate of the ADP-bound state, resulting in a non-equilibrium affinity (ultra-affinity) that can be enhanced beyond the maximum limit allowed by thermodynamic equilibrium, namely the affinity of the ADPbound state (De Los Rios and Barducci and De Los Rios, 2015).
More recently, the consequences of Hsp70 binding on the conformational ensembles of its substrates have also been investigated. Several lines of evidence indicate that the binding of Hsp70s to a polypeptide induces its expansion. Biochemical assays revealed that binding of Hsp70 increases the sensitivity of misfolded Luciferase to proteolysis and decreases its propensity to bind Thioflavin-T, strongly suggesting a loss of compactness (Sharma et al., 2010). Nuclear Magnetic Resonance (NMR) measurements have shown that Hsp70s destabilize the tertiary structure of several different substrates (Lee et al., 2015;Sekhar et al., 2015). Moreover, a single-molecule study based on Fö rster resonance energy transfer (FRET) spectroscopy quantified the considerable expansion of unfolded rhodanese in native conditions upon binding of multiple Hsp70 chaperones. In particular, this study revealed that the expansion is stringently ATP-dependent, because upon ATP exhaustion the system relaxes to the expansion values observed in the absence of chaperones (Kellner et al., 2014).
Despite these advances in the characterization of Hsp70 functioning, the mechanistic understanding of how the energy of ATP hydrolysis is used to expand a substrate has lagged behind. Our goal here is precisely to fill this gap between the molecular and functional characterization of Hsp70. To this aim, we first explore the structural and energetic features of Hsp70-bound rhodanese using Molecular Dynamics (MD) simulations. We next integrate this molecular information into a rate model that explicitly includes the Hsp70-rhodanese interactions and the chaperone ATPase cycle, thus elucidating how Hsp70s convert the chemical energy of ATP into mechanical work necessary to expand their substrates.

Structural and thermodynamic characterization of chaperone-substrate complexes
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 (Smith et al., 2014), which has been tailored to match experimental FRET data of intrinsically disordered proteins and satisfactorily reproduces the compactness of unfolded rhodanese in native conditions without any further tuning (see Materials and methods). Hsp70 chaperones were modeled with a structure-based potential to account for their excluded volume and they were artificially restrained to binding sites on the substrate. We identified six binding sites on the rhodanese sequence using two distinct bioinformatic algorithms (Rüdiger et al., 1997;Van Durme et al., 2009). 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 Figure 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. Consistently with FRET results (Kellner et al., 2014), 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 g >260 Å ), 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 (Jarzynski, 1997), thus allowing the determination of the conformational free energy DG of each distinct chaperone/substrate complex (see Figure 2figure supplement 1 and Materials and methods).
In Figure 2 (main), we report DG 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 DDG of adding a single chaperone (inset in Figure 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 DG as a function of R g is quantitatively captured by Sanchez (1979) theory for the coil-to-globule collapse transition in polymers (see Figure 2 and Materials and methods). Remarkably, the excellent agreement is not the outcome of a fitting procedure since all the parameters were extracted from experiments (see Appendix 2). 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.

ATP hydrolysis promotes multiple chaperone binding
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 (De Los Rios and . 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 Potential Energy (kcal/mol) Potential Energy (kcal/mol) Potential Energy (kcal/mol) Figure 1. Progressive expansion of substrate by multiple Hsp70 binding. Probability density maps of substrate potential energy and radius of gyration for representative Hsp70/rhodanese complexes with one (left), three (center) and six (right) bound chaperones. The different Hsp70 chaperones have been represented with different colors to ease their discernibility. The online version of this article includes the following source data for figure 1: Source data 1. Text files containing the data used to generate the probability density maps in Figure 1.
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 Figure 3). We took advantage of available biochemical data for determining the rate constants associated to all the relevant reactions (see Materials and methods). 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 (Schmid et al., 1994;Mayer et al., 2000;Kellner et al., 2014). This evidence was further corroborated by a recent NMR study (Sekhar et al., 2018) suggesting a conformational selection scenario where the energetic cost due to substrate expansion mostly affects the Hsp70/ rhodanese binding rate. Accordingly, we thus considered a substrate-independent unbinding rate constant k off , while we expressed the binding rate constant as Source data 1. Data from simulation results and Sanchez theory used to generate the plot in Figure 2 and the histogram in the inset of Figure 2. Figure supplement 1-source data 1. Simulation data used to generate the free energy plots; the work performed by the pulling force for each trajectory is also included in the subfolders.
where b ¼ 1=k B T, k B is the Boltzmann constant, T is the absolute temperature, k 0 on is the binding rate measured for a peptide substrate, and DDG 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 Figure 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 ATP hydrolysis (Kampinga and Craig, 2010;Hu et al., 2006).
The analytical solution of the model provides the steady-state probability of each binding configuration and allows the exploration of their dependence on external conditions. It is particularly instructive to investigate the system behavior as a function of the ratio between the concentrations 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 ' 10 À9 À 10 À8 ; Alberty, 2005), whereas in the cell ATP is maintained in excess over ADP by energy-consuming chemostats (½ATP=½ADP>1; Milo and Phillips, 2015). 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 Figure 4 (top panel) 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 <n>~4:9 (solid line in Figure 4), bottom panel). Further increase of the nucleotide ratio does not Figure 3. The relevant biochemical reactions of the Hsp70/rhodanese system can be described by a rate model. Each chaperone binding site on rhodanese (black dots) can be either free or occupied by an Hsp70 (yellow), which in turn can be either ADP-or ATP-bound. We depict here for the sake of clarity only a representative portion of the full model, which takes into account six binding sites. The reaction cycle is governed by the rates for chaperone binding/unbinding to the substrate (k ATP on ; k ADP on ; k ATP off ; k ADP off ) and for hydrolysis (k s h ), synthesis (k s s ) and exchange (k eff ex;DT ; k eff ex;TD Þ of nucleotides bound to the chaperones (see Materials and methods for further details). Importantly, the binding rate constants, k ATP on and k ADP on , take into account the conformational free energies, according to Equation (1). significantly change this scenario indicating an almost constant behaviour in large excess of ATP ð½ATP=½ADP>10Þ. It is important here to underscore that the binding of the chaperones in these conditions is a non-equilibrium effect, driven by the Hsp70-induced hydrolysis of ATP, and it is not a mere consequence of the excess of ATP over ADP or of the large difference between the substrate association rates to the ATP-and ADPbound chaperones. Indeed, if we neglect Hsp70 ATPase activity (k h h ; k h s ¼ 0) without changing any of the other model parameters, efficient chaperone binding is abolished (<n> ( 1, as shown in the bottom panel of Figure 4, dashed line). As a matter of fact, in such equilibrium scenario an excess of ATP over ADP actually slightly disfavors chaperone binding, because the Hsp70 affinity for the substrate is slightly lower in the ATPbound state than in the ADP-bound state (k ATP off =k ATP on >k ADP off =k ADP on , see De Los Rios and Barducci, 2014 for further discussion.) 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 ( ½ADP) or in non-equilibrium conditions with ATP in large excess over ADP (½ATP=½ADP>10). In order to probe the robustness of our results with respect to inaccuracies in the molecular model, we also took into account normally distributed errors on the conformational free energies DG i .
The results are reported as histograms in Figure 5 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 the formation of rhodanese-DnaK complexes is strictly dependent on the hydrolysis of ATP and that ADP cannot trigger the expansion of the substrate (Kellner et al., 2014). 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 singlemolecule FRET experiments in excess of ATP (Kellner et al., 2014). 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 (Kellner et al., 2014). In equilibrium conditions, namely when ½ATP=½ADP ( 1, the calculated FRET efficiency is~0.8 for all considered pairs of fluorescent dyes (inset of Figure 5, blue circles) and it matches the experimental results for the compact unbound rhodanese (~0.8). A dramatic difference is instead 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). 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 (Kellner et al., 2014). This agreement corroborates the prediction of the DnaK binding sites on the rhodanese sequence and the overall reliability of our model.

Energy balance and thermodynamic efficiency
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 favorable 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, DG Swell (Figure 6, top). This quantity measures the excess probability of each complex with respect to equilibrium conditions weighted by its corresponding conformational free energy DG i .
where p i ½ATP ½ADP 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 DG Swell and the free energy of hydrolysis of ATP DG h , The ratio DG Swell =DG h reports on the effectiveness of the transduction process. We plot in Figure 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 DG Swell , attain plateaus for ½ATP ) ½ADP (see 4, top panel), whereas DG h increases monotonically with the nucleotide ratio. The regime where transduction is maximally efficient 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. In this two-step mechanism, the amount of energy available for the mechanical expansion is limited by that provided by non-equilibrium Hsp70 binding, which does not explicitly depend on the overall conformational properties of the substrate. To further dissect the energetic determinants of Hsp70 functioning and obtain more general conclusions, we thus analyzed the energy balance of chaperone binding to a model substrate, such as a peptide with a single binding site. To this aim, we focused on a simplified reaction cycle, which essentially corresponds to a single triangle within the overall scheme in Figure 3 and does not imply any conformational free-energy variation upon Hsp70 binding. We report in Figure 6   Source data 1. Molecular simulation data used to generate the plots in Figure 5.
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 orders of magnitude lower than its equilibrium counterpart, as already discussed in De Los Rios and . Here, we 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, DG h , as previously done in the case of DG Swell . Interestingly, also in this case the energy ratio is maximal in cellular conditions (red dashed line in Figure 6, bottom panel), suggesting that the optimality of the overall expansion process does not depend on specific features of the substrate but it is a direct consequence of the intrinsic kinetic parameters of Hsp70 chaperones.

Discussion
Integrating molecular simulations, polymer theory, single-molecule experimental data and non-equilibrium rate models, we have developed a comprehensive framework that provides a quantitative picture of Hsp70-induced expansion of substrate proteins and offers a broad insight into the cellular functioning of this versatile chaperone machine. We relied on molecular simulations for characterizing the structural and thermodynamic features of the complexes formed by the bacterial chaperone DnaK and its unfolded substrate rhodanese. Notably, we investigated a large variety of possible chaperone-substrate complexes for determining their conformational free energy as a function of stoichiometry and chaperone binding patterns. This computational strategy based on an enhanced-sampling protocol confirmed that excluded-volume interactions upon chaperone binding can greatly perturb the conformational ensemble of the unfolded substrate leading to its expansion. Remarkably, simulation results were found to be in excellent agreement with the predictions of Sanchez theory for globule to coil transition, thus providing another example of how polymer theory can be successfully used to decipher the behaviour of disordered proteins (Sherman and Haran, 2006;Hofmann et al., 2012;Schuler et al., 2016). We then combined conformational free energies with available biochemical data to develop an analytical rate model of the chaperone/substrate reaction cycle, which included both chaperone binding/unbinding and nucleotide hydrolysis/exchange processes.
This model fully takes into account non-equilibrium effects due to ATP hydrolysis and represents a natural extension of the ultra-affinity framework originally developed for peptide substrates with a single Hsp70 binding site (De Los Rios and . We could thus investigate the population of each complex and the average structural properties of the system as a function of the ATP/ ADP nucleotide ratio, which measures how far the system is from thermodynamic equilibrium. The reliability of the model was corroborated by a quantitative comparison with recent singlemolecule FRET data, indicating that our non-equilibrium framework accurately captures the salient features of the ATP-dependent expansion. We then used this unprecedented access to the thermodynamics details of this complex molecular process to compare the free-energy cost associated with substrate swelling with the chemical energy released by ATP-hydrolysis. Remarkably, this analysis revealed that energy transduction is maximally efficient for ATP/ADP values in cellular conditions. This result hints at the possibility that Hsp70 chaperones have been tuned by evolution to optimize the conversion of chemical energy into mechanical work for substrate expansion. Further analysis indicated that this optimality is likely inherited from the intrinsic properties of Hsp70 chaperones, which can convert up to 20% of the ATP chemical energy into non-equilibrium, excess binding energy at physiological conditions ( Figure 6, bottom panel).
From a broader perspective, the ATP-driven action of Hsp70s induces a non-equilibrium redistribution of their protein-substrates over their structural ensemble. In particular, thanks to the fine-tuning of the process by co-chaperones (J-domain proteins and Nucleotide Exchange Factors), the expansion process highlighted here, followed by substrate release, may result in the enhancement of the native state population beyond the predictions of thermodynamic equilibrium, as recently observed even under otherwise denaturing conditions (Goloubinoff et al., 2018). Consistently, Zhao and coworkers have recently observed that Hsp70 chaperones crucially contribute in vivo to the solubility and functionality of a sizeable fraction of the E. coli proteome that, in their absence, would instead spontaneously misfold and aggregate (Zhao et al., 2019). Remarkably, a similar effect has been observed in vitro for the GroEL chaperonin (Chakrabarti et al., 2017;Goloubinoff et al., 2018), hinting at the possibility that multiple chaperone families might reshape the equilibrium conformational distribution of proteins through energy-consuming processes. These considerations might have important consequences for our ability to translate results from in vitro experiments to the active cellular context (Bershtein et al., 2013). Likewise, they raise fundamental questions about the evolution of protein sequences: indeed, since chaperones are ubiquitous and very much conserved across the different kingdoms of life, their ability to favor native states might have partially relieved the selection pressure for strong equilibrium thermodynamic stability, thus allowing evolution to proceed faster and to be more tolerant for slightly destabilizing mutations, as suggested in Rutherford and Lindquist (1998); Tokuriki and Tawfik (2009).
Besides the unfolding of non-native substrates discussed in this work, Hsp70s are highly versatile machines that play a fundamental role in a variety of diverse cellular functions such as protein translocation, protein translation, and disassembly of protein complexes. All these processes share basic analogies from the mechanistic point of view: in all these cases, Hsp70 binding to flexible substrates in constrained environments requires the energy of ATP hydrolysis (ultra-affinity) and results in the generation of effective forces due to excluded volume effects (entropic pulling), which ultimately drive protein translocation into mitochondria (De Los Rios et al., 2006;Assenza et al., 2015), clathrin cage disassembly (Sousa et al., 2016) and/or prevention of ribosome stalling (Liu et al., 2013). Here, by detailing how energy flows from ATP hydrolysis to mechanical work due to entropic pulling, we have elucidated a general force-generating mechanism of Hsp70 chaperones. This mechanism does not rely on any power-stroke conformational change but it rather depends on the efficient conversion of ATP chemical energy into ultra-affinity.

Molecular model
In all the simulations, rhodanese and Hsp70 were coarse grained at the single-residue level as collections of beads centered on the C a atom of each amino acid. The unfolded state of bovine rhodanese (PDB:2RHS) was modeled according to the force field for disordered proteins from Smith et al. (2014). Two-and three-body bonded interactions along the substrate backbone were included via harmonic potentials, namely V bond ¼ k l P b r b À l ð Þ 2 =2 and V bend ¼ 1 2 k P a a À 0 ð Þ 2 , respectively. In the previous formulas, r b denotes bond lengths; a the bend angles; l ¼ 3:9 Å ; ðk l =k B TÞ À1=2 ¼ 0:046 Å ; 0 ¼ 2:12 rad; ðk =k B TÞ À 1 2 ¼ 0:26; and k B T is the thermal energy. Four-body bonded interactions were implemented as Fourier terms, , where f d is the torsion angle and A 1 ¼ 0:705, A 2 ¼ À0:313, A 3 ¼ À0:079, A 4 ¼ 0:041, B 1 ¼ À0:175, B 2 ¼ À0:093, B 3 ¼ 0:030, B 4 ¼ 0:030. The steric repulsion was implemented through a Weeks-Chandler-Andersen potential, In the previous formula, r ij is the distance between beads i and j, while s ¼ 4:8 Å . The hydrophobic part of the potential is specific to the interacting residues and is modeled as the attractive part of the Lennard-Jones potential, In the previous formula, h ¼ 0:7722 k B T sets the overall strength of the hydrophobic interactions, while ij depends on the residues i and j involved in the interaction, and is defined as the geometric mean of their hydrophobicities, ij ffiffiffiffiffiffi ffi i j p . The values of the hydrophobicities considered are based on a shifted and normalized Monera hydrophobicity scale (Smith et al., 2014). Electrostatic interactions were neglected based on control FRET experiments (see section 5.3). Without further tuning, this force field gives a radius of gyration of unbound rhodanese equal to R g ¼ ð23:3 AE 0:1Þ Å , which is in good agreement with the experimental value R g ¼ ð20:1 AE 0:8Þ Å (Kellner et al., 2014;Hofmann et al., 2014).
We modeled Hsp70 by means of a simple structure-based potential (Assenza et al., 2015) built on the conformation of ADP-bound Hsp70 (Bertelsen et al., 2009, PDB:2KHO). We described both the NBD (residues 4-680) and the SBD (residues 690-603) as rigid bodies whereas the interdomain linker (residues 681-689) was modeled according to the potential for flexible proteins described above. Importantly, non-bonded interactions of Hsp70 residues were limited to excluded volume effects and described by WCA potential (see Equation (5)). Electrostatic interactions were not explicitly included due to their marginal role in Hsp70/rhodanese complexes evidenced by FRET experiments (Appendix 1).
The binding sites for DnaK on the substrate were identified by applying the algorithms by Rüdiger et al. (1997) and Van Durme et al., 2009 on rhodanese and selecting only fragments for which at least partial consensus between the two predictions was obtained. Following this procedure, we identified six binding sites roughly centered on residues 10,118,131,162,188,260 of the rhodanese sequence. The residues of the binding site were aligned to a SBD-bound peptide reported in the literature (PDB: 1DKX Zhu et al., 1996) and constrained to move rigidly with the corresponding SBD, thus ensuring that each chaperone was irreversibly bound to the substrate. Following this procedure, 2 6 = 64 different chaperone/substrate complexes were built depending on the occupancy of each binding site.

Simulation protocols
All the simulations were performed with a version of LAMMPS (Plimpton, 1995) patched with the open-source, community-developed PLUMED library (Bonomi et al., 2019), version 2.1 (Tribello et al., 2014). The temperature T ¼ 293 K was controlled through a Langevin thermostat with damping parameter 16 ns -1 . The time step was set equal to 1 fs, and each residue had a mass equal to 1 Da.
In order to obtain conformational properties, for each of the 64 chaperone/rhodanese complexes we performed at least 10 independent simulations of 2 Á 10 7 timesteps. To ensure that full equilibration was achieved, only the last 10 7 timesteps of the obtained trajectories were considered for analysis. Statistical errors on the computed quantities were estimated as standard errors of the mean computed across independent realizations and are smaller than the size of symbols reported in the figures. The FRET efficiency E for a given couple of dyes was computed starting from the distance r separating the corresponding amino acids as where r 0 = 54 Å , as in Kellner et al. (2014). For each realization, the time average of E was computed. The final values employed to compute the results reported in the inset of Figure 5 in the main text were obtained as the average between independent realizations. The conformational free energies were computed by means of steered simulations, where for each complex rhodanese was pulled from equilibrium until an elongated conformation was obtained (Figure 2-figure supplement 1, top panel). Due to the large intermolecular distances, the effect of chaperones on the conformational properties of fully-stretched rhodanese is negligible, so that this state can be used as a reference to compute the free-energy differences between different chaperone/rhodanese complexes. The pulling was implemented by adding a harmonic potential acting on the radius of gyration R g of rhodanese. The equilibrium position of the harmonic trap was increased at a constant pulling speed v = 10 -5 Å /fs from the equilibrium value up to R g ¼ R fin g . For each chaperone/substrate complex, 100 independent pulling simulations were performed, starting from uncorrelated initial snapshots extracted from the equilibrium distribution. For each realization, the work W performed by the bias potential during the steering process was measured. The free-energy difference dG between the equilibrium starting point and the reference state (corresponding to R fin g ) was then computed via the Jarzynski (1997) equality: where < . . .> denotes statistical average. The error on dG was estimated according to the bootstrap method. The quantity DG considered in the main text was finally computed as DG ¼ dG 0 À dG, where dG 0 corresponds to the case of rhodanese alone (Figure 2-figure supplement 1, bottom panel). The uncertainty on DG was estimated by propagating the error bars on dG and is always smaller than the size of symbols. In order to enhance the robustness of the results, the final values reported in the main text were obtained as a further average over the values of R fin g within the range 260 Å R fin g 290 Å .

Rate model
For the kinetic model we consider a system in which each of the six binding sites can either be occupied by a chaperone in the ATP or ADP state, or it can be free, so that in total there are 3 6 = 729 possible configurations. The concentration c i of each state evolves in time according to a system of rate equations where k ij is the transition rate from state i to state j. The first term in the right hand side (r.h.s.) of Equation (9) represents the total flux of molecules from the other states toward state i, while the second term in the r.h.s. of Equation (9) accounts for the flux of molecules from state i to any other state. We focused on the steady-state, when the concentrations of the various states do not change over time, which is defined by Here, we provide a list of the relevant reactions that must be taken into account and of their corresponding rates. Each configuration is labelled by means of six symbols: 0 for empty sites, T for sites occupied by an ATP-bound chaperone and D for sites occupied by an ADP-bound chaperone (e.g. ð0; T; 0; D; 0; 0Þ, where the first, third, fifth and sixth Hsp70 binding sites are unoccupied, the second binding site is associated to a chaperone in the ATP-bound state while the fourth binding site is associated with a chaperone in the ADP-bound state). With this notation, the rates corresponding to every reaction are easily determined. Examples of the reactions that need to be considered are We further provide, as an example, the equation for a precise configuration, say ð0; T; 0; D; 0; 0Þ (here the label stands for the concentration of the configuration). The two binding sites that are occupied can undergo chaperone unbinding, ATP hydrolysis/synthesis or nucleotide exchange. The remaining unoccupied binding sites can bind either an ATP-or an ADP-bound chaperone. We thus have d dt ð0; T; 0; D; 0; 0Þ ¼ Àð0; T; 0; D; 0; 0Þ Ã ðk eff ex;DT þ k s s þ k eff ex;TD þ k s h þ k atp off þ k adp off Þþ þð0; 0; 0; D; 0; 0Þ½Hsp70 Á ATPk atp on e ÀbDDG þ þð0; T; 0; 0; 0; 0Þ½Hsp70 Á ADPk adp on e ÀbDDG þ þðT; T; 0; D; 0; 0Þk atp off þ þðD; T; 0; D; 0; 0Þk adp off þ þð0; T; T; D; 0; 0Þk atp off þ þð0; T; D; D; 0; 0Þk adp off þ þð0; T; 0; D; T; 0Þk atp off þ þð0; T; 0; D; D; 0Þk adp off þ þð0; T; 0; D; 0; TÞk atp off þ þð0; T; 0; D; 0; DÞk adp off : Below we further detail the rates of our model. It is possible to move from an ATP-state to an ADP-state either via hydrolysis/synthesis or via nucleotide exchange. In the case of exchange, effective constants are used, which take into account the unbinding of one nucleotide species and the binding of the different one. The effective exchange rates are thus a function of the ratio [ATP]/[ADP] (see also De Los Rios and : where k þD , k þT , k ÀD and k ÀT are the binding and unbinding rates for ADP and ATP respectively. The rates of binding between the chaperone and single peptides have been previously determined experimentally (Mayer et al., 2000), and they were corrected in order to take into account the conformational change of the full polypeptide substrate upon binding, as we illustrated in the main text.
Substrate binding enhances the chaperone ATPase activity. Furthermore, the stimulation of ATP hydrolysis always takes place in cooperation with JDP co-chaperones. In our model, we did not consider them explicitly but their contribution was implicitly included through the choice of the rate constants.
In particular, the hydrolysis rate in the absence of the substrate, k h , is much smaller than the same rate in the presence of the substrate, k s h (k h ( k s h ). We assumed that the ratio between the rate of hydrolysis k h and the rate of synthesis k s is not altered by the substrate: The substrate binding/unbinding rates, the rates of nucleotide exchange and the hydrolysis and synthesis rates are collectively constrained by thermodynamic relations. Indeed, when the ratio between the concentrations of ATP and ADP is equal its equilibrium value (when the spontaneous hydrolysis and synthesis reactions are at steady state and compensate each other), detailed balance must be satisfied (Ge et al., 2012). As a consequence, for every closed cycle in the reaction network the product of the rates in one direction must be equal to the product of the rates in the opposite direction. Therefore, if k atp on , k adp on , k atp off and k adp off are the rate of substrate binding and unbinding from a chaperone in the ATP and ADP states, we must have k atp on k s h k adp off k s k adp on k s s k atp off k h ¼ k atp on k adp off k adp on k atp off ¼ 1 : Remarkably, taking the rates as provided in Mayer et al. (2000); Hu et al. (2006); Kellner et al. (2014), this relation is not satisfied, and we had thus to modify them. We thus calculated the product in the formula above and then corrected the rates in the following way: k atp on k adp off k adp on k atp off ¼ r k atp on ; k adp off ! k atp on =r 1=4 ; k adp off =r 1=4 (17) k adp on ; k atp off ! k adp on Ã r 1=4 ; k atp off Ã r 1=4 : The concentration of free chaperones in the ATP and in the ADP states was obtained, at the leading order, by solving a three-state system whose reactions have the form Hsp70 þ ADP ÀÀ* )ÀÀ Hsp70 Á ADP ÀÀ* )ÀÀ Hsp70 Á ATP ÀÀ* )ÀÀ Hsp70 þ ATP: Since we worked in the assumption of excess of chaperones in the system, once these concentrations were obtained, they remained fixed once for all, without being considered as a variable of the biochemical network.
We report in the following table the rates used in the model. Parameters of the model (Mayer et al., 2000;Hu et al., 2006;Kellner et al., 2014) k atp off 2:31 s À1 k adp off 2 Ã 10 À3 s À1 k ÀT 1:33 Ã 10 À4 s À1 k ÀD 0:022 s À1 k atp on 1:28 Ã 10 6 M À1 s À1 k adp on 10 3 M À1 s À1 k þT 1:3 Ã 10 5 M À1 s À1 k þD 2:67 Ã 10 5 M À1 s À1 k h 6 Ã 10 À4 s À1 k s h 1:8 s À1 To test the robustness of the model for the radius of gyration, the average FRET efficiency and the free-energy DG swell , 100 realizations were implemented, taking each time the values DG i from a Gaussian distribution with s ¼ 0:3 kcal/mol and mean equal to the value obtained with the MD simulations. Figures 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 (Pettersen et al., 2004).