Article pubs.acs.org/JPCB
Relation between Protein Intrinsic Normal Mode Weights and PreExisting Conformer Populations Beytullah Ozgur,†,∥ E. Sila Ozdemir,‡,†,∥ Attila Gursoy,*,§,† and Ozlem Keskin*,‡,† †
Center for Computational Biology and Bioinformatics, ‡Chemical and Biological Engineering, and §Computer Engineering, College of Engineering, Koc University, 34450 Istanbul, Turkey S Supporting Information *
ABSTRACT: Intrinsic fluctuations of a protein enable it to sample a large repertoire of conformers including the open and closed forms. These distinct forms of the protein called conformational substates preexist together in equilibrium as an ensemble independent from its ligands. The role of ligand might be simply to alter the equilibrium toward the most appropriate form for binding. Normal mode analysis is proved to be useful in identifying the directions of conformational changes between substates. In this study, we demonstrate that the ratios of normalized weights of a few normal modes driving the protein between its substates can give insights about the ratios of kinetic conversion rates of the substates, although a direct relation between the eigenvalues and kinetic conversion rates or populations of each substate could not be observed. The correlation between the normalized mode weight ratios and the kinetic rate ratios is around 83% on a set of 11 non-enzyme proteins and around 59% on a set of 17 enzymes. The results are suggestive that mode motions carry intrinsic relations with thermodynamics and kinetics of the proteins.
■
INTRODUCTION
troscopy, and single-molecule FRET are among the experimental methods to study motions.3,8 It has been shown that a protein does not possess a single, static structure but rather samples a repertoire of substates at equilibrium around its native state. The so-called “pre-existing equilibrium model”9−13 aims to explain the conformational rearrangements at equilibrium. The model is based on the energy funnel landscape theory14,15 in which a protein exhibits different conformations or substates separated by energy barriers obeying a statistical thermodynamics distribution. The heights of the barriers reflect the time scales of the conformational exchanges. The thermal fluctuations provide innate flexibility for the protein, and these motions enable the protein to span an ensemble of substates directing protein function and behavior such as ligand binding, function inhibition, or catalysis. The substate with the binding pocket most complementary to the ligand binding will be selected for binding. Ligand binding simply shifts the distribution of the populations of each conformer. The conformer selected initially may not have the lowest energy; however, in the bound form it will be the most stable which will result in a “population shift” toward this conformer.16−18 Pre-existing equilibrium model is supported by several experimental studies including nuclear magnetic resonance (NMR) relaxation and single-molecule
Intrinsic dynamics of a protein are related to its biological function. The time-dependent, dynamic nature of proteins highlights the importance of protein motions. Proteins tend to undergo conformational changes in order to properly function. For example, ligands can induce conformational changes to modulate binding affinities, or chemical modifications (such as phosphorylation) can stabilize/destabilize the conformations to control different cascades in pathways.1,2 Proteins exhibit a broad spectrum of motions ranging between fast scale motions in picoseconds to nanoseconds and slow scale motions in microseconds to milliseconds. The former motions are coupled to local, small amplitude fluctuations between close energy states that have small energy barriers around kT.3 Side chain rotations and small amplitude backbone fluctuations are examples of fast scale motions. NMR relaxation experiments4−6 allow investigation of fast scale dynamics in terms of bond fluctuations. Infrared or fluorescence correlation spectroscopy can also provide insights about local fluctuations.3 On the other hand, slower scale motions are collective, including largeamplitude motions of secondary structure elements, subunits, or domains fluctuating between energetically distinct substates separated by energy barriers of several kT.3,7 These slow-scale fluctuations enable a protein to undergo relatively more stable states in which interconversions between these states carry importance about the protein function including enzyme catalysis, signal transduction, and biomolecular interactions. Hydrogen−deuterium exchange, X-ray diffraction, NMR spec© XXXX American Chemical Society
Special Issue: Klaus Schulten Memorial Issue Received: October 14, 2016 Revised: February 9, 2017 Published: February 9, 2017 A
DOI: 10.1021/acs.jpcb.6b10401 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B experiments.19−25 These techniques provide strong evidence about the conformational ensemble dynamics of several proteins including RNase,16,26,27 adenylate kinase (ADK),17,28 calmodulin,29−31 and cyclophilin A.32 Interestingly, an IgE antibody Spe7 was found to coexist as two distinct substates (conformers) in its native state where each conformer was bound to a different ligand. The structures of these two conformers were also revealed by crystallographic methods.33−35 Such studies help to establish an ensemble based approach for elucidating protein function. Characteristic motions of the cyclophilin A (CypA) cis−trans isomerase were studied using NMR spectroscopy,32 and this study elucidated a relation between enzyme flexibility and catalysis. The conformational changes taking place during the turnover were found to pre-exist in the free enzyme with frequencies similar to catalytic turnover rates, implying the catalysis to be an inner property of the free protein.32 In the same study authors concluded that CypA dynamics associated with catalysis is a built-in property of the enzyme that is also apparent in the free CypA. For the cytochrome P450, the ensemble dynamics and interconversion rates were determined by crystallographic techniques36,37 and also the study36 presented the interconversion rates between corresponding substates. NMR relaxation and dispersion experiments also determine the relative populations of the ensemble substates and chemical shift differences between corresponding substates.5,38 Characterization of catalytic motions of free RNase16 revealed three substates pre-existing together in conformational ensemble including free enzyme, enzyme−substrate complex, and enzyme−product complex. Such studies established these innate slow scale motions to be relevant to functional motions. The ligand binding acts to stabilize the pre-existing conformer, altering the initial and final probability distributions of substates. All these experimental methods help to characterize slow motions in terms of ensemble structures, kinetic and thermodynamic properties. Still, low stability of the ensemble intermediates prevents exact determination of dynamics occurring in a protein. Computational approaches can complement understanding of the lowly populated substates. Computational approaches such as molecular dynamics (MD)39 simulations reach fast scale motions and midscale motions using long simulation times; even slow motions have become accessible with recent advances.40−42 In an MD study, Wang and Lu studied conformational dynamics of adenylate kinase using a coarse-grained two-well model in order to overcome the computational bottleneck.43 In that study, the authors determined the opening and closing pathways of ADK enzyme and predicted the kinetic rates consistent with experimental findings. Further studies have been performed to configure more accurate and comprehensive conformational change models of ADK applying new MD techniques.44,45 MD coupled with adaptive biased force and thermodynamic integration free energy methods were used to reproduce the partially closed apo state of maltose binding protein (MBP) and suggested that sugar binds more tightly to the partially closed than to the open apo state.46−48 Also, kinetics of conformational change between apo and closed states of cyclophilin A and NtrC proteins have been studied using MD simulations.49−51 Several other computational methods have emerged to predict collective, slow-scale motions of conformational dynamics including coarse-grained normal-mode analysis.
Elastic network models (ENMs) depend on the basic assumption of normal-mode analysis where the energetic fluctuations of the protein have been considered to follow a harmonic pattern. ENMs are based on the intercontact topology of a protein and have the ability to define the important motions.52 Gaussian network model (GNM)53,54 and anisotropic network model (ANM)55 are the two forms of the ENMs where the former one supplies normal mode fluctuations in an isotropic manner containing no information about the directionality, whereas the latter one provides fluctuations in an anisotropic manner supplying direction information on normal modes. These two methods have been employed for a broad range of studies including predictions of Debye−Waller factors 53−56 and conformational transitions.57−64 The results obtained in these studies provide important insights about the conformational dynamics. Some of these normal modes, also called “soft modes”65 of several proteins, were found to be strongly correlated with the large amplitude conformational changes of these proteins observed upon ligand or protein binding.60,61,66−71 Soft modes can be defined as pre-existing valleys on the conformational energy landscape between different substates.72 These modes are necessary to describe the structural changes between two distinct forms of the same protein or/and protein complexes (i.e., open and closed forms).2,73 Besides, mode frequencies provide insights about the size of conformational change.71 Also, small size and short-lived atomic fluctuations which are experimentally difficult to detect can be reproduced and hypothetical transition pathways between these conformational states can be detected by ENM studies.63,74−76 Krebs et al. showed that approximately 3800 motion pairs could be described by evaluating only two of the lowest energy normal modes.77 The results obtained in these studies provide important insights about the conformational dynamics. ENMs are also helpful to explain conformational changes in proteins upon ligand binding.78−80 These changes can be predicted by slow modes defining global motions as well as some fast modes coupled with them. The energy transfer can occur between coupled modes. The motion of one of the modes drives the motions of the other mode, resulting in conformational changes in a protein.81 Bahar and colleagues showed that by use of ENMs of unbound (apo) proteins, conformational changes experimentally observed upon ligand binding can be reconfigured.61,65 In their related studies, it has been concluded that the motions exhibiting the largest conformational change calculated by ENMs for the monomer of unbound proteins correlate with the experimentally observed structural changes upon binding.61,65,82,83 Also Wako and colleagues reported that a full-atomic system to predict conformational changes of a protein from its apo to its holo form can be represented by a linear combination of the displacement vectors of atoms calculated from apo form using ENMs.84 Recent studies also have shown that ANM and GNM methods can be useful for docking tools by predicting conformational changes upon ligand or protein binding.78,85,86 Coarse grained methods are also frequently combined with molecular dynamics to give insights about the mechanisms of large scale motions.87,88 Gur et al. reported a close similarity between collective motions of proteins exhibited by MD simulations and those predicted by ANM.89 In this computational study, they concluded that the lowest frequency ANM modes facilitate transition between the most probable substates and successfully capture the global dynamics of proteins. B
DOI: 10.1021/acs.jpcb.6b10401 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B
Figure 1. Schematic representation of the proposed model. Conformer 1 (S1) and conformer 2 (S2) are in equilibrium with forward and backward rate constants of k1 and k−1. Normal mode analysis provides modes with respective normalized weights.
equilibrium models. We assume that different conformations (substates) of a protein can be analyzed by vibrational modes. Figure 1 schematically illustrates our approach. The two conformers (substates S1 and S2) of a protein are at equilibrium with forward and reverse rate constants as k1 and k−1 (S1 ⇌ S2) with an equilibrium constant of K (a dimensionless constant assuming first order reaction for both forward and reverse reactions). If the transition free energies of the two substates are ΔE1 and ΔE−1, then
In this study, we focus on modal analysis of pre-existing conformers and investigate the relation between modes and thermodynamics/kinetic properties of ensembles. We collected a set of proteins that have both an open and a closed conformation. We further gathered kinetic/thermodynamic data of these proteins from experimental studies. We performed parameter free ANM (pfANM52) analysis for each of two preexisting conformations of the same protein. We propose that frequencies (square root of eigenvalues) of the modes provide information about the kinetic rates of conformations. The frequency of the normal mode that drives one conformation to another is related to the interconversion rate between the two conformations (substates) and thermodynamic populations of two substates. We show that the ratios of interconversion rate constants of conformers are correlated with the ratios of the normalized weights of the modes.
■
k1 = A1 exp{−ΔE1/(RT )} k −1 = A −1 exp{−ΔE−1/(RT )}
(1)
where A1 and A−1 are pre-exponential factors and R and T are the gas constant and absolute temperature. The energy barrier between the two states (ΔΔE = ΔE1 − ΔE−1) will be proportional to the relative probability of each state. If the normal modes of each conformer are found, with respective normalized weights (eq 8), then we hypothesize that these pre-existing substates can be studied by vibrational mode analysis and there exist some modes that will drive conformer 1 to conformer 2 (say ith mode of conformer 1, S1, and let the associated normalized weight of this mode be wi) and similarly conformer 2 (S2) to conformer 1 (S1) (say jth mode of conformer 2, S2, and the associated normalized weight of this mode is wj); then there exists a relation between the ratios of the normalized weights of the modes and the ratios of the rate constants of these conformers (α is the correlation coefficient of this relation, K is the equilibrium constant) as wj k −1 1 = =α k1 K wi (2)
THEORETICAL MODEL AND PROTEINS STUDIED
Advances in biophysical techniques allow investigation of protein dynamics in detail within the time scales of functional motions. These studies elucidate time-dependent behavior of proteins enabling the study of the so-called energy landscape of proteins together with kinetic and thermodynamic properties of the landscape. Related ensemble substates from native state structures of proteins have been determined using NMR or Xray crystallography. Such studies accumulate in literature, supplying findings about the thermodynamics (relative populations of substates) and kinetics (time scales of interconversions between different substates) of the protein in detail. As aforementioned, normal-mode analysis can explain the behavior of proteins undergoing conformational changes between several substates of the same protein albeit population inference.1,59,71,90 Kinetic and transition state theories state that chemical reactions pass through an unstable structure called the transition state, which lies between the chemical structures of the reactants and products. The theory indicates that kinetic rate constants are proportional to the frequency of the vibrational mode responsible for converting a transition state to the product. Here, we make an analogy for pre-existing
In principle, each mode contributes toward the overall motion of a protein proportional to its weight; thus the weighted average of all modes will provide the overall motion of the protein. In pre-existing equilibrium model, a protein at its apo-form can sample all its accessible conformations. Since some modes may be responsible for driving the protein from its apo-form to its accessible conformations, the normalized weight C
DOI: 10.1021/acs.jpcb.6b10401 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B Table 1. Proteins Used in This Studya protein name
no. of residues
no. of domains
PDB I (open)
PDB II (closed)
RMSD (Å)
KL
lit. data
ref
adenylate kinase Spe7 34E4 Fab antibody MAD2 U2AF TF calmodulin TIM NtrC (at 25 °C) MBP thrombin RNase DHFR cholesterol oxidase CypA glucokinase Ltn PKA LDH AChR YopH HLADH arginine kinase HRas CPR FFDomain CheA Erk2
214 120 224 187 195 148 247 123 366 251 124 154 500 164 424 60 334 330 370 278 374 347 166 603 56 179 351
1 1 2 1 2 1 1 1 1 1 1 1 2 1 2 1 1 1 1 1 2 1 1 1 1 1 1
4ake (A) 1oaq (H) 1y0l (H) 1duj (A) 2yh1 (A) 1cfc (A) 1ypi (A) 1dc7 (A) 2v93 (open) (A) 1sgi (B) 7rsa (A) 5dfr (A) 1coy (A) 2cyh (A) 1v4s (A) 1j9o (A) 4nts (A) 4cuj (A) 4aq9 (A) 1ypt (A) 1ju9 (A) 3m10 (A) 4efl (A) 3es9 (A) 2kzg (A) 1i5d (A) 1erk (A)
1ake (A) 1ocw (H) 3cfj (H) 1s2h (A) 2yh0 (A) 1cfd (A) 7tim (A) 1dc8 (A) 2v93 (closed) (A) 1rd3 (B) 1u1b (A) 1drb (A) 3cox (A) 1oca (A) 1v4t (A) 2jp1 (A) 1atp (E) 4cuk (A) 4aq5 (A) 2i42 (A) 1hld (A) 1m15 (A) 1gua (A) 3qe2 (A) 1uzc (A) 1b3q (A) 2erk (A)
5.31 1.23 1.93 12.06 12.37 4.98 0.46 1.78 4.97 1.25 0.68 0.79 0.24 0.64 9.70 9.40 2.99 0.78 0.52 0.57 0.98 3.11 2.83 4.20 1.37 1.70 0.69
0.37 0.32 0.21 0.56 0.68 0.78 0.07 0.46 0.62 0.12 0.36 0.51 0.24 0.51 0.32 0.63 0.28 0.09 0.45 0.14 0.55 0.52 0.38 0.36 0.67 0.24 0.14
6500/2000 17/58 0.5/1.55 1/2 3/1 1/1 46700/2500 1900/11600 226/197 224/34 1615/85 0.064/0.038 0.13/0.0417 1080/60 7.16/0.8 0.28/0.38 953/66 3020/580 1/0.000003 0.064/0.0053 620/64 652/28 250/135 1/1 1920/22 13/0.25 240/60
17, 91 34 92 93 94 30 95 96 97 98 16 99 100 101 102 103, 104 105 106 107 108 109 110 111 112 113 114 115
a
The second column presents the length of each protein after the nonmatching residues are removed. The third column gives the number of domains of the proteins. The fourth and fifth columns present PDB structures used for constructing pfANM Hessian matrix and determining experimental displacements. The seventh column gives the collectivity of the corresponding motion between the two PDB structures. The selected chain identifiers of proteins are shown in parentheses. PDB structures for adenylate kinase, RNase, TIM, and NtrC include open and closed structures of each protein as PDB I and PDB II structures, respectively, except Fab antibody. The nonbinding form of the Fab antibody is represented by 1y0l, and 3cfj is the ligand bound form of the antibody. Spe7 includes two unbound conformer structures named Ab1 (1oaq) and Ab2 (1ocw). The literature data column shows kinetic rate constants (s−1) of the substates. The last column presents the references of experimental findings about the kinetic conversion rates between substates.
of the mode will be proportional to the kinetic rate constant of the selected conformer. Data Set of Proteins with Two Substates. In this study we collect experimental data for 27 different proteins that have different substate structures together with the interconversion rates between these structures. Most of the proteins have at least two substate structures determined by either NMR or Xray crystallography. One of the states corresponds to the closed state, whereas the other is the open state. Most of the substate structures are without a ligand such that the closed states correspond to the less populated but still coexisting conformation with the open substructure of the same protein. For some proteins (especially enzymes), the structures and kinetic data are measured with ligands, so in these cases, the most populated structures are the closed states. As a result, we have a list of 54 PDB structures (two for each protein) and their respective kinetic or thermodynamic data. The list of proteins is presented in Table 1 (the first column) with the referred substates (fourth and fifth columns), the rootmean-square deviation (RMSD) between the two substates and the collectivity (see Methods). The proteins span a wide range of functional classes: enzymes, antibodies, binding and signaling proteins. The RMSD between the two substates ranges from 0.24 to 12.37 Å. The superimposed regions of the proteins to calculate RMSDs are given in Table S1 in Supporting
Information. Small molecules bound to the proteins are listed in Table S2. Choice of Normal Modes Driving the Proteins. We perform modal analysis of the two substates of each protein in our data set using pfANM (see Methods). The slow motions, as aforementioned, can lead to conformational exchanges between distinct states. Then, the modes that drive one substructure to another are determined. We analyze individual modes in order to determine similarity between experimental displacements and normal modes using two similarity scores which are mode correlations (corr) and overlaps (I) (see Methods). We adopt two strategies to choose which modes to use in our analysis: 1. The modes with the best correlations: We simply choose the mode with the highest correlation (corr) value among the slowest 10 modes for each conformer. The list of modes and their correlations is given in Table S3. The associated normalized weights of the modes, for both open and closed states, are given in the last column of the table. 2. The modes with the best overlaps: We simply choose the mode with the highest overlap (I) value among the slowest 10 modes for each conformer. The list of modes and their overlaps are given in Table S4. D
DOI: 10.1021/acs.jpcb.6b10401 J. Phys. Chem. B XXXX, XXX, XXX−XXX
Article
The Journal of Physical Chemistry B Before providing the results of our analysis, first we give a brief summary on how we collect the experimental data for each protein. Adenylate kinase (ADK) performs reversible catalysis of ATP and AMP into ADP. The protein interconverts between an open and closed conformation directed by the motions of two ligand binding domains, the ATP lid and AMP binding domain (AMPbd). Once the catalysis is accomplished, the enzyme releases its products and adopts the open form.28,116 The crystal structures of ADK from several organisms have been determined both in free and in complex forms.117 NMR studies and comparison of the catalytic activities of the enzyme in two organisms, one mesophilic and the other thermophilic, indicate that reduced activity of the thermophilic enzyme is related to lid-opening rate of the thermophilic enzyme and the enzyme exists at least at two different conformations.17 The study also provided the steady-state catalytic reaction constants and the rate constants for opening/closing motions. Open ADK is reported to be converted to closed ADK with a rate constant of 6500 s−1, whereas a rate constant of 2000 s−1 is observed for the reverse reaction with NMR.91 Monoclonal antibody, Spe7, exists as two free conformers, Ab1 and Ab2 (pre-existing at equilibrium), within solution, and the population distribution of each conformer depends on the type of ligand as determined by X-ray crystallography.33,34 Preexisting equilibrium kinetic studies showed that the relative populations are 22% and 78% for the closed (Ab2) and open (Ab1) forms. The forward and reverse rate constants (k1 and k−1, respectively) are 17 s−1 and 58 s−1 for this interconversion reaction. RMSD between two structures is 1.23 and the collectivity of the motion is calculated as 0.32. Heavy chains are considered here. 34E4 Fab antibodies are catalytic antibodies that catalyze the conversion of benzisoxazoles to salicylonitriles.118 It follows a conformational diversity due to the presence of pre-existing conformers. There are two free forms of the enzyme in equilibrium interconverting between each other, and only one of these two forms is eligible to bind to its substrate. The substrate binding and nonbinding forms of the enzyme coexist together in the solution favoring the nonbinding form in a ∼3:1 ratio, forward and reverse rate constants of 0.5 and 1.55 s−1, respectively.92 Heavy chains are compared here. Mitotic arrest def iciency 2 (MAD2) is a protein involved in the spindle checkpoint that has two natively folded states. Mad2 adopts two topologically and functionally distinct native folds without ligand binding and chemical distributions. In the mixture of two alternative conformations of MAD2 (open and closed forms), closed-Mad2 and open-Mad2 were found in a 2:1 ratio, since concentrations (populations) are related to the rate constants (k1/k−1 = Population2/Polulation 1), the equilibrium rate constant (K [eq 2]) is 2.93 The RMSD between two structures has been calculated as 12.06 with a collectivity of 0.56. U2AF transcription factor has a large, multidomain 65 kDaA subunit (U2AF65) which has a crucial role in the assembly of pre-mRNA splicing complexes. RNA sequence is recognized by the tandem RNA recognition motif (RRM) domains (RRM1− RRM2) of U2AF65.119 The equilibrium rate constant (K) between open and closed conformations of unbound U2AF65 is found to be 3.94 The RMSD between active (closed) and inactive (open) forms of the enzyme is calculated as 12.37, and the collectivity of the experimental displacement is 0.68.
Calmodulin is a messenger protein that is activated by Ca2+ binding in eukaryotes. The structures of the unbound (apo) and bound (Ca2+ loaded) states of calmodulin have revealed that a closed to an open conformation change occurs upon Ca2+ binding. By performing 15N spin relaxation analysis, an open population of 50% and a rate constant of 2.7 × 104 s−1, which is the same for both closed → open and open → closed, have been obtained.30 The RMSD between two structures has been calculated as 4.98 with a collectivity of 0.78. Triosephosphate isomerase (TIM) catalyzes isomerization of a ketone to an aldehyde. Motion of loop 6 in the active site of TIM repositions the key residues for catalysis.120 The steadystate kinetic rates for open to closed and closed to open forms of the enzyme have been found as k1 = 46 700 s−1, k−1 = 2500 s−1, respectively.95 Modal analysis of the open and closed conformers of the protein was performed, and RMSD between two structures has been calculated as 0.46 with a collectivity of 0.07. Nitrogen regulatory protein C (NtrC) activation occurs when a phosphate is transferred from NtrB to D54 in NtrC.42 NMR studies suggest that the equilibrium is far shifted toward the active conformation upon phosphate binding. However, the population of active species in wild-type, unphosphorylated NtrC can be estimated to be between 2% and 10%.121 Kinetic rates from open to closed form and from closed to open form of the enzyme have been found as k1 = 1900 s−1, k−1 = 11 600 s−1, respectively.96 Experimental conditions of chosen protein structures from PDB to analyze population of active and inactive species were defined at 25 °C. The RMSD between two structures has been calculated as 2.52 with a collectivity of 0.46. MBP is a maltose binding protein, and it has an NMR structure with a predominantly open form coexisting in rapid equilibrium with a minor closed species.69 Pre-existing equilibrium kinetic studies showed that the relative populations are 95% and 5% for the open and closed forms.69 Also molecular dynamics (MD) analyses give forward and reverse kinetic rates between open and closed forms as k1 = 226, k−1 = 197.97 The RMSD between two structures has been calculated as 4.97 with a collectivity of 0.62. Thrombin is an important protein in coagulation cascade and pre-exists as slow form (free thrombin) and fast form (Na+ bound) states in solution.98 Stopped-flow fluorescence measurements showed that the free enzyme exists in equilibrium between two forms, E* and E, that interconvert with kinetic rate constants k1 = 224 s−1 (E* to E) and k−1 = 34 s−1 (E to E*) at 25 °C.98 Of these forms, only E can interact with Na+. While the E* form is negligibly populated (
