Large-scale motions of biomolecules involve linear elastic deformations along low-frequency normal modes, …

• Abstract
• Methods
• Results
• Discussion
• Conclusion
• Acknowledgements
• Footnotes
• References
• Figures

The Linear Elastic Model. We first study the deformation of the kinase starting from the open structure. The closed structure was fitted to the open structure, giving a rmsd between the two structures of 7.2 Å. The low-frequency modes were determined by using a coarse-grained model, the Tirion potential (19). Among the normal modes for the initial structure, normal mode 1 (the lowest-frequency mode) is most relevant, having the highest overlap to the conformational change. By generating some deformed structures after mode 1, the energy surface of the open form is evaluated (Fig. 1). If the Tirion potential were perfectly harmonic, this would give a quadratic function of the displacement from the open form with a curvature dependent on the frequency of normal mode frequency. The Tirion energy surface and this ideal harmonic description only agree in the vicinity of the energy minimum structure. The Tirion energies are larger than expected from the standard normal mode analysis, because a single normal mode does not accurately represent the conformational changes during the deformation. Examining the open and closed forms suggests the conformational change resembles the bending of a rod. Obviously, when bending a rod through a finite angle, a single normal mode only indicates the initial direction of the transformation. Bending is kinematically nonlinear in a Cartesian basis. This kinematic nonlinearity, well known in macroscopic elastic theory, is the main shortcoming of the strictly linear model.

Nonlinear Elastic Models. To extend the model to include the bending nonlinearity, we explore a nonlinear path by using combinations of normal modes that adiabatically change as we move from one global structure to the other as a basis. We used the combinations of one or three normal modes that are most relevant to conformational change (high overlap modes). The procedure to find the best combination of modes is iterative (see Movie 1, which is published as supporting information on the PNAS web site, www.pnas.org). For each structure on the pathway, the strain energy is first calculated by using the Tirion potential with a spring network defined from the open form structure. Then the spring network of connectivities is redefined at each step of the calculation before determining a new set of normal modes for the distorted structure. Following the energy along the resulting nonlinear path made by using only one mode gives much lower energy than when the energy was strictly computed along the linear path moving along the first low-frequency normal mode (Fig. 1). The resulting adiabatic approximation surface agrees well with the quadratic approximation. Thus, although the harmonic approximation as a function of Cartesian coordinates is not accurate, the harmonic approximation is quite good when the energy is written as a function of a nonlinear reaction coordinate.

Strain Energy Is Localized. The strain in the molecule when it is deformed along the reaction path is not uniform. We assign a strain energy to each protein residue as the sum of the atomic strain energies. The strain intensity for all residues is determined for each structure along the nonlinear conformational change. Fig. 2a shows the spatial pattern of residue strain as the molecule is deformed. When the rmsd to the open form becomes >3 Å, some residues become particularly highly strained. These must rearrange to allow the transition to come to completion. The residue strain energies as a function of the sequence are shown in Fig. 2b. Fig. 2c shows the 3D structure of the protein with the residues colored according to their strain energies; red corresponds to high strain and blue corresponds to low (see Movie 2, which is published as supporting information on the PNAS web site). The T1 region, residues 10–12, is under high strain. The helix 6, residues 110–120, and helix a7, residues 160–170, are also under high strain. There is a clear correlation between the high-strain energy regions along the nonlinear conformational conversion path and what have been called hinges (15, 29).

The Cracking Model. The observations that some particular contiguous residues are under high strain leads us to hypothesize that these special regions may crack, i.e., unfold partially, during the conformational change. To include this possibility, a residue will be allowed to unfold if the strain energy of residue is higher than the difference of local folding energy and the structural entropy to be gained by locally unfolding. Under this local approximation, the free energy per residue can be written as

For simplicity, we assume that all residues have identical folding energy, E_o^residue, at the initial state and that their entropy change by unfolding, S^residue, is also identical. Experimentally determined stabilities are used to set these parameters. It is easy to incorporate residue specific energies for these thresholds as in the work of Munoz and Eason (10). The unfolding free energy of the adenylate kinase from Escherichia coli has been reported to be »4 kcal/mol for guanidine hydrochloride denaturation (30) and »10 kcal/mol for urea denaturation (31). Heat denaturation of adenylate kinase from Saccharomyces cerevisiae is reported to give a stability of »4 kcal/mol (32). Dividing the total unfolding free energy by the number of residues (»200) yields DG^residue = - TS^residue - E_o^residue  of »0.02–0.05 kcal/mol. This simple model of dividing the unfolding energy evenly neglects the fact that a good share of the binding in the folded state comes from nonadditive interactions such as hydrophobicity and side-chain conformations that enhance the cooperative aspect of denaturation. The unfolding energy for a segment should therefore be larger than the equilibrium value spread uniformly throughout the chain. Thus we examine the result by using two values of the threshold, one with and the other with 0.1 kcal/mol for unfolding. The cooperativity coming from nonadditive interactions also requires that several contiguous residues must be simultaneously under high strain to unfold (33). Thus, as in unfolding calculations (10), we require that at least a number N_min of contiguous residues should have strain energies exceeding the threshold, where , to be treated as being cracked. Models of folding kinetics fit experiment well by using N_min = 5 (10).

Fig. 3 shows the strain energy along the nonlinear conformational change path when cracking is allowed. Obviously partial unfolding gives a lower energy barrier. The extent of cracking strongly depends on the threshold . Cracking leads to a critical yield stress so the energy becomes an approximately linear function of deformation rather than following the quadratic dependence of the elastic regime. Such linear relations have also been invoked to explain the efficiency of motor proteins by Bustamante et al. (34). The cracking model suggests the linear behavior may be much more general. The critical yield stress found here, 140 pN, is comparable to values seen in single-molecule pulling experiments. As it becomes clear in the specific analysis exhibited later in this article, cracking leads to linear behavior of the rate versus driving force over a large range of reaction driving forces, a situation that would not be true if only elastic deformations were allowed. Studies on allosteric motions in hemoglobin by Eaton et al. (35) have also exhibited such a linear free energy dependence supporting the general applicability of the cracking mechanism.

Our theoretical analysis predicts that for adenylate kinase, the structure that corresponds to the beginning of cracking (at rmsd = 4 Å along the three-mode iterative path from the open form) shows high strain energy in 6. This helix would likely become partially unfolded during conformational change and figure prominently in the transition-state ensemble for allosteric change. Evaluating whether a mere sliding of the helix without losing contacts is energetically competitive with unfolding would require a much more detailed atomistic calculation.

Free Energy Barrier. A full understanding of the allosteric transition requires consideration of the product energy surface to describe the final approach to the closed structure. This is obtained by using the same strategy describe above but starting for the closed structure (16) deforming toward the open structure (15). The crossing of these two surfaces locates approximately the transition state for the motion (Fig. 3). Obviously small-scale rearrangements are needed to complete the motion from one surface to the other. To locate the transition region we must know the relative stability of the closed and open forms. The conformational change in kinase doubtless is affected by the binding of ligand. The x-ray crystal structure of the closed form is a complex with the inhibitor Ap₅A, which is a bisubstrate analog inhibitor that connects ATP and AMP by a fifth phosphate, thus mimicking both substrates. Before a ligand binds, the open form is more stable, while upon ligand binding, the closed form becomes energetically favorable.

Sanders et al. (36) have experimentally determined the standard binding free energy, DG°, of ATP and AMP, obtaining –6.3 and –4.6 kcal/mol, respectively. If we imagine the ligands are in rapid binding equilibrium, from these data, depending on the concentration of substrate, the effective binding free energy of Ap₅A can be varied in the range of a few kcal/mol. Thus, we use the reaction free energy for the allosteric transition accompanied by binding, DG_eq, of –3 kcal/mol and 0 kcal/mol as reasonable gaps when drawing Fig. 3. As discussed below, however, this topic needs further investigation. The threshold of the ligand bound, the closed state, , is related to the threshold of the open state, , in such a way that the open and the closed states have the same energy when completely unfolded, i.e., , where N_res is the number of residues.

To find where the binding event would occur, the energy profile from the closed form is superimposed on the energy profile from the open form (Fig. 3). The two profiles are actually defined with different coordinates, i.e., rmsd to open and closed forms, respectively. Thus, a mapping between the two coordinates is used as described in Methods.