ATP regulates the function of many proteins in the cell by transducing its binding and hydrolysis energies into protein conformational changes by mechanisms which are challenging to identify at the atomic scale. Based on molecular dynamics (MD) simulations, a method is proposed to analyze the structural changes induced by ATP binding to a protein by computing the effective free-energy landscape (FEL) of a subset of its coordinates along its amino-acid sequence. The method is applied to characterize the mechanism by which the binding of ATP to the nucleotide-binding domain (NBD) of Hsp70 propagates a signal to its substrate-binding domain (SBD). Unbiased MD simulations were performed for Hsp70-DnaK chaperone in nucleotide-free, ADP-bound and ATP-bound states. The simulations revealed that the SBD does not interact with the NBD for DnaK in its nucleotide-free and ADP-bound states whereas the docking of the SBD was found in the ATP-bound state. The docked state induced by ATP binding found in MD is an intermediate state between the initial nucleotide-free and final ATP-bound states of Hsp70. The analysis of the FEL projected along the amino-acid sequence permitted to identify a subset of 27 protein internal coordinates corresponding to a network of 91 key residues involved in the conformational change induced by ATP binding. Among the 91 residues, 26 are identified for the first time, whereas the others were shown relevant for the allosteric communication of Hsp70 s in several experiments and bioinformatics analysis. The FEL analysis revealed also the origin of the ATP-induced structural modifications of the SBD recently measured by Electron Paramagnetic Resonance. The pathway between the nucleotide-free and the intermediate state of DnaK was extracted by applying principal component analysis to the subset of internal coordinates describing the transition. The methodology proposed is general and could be applied to analyze allosteric communication in other proteins.
Biomolecular machines, as motor proteins [1] and ATP-dependent molecular chaperones [2], undergo a major conformational change upon ATP binding and hydrolysis. In many cases, the ATP molecule regulates the protein function by an allosteric mechanism [3]–[7]: binding of ATP in the ATPase domain of these proteins transmits a conformational change to a distant binding site where the substrate binding propagates an analogous conformational change towards the nucleotide-binding site. The exact biophysical characterization of the mechanisms of propagation of a conformational change induced by ligand binding through a protein remains a major challenge in structural biology [8], [9]. In biomolecular machines, as the ATP-dependent molecular chaperones [10], crystallographic data of the different nucleotide-bound states do not provide generally a complete description of the mechanism by which the protein undergoes a conformational change, because the pathway between the different conformational states and the dynamical information are missing [5]. In recent years, NMR, among other techniques, contributed to the characterization of allostery and protein dynamics [8], [9], [11] including nucleotide-dependent conformational changes in molecular chaperones [12], [13]. Conformational changes induced by a ligand can be mediated by perturbations of the mean protein conformation (enthalpic effect) as well as by changes of its fluctuations and dynamics (entropic effect) [5], [9], [14]–[16]. Both enthalpic and entropic effects are encompassed in the free-energy landscape (FEL) of a protein [17]. In solution, the FEL is best regarded as a multi-dimensional surface with multiple local minima separated by barriers [17]. Different protein-bound states correspond to different FELs, with different shapes and distributions of populations of conformers [4], [5]. Hence, the nucleotide binding and the nucleotide release in a biomolecular machine redistribute its conformational substates and modify its FEL [18]. In many cases, how the dynamics contributes to the conformational transition at the different time-scales of the protein motions and how the FEL is modified upon nucleotide binding remain elusive. Nowadays, standard all-atom molecular dynamics (MD) simulations [19]–[21] represent a powerful tool, complementary to experiments, for investigating the conformational sampling and nucleotide-dependent conformations of biomolecular machines [22]–[25]. The fully unbiased all-atom MD simulations permit to monitor the sequence of spontaneous conformational changes at the atomic level [26] and to translate the conformational fluctuations in terms of modified protein FEL. Although unbiased MD simulations of a fully solvated protein may extend to the millisecond range [27] with specific hardware [28], the current state-of-art of unbiased MD simulations correspond typically to 100 ns up to 1 µs time-scale for a system of several thousands of atoms in explicit solvent with the standard computational capabilities available to most of the research groups [29]–[31]. However, transition between active and inactive states in proteins occurs typically on the microsecond to the millisecond time-scales [32]. The complete conformational transition between different nucleotide-binding states of a molecular chaperone in fully unconstrained all-atom MD is thus still inaccessible [31]. To overcome this time-scale restriction, different computational strategies were developed to shorten artificially the time-scale of the conformational transitions. Targeted MD [33], accelerated MD [34], metadynamics [35] to cite a few allow the acceleration of the conformational sampling of all-atom MD by adding a bias potential. An alternative strategy consists to eliminate the fast degrees of freedom of the protein from the MD calculation by coarse-graining the protein structure. A recent application to the 70 kDa heat-shock protein (Hsp70), a key ubiquitous molecular chaperone [36], allow us to observe a short-life time full transition between an ATP-like state and an ADP-like state using a coarse-grained model where the solvent and the nucleotides were implicit [37]. However, current coarse-grained models of proteins lack the description of the crucial explicit interactions between the nucleotide and the protein which can be only described accurately at the atomic level. In spite the fact that all-atom unbiased MD simulations are not able to follow the full conformational change of a large ATP-dependent protein, as Hsp70 analyzed in the present work, MD simulations probing the structural fluctuations around the end points of a conformational transition are useful to understand the conformational changes. Normal mode analysis (NMA) has been also extensively applied to study the conformational transitions in proteins because the directions of the low-frequency vibrational modes provide the direction of the largest deformation at thermal equilibrium and can serve as collective coordinates to describe the conformational changes [38], [39]. Indeed, intrinsic low-frequency vibrational dynamics of proteins in different conformational states correlates with the structural changes between the different conformations induced by ligand binding or by protein-protein interactions [16], [40]–[44]. Aside the limited sampling issue of standard all-atom MD simulations of large conformational changes upon ligand binding, another issue is to extract from these data relevant information about the mechanism of propagation of a small perturbation at one ligand binding site to a different remote binding site. To get insight into the mechanisms, several analytical tools were applied to ATP-dependent conformational transitions such as the analysis of contact maps [45] and of the correlations between the residue motions [22]. In NMA, relevant information is gained by projecting the direction of the low-frequency modes of a given conformation on the shortest path between the two end structures of a conformational transition [1], [23], [46]–[49]. None of these approaches is easily related simultaneously on the FEL of the protein and to mutagenesis experimental data. In the present work, we addressed these two issues (MD sampling and relation between FEL and biochemical data) for the conformational change induced upon nucleotide binding to a molecular chaperone Hsp70 by using rather extensive unbiased MD simulations of the protein in three different nucleotide-binding states; nucleotide-free, ADP-bound and ATP-bound. Hsp70 is a key molecular chaperone which assists in the correct folding and refolding of proteins in all prokaryotic and eukaryotic organisms [2], [36], [50]–[53]. Hsp70 s contribute to other crucial cellular processes including the protein degradation, the translocation of peptides across the cell membrane, the formation of protein complexes, and the inhibition of the programmed cell death (apoptosis) [54]–[57]. Among the Hsp70 s, human Hsp70 (hHsp70) has attracted great interest because of its demonstrated implications in numerous misfolding diseases [58] (Parkinson, Alzheimer…) and in cancer [59]. Hsp70 s of all species share common structural features and it is hypothesized that they perform their biological functions in a similar manner. Most mechanistic information about hHsp70 was derived from comparisons with the highly homologous bacterial E. coli Hsp70 (DnaK) (see Table S1). It should be mentioned however that a high homology of sequence between two proteins does not guarantee that their allosteric mechanisms and residues involved are similar [60]. In the present work, the calculations and analysis were performed for DnaK because a large set of structural and biochemical data are available (see Results and Discussion sections for a comparison with the present results and their implications for hHsp70). In particular, the allosteric communication in Hsp70 s, i. e. how the binding of a ligand in the nucleotide-binding domain (NBD) modifies the conformation of the remote substrate-binding domain (SBD) of the protein, was extensively studied for DnaK [12], [61]–[69]. The structural model used in the present work is the two-domain structure of ADP-bound DnaK (residues 4–603, missing the C-terminal part) solved by NMR (PDB ID: 2KHO) [70]. Our first aim was to probe the conformational dynamics and the structural perturbations induced upon exchange of ADP by ATP and by the removal of ADP from this NMR model of E. coli DnaK on the micro-second time-scale at the atomic scale. The sampling issue of the unbiased MD simulations of E. coli DnaK was addressed by considering two trajectories with different initial conditions of the protein in each of its nucleotide-binding states and more important by restricting the data to a subset of the conformational degrees of freedom for which the convergence between the two trajectories was reached. This was achieved by projecting the FEL of the protein along its amino-acid sequence by computing the free-energy profiles (FEPs) along the coarse-grained dihedral angles (CGDAs) γ defined from four consecutive Cα atoms of the protein (Fig. S1) and describing the backbone motions [71], [72]. Only the FEPs which were similar in the two trajectories of the protein in the same nucleotide-binding state were kept for the analysis of the allosteric communication, we named these FEPs the DATA. No information can be exploited from the other FEPs: they may or not be relevant to the conformational fluctuations of the protein (NO DATA set). With this approach, we were able to quantify how each FEP varies along the amino-acid sequence depending on the nucleotide-binding state of the protein and to relate the FEP modifications to biochemical data of wild-type or mutant E. coli DnaK chaperones. Indeed, the comparison between the DATA sets of FEPs for E. coli DnaK in different nucleotide-binding states lead finally to a subset of only 27 key dihedral angles involved in the conformational nucleotide-dependent transitions. The fluctuations of these 27 CGDAs γ on their FEPs are subdiffusive [71]–[73]. The coupling between the CGDAs γ can be revealed by analyzing the correlation of their motions using the dihedral principal component analysis (dPCA) [74] allowing the definition of collective motions to which these degrees of freedom contribute the most, in a similar manner to principal component analysis in Cartesian coordinates [75], [76] or to the full-correlation analysis based on the information theory [77]. The study of the FELs was complemented by an analysis of the long-life time contacts between residues present in ATP-DnaK trajectories and not present in the trajectories of the protein in the other nucleotide-binding states. The paper is organized as follows. Results are presented in the next section and are compared to available experimental data. The structure of Hsp70 s and the results of the unbiased all-atom MD simulations for E. coli DnaK are shown. There, we described and compared to available experimental data an intermediate state found by MD upon ATP-binding, named ATP*-Hsp70, where the two domains of the protein are docked on a sub-microsecond time-scale. The construction of a map of persistent residue contacts in ATP-bound DnaK conformational dynamics is described and compared to available experimental data. Next, the methodology to cope with the limited sampling of the present unbiased MD simulations and to extract the key residues from the analysis of the FEPs along the amino-acid sequence is presented, as well as its application to the MD simulations of ATP-bound DnaK. Then, a comparison between the key residues extracted from this new methodology with biochemical and structural data on the allostery in Hsp70 s is discussed in detail. The paper ends with a short conclusion where the major results are summarized and with the Methods section.
In order to characterize the interdomain interactions leading to the ATP*-DnaK intermediate state (Fig. 2 C), as well as the interactions between the nucleotide and the NBD of E. coli DnaK, the contact maps (CMs) of ADP and ATP-bound DnaK and the contact map difference (CMD) ATP/ADP were computed for the MD trajectories of E. coli DnaK in these two different nucleotide-binding states (see Methods section for the calculations of CM and CMD). The CM of ATP-bound DnaK (Fig. S2 A), which is composed of 1994 contacts, contains only 6 interdomain contacts (contacts between NBD/linker, linker/SBD or NBD/SBD): R167-L390, R167-L391, L391-T417, L391-P419, L392-T416 and D393-N415 (Fig. 4 A). The CM of ADP-bound DnaK has no interdomain contacts (Fig. S2 B), which means that the interdomain contacts are unique to the ATP* intermediate state, as revealed by the CMD ATP/ADP (Fig. S2 C). In MD, this major difference about the interdomain contacts in ADP-bound and in ATP-bound DnaK arises from a different interaction between the NBD and the linker in these two nucleotide-binding states, as shown by the different distances R167-L390 and R167-L391 in these two states (Table 1). The comparison between MD results for ADP-bound and ATP-bound DnaK shows that the communication between the NBD, the linker and the SBD involves the residue R167 of the subdomain IA of the NBD, the residues L390 and L391 of the linker and the residues 415–419 of the subdomain β of the SBD. These residues were identified to be crucial for the interdomain communication in DnaK. Indeed, the mutation R167A induces a loss of allosteric coupling between the NBD and the SBD [64] and the mutations K414I [90] and N415G [91] eliminate and attenuate allostery, respectively. The residues L390, L391 and L392, which belong to the VLLL motif of the linker (highly conserved between the different species of Hsp70 [67]), was shown essential for E. coli DnaK' s in vivo functions [67]. In addition, it has been proposed that ATP binding to the ATPase domain causes a bending of the C-terminal part of the linker of E. coli DnaK at D393, in order to facilitate an insertion of the linker between the NBD and the SBD [64]. In the present MD simulations, the coupling between the NBD and the SBD confirms the insert of the linker upon ATP binding between the two domains through the interactions R167/L390 and D393/N415. The contacts between ATP and the NBD were identified (Fig. 4 B and Fig. 4 C). This network of contacts is formed by 12 residues which belong to the subdomains IA, IIA and IIB of the NBD: T11, T12 and N13 (IA); G195, G196, G197, G198 and T199 (IIA); I271 (IIB); G341, G342 and R345 (IIA). Among these 12 residues, the CMD ATP/ADP reveals that only 6 are unique to the ATP-bound state of E. coli DnaK and were not found in the APO and ADP-DnaK MD simulations performed in the present work (Fig. 4 B and C and Fig. S2 C): T11, G197, G198, T199, I271 and G341. The contacts depicted in Fig. 4 for the intermediate state ATP*-DnaK can be compared qualitatively with those observed in the very recent crystallized triple mutant ATP-bound DnaK [69]. Although the structure of the intermediate state found in MD and the crystallized open conformation of DnaK are quite different as discussed above, many contacts between ATP and the NBD found in ATP*-DnaK (Fig. 4 B and C) are also found in the experimental structure of ATP-bound DnaK (PDB ID: 4B9Q) [69]. Except D195-ATP, at a contact distance (see Methods section) of 4. 9 Å in the open X-ray structure, all contacts found in MD between the NBD and ATP (Fig. 4 B) occur also in the ATP-bound experimental structure. However, contacts between ATP and the NBD observed in the fully open structure are missing in ATP* intermediate state: in particular the residue K70, which is at 2. 8 Å from ATPγ in the experimental structure, is at about 4. 2 Å in ATP*-DnaK. A contact occurs however between K70 and E171 both in the ATP*-DnaK MD trajectory (average distance K70-E171 = 2. 4 Å) and in the experimental open conformation but not in the ADP-DnaK trajectories (average distance ADP K70-E171 = 8. 6 Å). Because the metastable position of the linker is different in ATP* compared with the position of the linker in the open conformation of DnaK (Fig. 1B), different interdomain contacts were found experimentally in particular R167-D481, I168-D481, N415-D326 and D393-N170. The analysis of the contact maps provides only a static representation of the MD data and of the influence of ATP binding on the conformations of E. coli DnaK. The modifications of the protein dynamics induced by the nucleotide and the correlations between the different structural modifications along the DnaK amino-acid sequence are not easily interpreted from this sole analysis. As stated in the introduction, the FEL holds the key to understand how a nucleotide modifies the properties and populations of the protein conformers. To provide the link between the FEL and the global conformational change observed in ATP-bound MD trajectories, we analyze next the free-energy profiles (FEPs) of the CGDAs γ (see Fig. S1 and Eq. 1) along the amino-acid sequence of E. coli DnaK. In order to identify which residues are involved in the nucleotide signal propagation along the backbone of E. coli DnaK, we examined which of the 597 CGDAs γ have a FEP depending on the nucleotide-binding state of the NBD. However, the FEPs of DnaK in a given nucleotide-binding state computed from the present unbiased MD trajectories are effective FEPs. The effective FEP differs from the actual FEP, which is an equilibrium thermodynamic property and should be computed from a large ensemble of trajectories of the protein in the same nucleotide state. However, the size of the system simulated (about 500000 atoms) and our computing capabilities limited the number of trajectories to only two trajectories of several hundred of nanoseconds per nucleotide-binding state of the NBD. Because of this limited sampling issue, only the FEPs which are similar in different trajectories in a given nucleotide-binding state were considered.
In summary, all-atom MD simulations of Hsp70-DnaK in different nucleotide-binding states in explicit solvent were performed. We observed the docking of the SBD of DnaK onto the lobe I of its NBD upon ATP-binding: the structure found is an intermediate plausible ATP-bound state, named ATP*, which is compared to several experimental data. A strategy was presented in which the FEL computed from unbiased MD simulations is represented by FEPs of CGDAs of the main chain along the amino-acid sequence of the protein. The FEPs can be quantitatively compared in different nucleotide-binding states by using a similarity index. The analysis of the FEPs allowed us to identify a small network of 27 CGDAs γ. The coupling between these 27 CGDAs was deciphered by using dihedral principal component analysis. The conformational change induced by ATP binding was represented by a pathway of minimum energy on the FES built from the two lowest dihedral principal components of the ATP-bound DnaK trajectories. The 27 coordinates correspond to 91 residues which are involved in the interdomain communication upon nucleotide binding to DnaK. Most of these 91 residues revealed by MD were shown experimentally to be relevant for the communication between the NBD and the SBD in the Hsp70 chaperone cycle. The present work also suggests 26 new key residues that could be tested experimentally. For example, the residues A503 to G506, defining γ504, belonging to the kink between SBD-α and SBD-β, seem to play a key role in the dynamics of the docking as well as for the intrinsic dynamic of the SBD. We found that the conformation of these residues influences the distribution of the distance E430/R547 explaining the structural origin of the variations of this distance observed in EPR upon ATP binding to E. coli DnaK [81]. The strategy developed to analyze the long-range effect of ATP binding on the structure and dynamics of E. coli DnaK is general and could be applied to decipher allosteric communication in other proteins.
All unbiased all-atom MD simulations in explicit water of E. coli DnaK using the NMR-derived structure (PDB ID: 2KHO in which the missing atoms of ARG and TYR were added) [70] were carried out with the GROMACS software package [105] and the GROMOS96 ffG43a1 force field [106], [107]. The time step used in all simulations was 0. 001 ps and the list of neighbors was updated every 0. 005 ps with the ‘grid’ method and a cut-off radius of 1 nm. The coordinates of all the atoms in the simulation box were saved every 2 ps. The initial velocities were chosen randomly. We used the NPT ensemble with a cubic box. The temperature and pressure were kept to the desired value by using the Berendsen method [108] and an isotropic coupling for the pressure (T = 300 K, τT = 0. 1 ps; P0 = 1 bar, coupling time τP = 0. 5 ps). The electrostatic term was computed by using the particle mesh Ewald algorithm [109] (with a radius of 1 nm) with the Fast Fourier Transform optimization (on 140 points for DnaK, with an order equal to 4 for the interpolation). The “cut off” algorithm was applied for the non-coulomb potentials with a radius of 1 nm. The system was warmed up for 40 ps and equilibrated for 600 ps with lower restraints, finishing with no restraints at 300 K. We performed 6 MD runs, 2 in each nucleotide-binding state (APO1, APO2, ADP1, ADP2 and ATP1, ATP2), as described in the following subsections. The stability of the protein structure in each MD trajectory was monitored by computing the Root-Mean-Square-Deviation (RMSD) between the molecular structure every 2 ps and the initial structure (as done before [23], Fig. S10 A). To reduce the computational cost, the APO and ADP-bound DnaK MD runs were stopped when the RMSD was converged to a constant value with a standard deviation lower than 0. 5 Å (Fig. S10 B) during at least 50 ns, resulting in six MD runs of different durations. The 2 MD runs of ATP-DnaK were continued after convergence of the RMSD up to 0. 5 µs in order to explore a possible opening of the SBD which did not occur however on that time-scale. The conformational fluctuations of the main chain of E. coli DnaK were monitored by recording the fluctuations of the CGDAs γ [71]. The CGDA γ for a residue i is the dihedral angle formed by the vectors (virtual bonds) joining 4 successive Cα atoms (i−1, i, i+1 and i+2) along the amino-acid sequence (Fig. S1 A) [117]. The first dihedral angle along the sequence is γ2 and the last is γN-2, where N is the total number of residues. Each CGDA γi is a local probe of the FEL along the amino-acid sequence of the protein [71]. The FEP of each γi is computed from the usual Boltzmann formula: (1) where kB is the Boltzmann constant, T the temperature (300 K) and p (γi) is the probability distribution function (PDF) of the CGDA γi computed from the MD trajectory (Fig. S1 B) [71]. Note that the experimental NMR structure of E. coli DnaK is comprised of 600 residues, starting from residue 4 to residue 603 [70], and is described by 597 angles γi, from index i = 5 to index i = 601. To quantify the convergence between two FEPs of a same CGDA γi in two MD trajectories with the same nucleotide-binding state but with different initial conditions, we computed the similarity index H between their associated PDFs [72]: (2) where p1 (γi) and p2 (γi) are two PDFs of the same CGDA γi computed from two MD trajectories 1 and 2. The similarity index H varies between 0 (dissimilar) and 1 (identical). Two FEPs, v1 (γi) from the MD trajectory 1 and v2 (γi) from the MD trajectory 2, are considered similar (“converged”) if the index H is larger than 0. 5. Distances between all pairs of atoms of the 600 residues were calculated every 2 ps for each MD run (APO1, APO2, ADP1, ADP2, ATP1 and ATP2). The distance matrix between the residues was built by retaining only the smallest atomic distance between the atoms of each residue pair. If this smallest distance was smaller than 4 Å, the residues were said in contact (contact matrix element is 1) and not in contact otherwise (contact matrix element is 0). For each nucleotide-binding state (APO, ADP and ATP), only contacts with at least 50% of occupancy in the two MD runs of E. coli DnaK in the same nucleotide-binding state (for example ATP-bound DnaK) and for which the average distance was similar in the two MD runs (for example MD runs ATP1 and ATP2) within 0. 5 Å were finally considered, defining the contact map (CM, Fig. S2 A and B) of a nucleotide-binding state (for example ATP). In addition, the difference between the contact maps, the contact map difference (CMD, Fig. S2 C) was computed. For example, the CMD ATP/ADP is the difference between the contact matrices ATP and ADP defined above, and thus reveals the atomic contacts present in the ATP nucleotide-binding state of DnaK, but not in its ADP nucleotide-binding state. In dPCA [74], the correlated internal motions of a protein are quantified by a covariance matrix Cij [Eq. (3) ] of the Cartesian components of the two-dimensional vectors u (t) [Eq. (4) ] representing the CGDAs: (3) where <…> denotes the average over all sampled conformations and (4) The diagonalization of the covariance matrix Cij gives 2N (= number of CGDAs γ eigenvalues λk, ordered by decreasing value: λ1>λ2>…>λ2 and eigenvectors ek = {e1k, e2k, …, e27k}T, with eik = {eik (X); eik (Y} where eik (X) is the Cartesian component x of the projection of the eigenvector ek on the ith CGDA (similarly for eik (Y) ). The contribution of each CGDA γi to a mode k is the so-called influence [74]: (5) The mean-square fluctuation (MSFi) of each vector ui can be decomposed into collective modes: (6) The eigenmodes with the largest eigenvalues λk (named slow modes) correspond to the collective modes contributing the most to the MSF of the protein [see Eq. (6) ]. Therefore, λ1 is the eigenvalue of the collective mode which contributes the most to the structural fluctuations of the CGDAs γi in an MD run [72], [74]. Finally, the projection of the MD trajectory on the kth eigenvector ek is named the dihedral principal component dPCk, (7) More details can be found in Refs. [72], [74], [118], [119]. The effective two-dimensional FES shown in Fig. 8 was computed from the two-dimensional probability density of dPC1 and dPC2 computed from the covariance matrix of the subset of the 27 dihedral angles involved in the conformational change of DnaK (all shown in Fig. 7). The two-dimensional FES computed from the covariance matrix of all dihedral angles is indeed irrelevant due to a lack of convergence sampling [94]. Indeed, some dihedral angles are in a diffusive regime. As shown in Refs. [93], [94], the diffusive regime of the internal degrees of freedom of a protein can be measured by the cosine content (CC) of the principal components. A large value (close to 1. 0) of the CC means that some of the degrees of freedom are fluctuating as the variables of free random walks on the time-scale of the present MD simulations. The FES built on principal components with CC close to 1. 0 shows artifacts, i. e. basins which are signatures of the cosine shape of the principal components and not of an actual minimum of the FES [93], [94], [119]. We computed the cosine content of each principal component [92]. The values of the cosine content computed for dPC1 and dPC2 calculated from the covariance matrix in the space of the 27 dihedral angles having a converged FEP are very small (CC1 = 0. 11 and CC2 = 0. 16 for ATP1 trajectory for example). On the contrary, the cosine content of the principal component dPC1 computed in the space of all dihedral angles is very large (CC1 = 0. 69 for ATP1 trajectory for example). The Fig. S11 A shows an example of a free-energy surface (FES) generated by projecting a MD trajectory on the essential plane defined by two principal components [Eq. (7) ]. In this example, the FES has two minima A1 and A2. The pathway that connects the two minima A1 and A2 by traveling along the “valleys” of the FES and going over the “passes” (which are saddle-points on the surface and correspond to the transition states) in order to minimize the free-energy difference while moving, is commonly referred to the Pathway of Minimum Energy (PME) [120]. To find the PME for the transition A1→A2, we applied the following methodology. First, an initial guess of the path is made, here the straight interpolation line d1 from A1 to A2. Secondly, we move from A1 along the direction d1 by 0. 1 (which is the bin size defined for the calculation of FES) and draw a segment p1 perpendicular to the direction d1. Note that the size of the segment p1 is chosen to belong to a circle of diameter A1A2. Finally, we search along the segment p1 the bin with the minimum of free-energy, named s1. This bin is also considered as the first step of the PME. We repeat the algorithm by constructing a straight interpolation line d2 from s1 to A2. The algorithm ends after N steps when the distance between sN and A2 is smaller than the size of the bin defined for the calculation of FES. The PME corresponds to the ensemble {sj}j = 1, . . , N and to the corresponding free-energy profile VPME (sj) (Fig. S11 B).
The precise biophysical characterization of the mechanisms of the protein conformational changes controlled by a nucleotide remains a challenge in biology. Molecular dynamics simulations of proteins in different nucleotide-binding states contain information on the nucleotide-dependent conformational dynamics. However, it is difficult to extract relevant information about the conformation-induced mechanism from the raw molecular dynamics data. Herein, we addressed this issue for the major ATP-dependent molecular chaperones Hsp70 s, which contribute to crucial cellular processes and are involved in several neurodegenerative diseases and in cancer. To function, Hsp70 undergoes several conformational changes controlled by the state of its nucleotide-binding domain. We demonstrated that the analysis of the effective free-energy landscape of the protein projected along the amino-acid sequence and computed from the molecular dynamics simulations of Hsp70 in different nucleotide-binding states, holds the key to identify the key residues of the conformational induced pathway. Identification of the key residues involved in the propagation of the structural changes induced by ATP binding offer alternative druggable specific sites other than the ligand binding clefts. The methodology developed for Hsp70 is general and can be adapted to any ligand induced conformational change in proteins.
Decipher the Mechanisms of Protein Conformational Changes Induced by Nucleotide Binding through Free-Energy Landscape Analysis: ATP Binding to Hsp70