Article Information

• PDF (143 kb)
• Full Text with Large Figures
• Cited by in Scopus (8)

PubMed

• Articles by Jin Wang
• Articles by Erkang Wang

Related Articles

• Bistability Explains Threshold Phenomena in Protein Aggregat...

  Bistability Explains Threshold Phenomena in Protein Aggregation both In Vitro and In Vivo Biophysical Journal, Volume 90, Issue 3, 1 February 2006, Pages 886-895 Theodore R. Rieger, Richard I. Morimoto and Vassily Hatzimanikatis Abstract Neurodegenerative disease can originate from the misfolding and aggregation of proteins, such as Amyloid-, SOD1, or Huntingtin. Fortunately, all cells possess protein quality control machinery that sequesters misfolded proteins, either refolding or degrading them, before they can self-associate into proteotoxic oligomers and aggregates. This activity is largely performed by the stress response chaperones (i.e., Hsp70). However, the expression level of molecular chaperones varies widely among cell types. To understand the potential consequence of this variation, we studied the process of protein aggregation in the presence of molecular chaperones using mathematical modeling. We demonstrate that protein aggregation, in the presence of molecular chaperones, is a bistable process. Bistability in protein aggregation offers an explanation for threshold transitions to high aggregate concentration, which are observed both in vitro and in vivo. Additionally, we show that slight variations in chaperone concentration, due to natural fluctuations, have important consequences in a bistable system for the onset of protein aggregation. Therefore, our results offer a possible theoretical explanation for neuronal vulnerability observed in vivo and the onset of neurodegenerative phenotypes in neurons lacking an effective heat-shock response. Abstract | Full Text | PDF (209 kb)

• The Membrane Potential and its Representation by a Constant ...

  The Membrane Potential and its Representation by a Constant Electric Field in Computer Simulations Biophysical Journal, Volume 95, Issue 9, 1 November 2008, Pages 4205-4216 Benoît Roux Abstract A theoretical framework is elaborated to account for the effect of a transmembrane potential in computer simulations. It is shown that a simulation with a constant external electric field applied in the direction normal to the membrane is equivalent to the influence of surrounding infinite baths maintained to a voltage difference via ion-exchanging electrodes connected to an electromotive force. It is also shown that the linearly-weighted displacement charge within the simulation system tracks the net flow of charge through the external circuit comprising the electromotive force and the electrodes. Using a statistical mechanical reduction of the degrees of freedom of the external system, three distinct theoretical routes are formulated and examined for the purpose of characterizing the free energy of a protein embedded in a membrane that is submitted to a voltage difference. The W-route is constructed from the variations in the voltage-dependent potential of mean force along a reaction path connecting two conformations of the protein. The Q-route is based on the average displacement charge as a function of the conformation of the protein. Finally, the G-route considers the relative charging free energy of specific residues, with and without applied membrane potentials. The theoretical formulation is illustrated with a simple model of an ion crossing a vacuum slab surrounded by two aqueous bulk phases and with a fragment of the voltage-sensor of the KvAP potassium channel. Abstract | Full Text | PDF (758 kb)

• Mathematical Modeling of the Eukaryotic Heat-Shock Response:...

  Mathematical Modeling of the Eukaryotic Heat-Shock Response: Dynamics of the hsp70 Promoter Biophysical Journal, Volume 88, Issue 3, 1 March 2005, Pages 1646-1658 Theodore R. Rieger, Richard I. Morimoto and Vassily Hatzimanikatis Abstract The heat-shock response in humans and other eukaryotes is a highly conserved genetic network that coordinates the cellular response to protein damage and is essential for adaptation and survival of the stressed cell. It involves an immediate and transient activation of heat-shock transcription factor-1 (HSF1) which results in the elevated expression of genes encoding proteins important for protein homeostasis including molecular chaperones and components of the protein degradative machinery. We have developed a mathematical model of the critical steps in the regulation of HSF1 activity to understand how chronic exposure to a stress signal is converted into specific molecular events for activation and feedback regulated attenuation of HSF1. The model is utilized to identify the most sensitive steps in HSF1 activation and to evaluate how these steps affect the expression of molecular chaperones. This analysis allows the formulation of hypotheses about the differences between the heat-shock responses in yeast and humans and generates a model with predictive abilities relevant to diseases associated with the accumulation of damaged and aggregated proteins including cancer and neurodegenerative diseases. Abstract | Full Text | PDF (222 kb)

• …more

Copyright © 2005 The Biophysical Society. All rights reserved.
Biophysical Journal, Volume 89, Issue 3, 1612-1620, 1 September 2005

doi:10.1529/biophysj.104.055186

Biophysical Theory and Modeling

Quantifying Kinetic Paths of Protein Folding

Jin Wang^*, †, ,  , Kun Zhang^*, Hongyang Lu^* and Erkang Wang^*

^* State Key Laboratory of Electroanalytical Chemistry, Changchun Institute of Applied Chemistry, Chinese Academy of Sciences, Changchun, Jilin 130021, People’s Republic of China
^† Department of Chemistry and Physics, State University of New York, Stony Brook, New York 11794

Address reprint requests to Jin Wang.

Protein folding is one of the most important issues in modern molecular biology. Studying the dynamics is essential in understanding how the protein folds. The question is how the many conformational degrees of freedom can converge to the native state in a finite biological time scale (millisecond to second) instead of cosmological timescale 1. The energy landscape theory of protein folding resolves this issue naturally by assuming there is a bias or funnel toward the native state 2,3,4. This bias is believed to be from the natural evolution. Superimposed on the funneled landscape are the local traps. The slope of the funnel must be steep enough to overcome the traps to reach the folded state. The energy landscape theory is successful in explaining many experiments 5 at both the qualitative and quantitative levels.

According to the energy landscape theory, in general at the initial stage of folding, there are multiple paths toward native state. The discrete paths emerge only when the landscape becomes rough and local traps are important at late stage of folding. Searching for kinetic paths has been a central issue for the folding experimental community for many years 6,7,8,9,10,11,12,13,14,15,16,17,18. Unfortunately, most of the current kinetic folding studies are formulated in terms of the rate dynamics giving only the end results, rather than the paths that represent the full intermediate histories connecting the initial and final ends. It is, therefore, important and natural to formulate the theory in terms of path language. Such a formulation would help to resolve the challenging kinetic path issue of the folding problem and provide a direct tool and language for the theoretical and experimental community to understand each other better. Another advantage of using paths is that the direct integration over paths is normally easier computationally than solving differential equations locally in microscopic details.

Path integral methods since first appearing 19,20 have been successfully applied to many areas in physics 20,21,22 and chemistry 23,24,25. There are so far very limited studies on folding paths. Wang et al. 4 have studied a downhill folding process (very steep funnel) without activation barrier. It is shown that there exists a multiple path to discrete path transition at a temperature higher than the thermodynamic glassy trapping temperature. The relevance to single molecule dynamics is studied 26,27. Olender and Elber and Elber et al. 28,29 studied peptide folding with atomic level simulations and identify some key paths. The purpose of this study is to formulate a diffusive path integral framework for the general case where there exists activation free-energy barriers on the folding landscape, and to identify and quantify the dominant path contributions to the kinetics.

For mathematical simplicity, we study the protein folding problem not at atomic level but at the coarse-grained level—the residue-residue level. This will reduce significantly the computational unforeseeable tasks without the loss of too many important universal features and serve as a guiding force for the more detailed atomic level investigations.

Let us turn to a model Hamiltonian that describes protein folding. To first order approximation, we assume that the energetics that favors bringing two or multiple residues close together from the protein is due mainly to the short-range (in space) hydrophobic driving force. The form of the interactions is −ϵ_ijk…p(α_i, α_j, α_k, …α_p, r_i, r_j, r_k, …, r_p), where ϵ_ijk…p is the multibody coupling strength, r_i is the position of the ith residue, and α_i represents the physical properties of the residue i, for example, hydrophobic charges, etc. Here, we also assume that the environmental solvent effects are already averaged out, resulting in the multibody cooperative hydrophobic interactions among residues upon folding.

We may write down the Hamiltonian energy function of a polypeptide sequence as:

(1)

Suppose there exists a native configurational state n of energy E_n. We can find the probability that configuration a has energy E_a, given that a has an overlap Q with n, where Q is the fraction of native contacts of state a: and N is the total number of native contacts. Q can be used as an order parameter or a reaction coordinate for the physical folding process that measures how close the states are toward native state. Note that for Q=1, the state is in the native folded state and for Q=0, the configurations are in totally nonnative unfolded states.

The conditional probability is obtained directly by averaging over the Gaussian distribution of contact energy ϵ_ijk…p (). By approximating the cooperative multibody interactions σ_ijk…p in the Hamiltonian into the factorization of pair interaction terms σ_ijσ_jk… through a suitable decomposition law such as in the superposition approximation in the theory of fluid, the expression can be simplified as: where m is the order of the interactions (m=2 for two-body interactions, m=3 for three-body interactions, and m=p for p-body interactions), is the average mean energy, and Δϵ is the effective width of the energy distribution per contact.

The configurational entropy S_tot as a function of similarity Q with a given state is treated in details by the previous studies 32,33.

Given the S_tot(Q) and conditional probability distribution obtained earlier, the average numbers of states of energy E and overlap Q with native state n is: This is effectively the microcanonical ensemble description of the thermodynamics. At each stratum of the order parameter or reaction coordinate Q, the set of states is modeled by a random energy model. By the thermodynamic relation of we can obtain the energy and entropy of the biomolecular folding as: and where s_tot(Q)=S_tot(Q)/N. The entropy vanishes at a characteristic temperature: which signals the trapping of the polypeptide chain into a low-energy conformational state within the stratum characterized by Q. Notice that when Q=0 (nonnative unfolded states),

From the thermodynamic expression of the energy and the entropy given above, we can easily obtain the expression for the free energy per contact as 33: where The free energy is composed of three terms, the entropy, the native driving force, and roughness contribution of the energy landscape. In the parameter space in (δϵ_n, Δϵ, T), the expression above can have a double minimum structure in the reaction coordinate Q with one minimum at low Q corresponding to the nonnative states separated by a barrier from another minimum at high Q corresponding to the native folded state. As the cooperativity measured by multibody interaction order m increases, the free-energy minimum of nonnative states and native folded state shift toward Q∼0 and Q∼1, respectively. To the extent that this approximation is good (m →∞), we can equate the free energies of the nonnative states and native folding state to obtain the folding transition temperature (F(Q=0)=F(Q=1)):

Take the ratio of folding temperature and trapping temperature, we obtain:

(2)

Therefore, maximizing the ratio of the energy gap (or the slope) versus the roughness of the underlined energy landscape becomes the criterion for the thermodynamic stability of folding, implying a funneled energy landscape.

Under the free-energy profiles, the equation of motion for native contact Q formation can be formulated as:

(3)

When taking into account the combination of multibody interactions (up to the six-body interactions because the order of the hydrophobic multibody interactions beyond two-body interactions is typically ranging from three or four up to six), the free energy becomes:

(4)

where α, c₁, c₂, c₃, 1−α−c₁−c₂−c₃ are the coefficient mimicking the relative importance of the order of multibody 2,3,4,5,6 interactions. N is the length of the polypeptide chain. S_tot(Q) ≈ S₀(1−Q)−QLog(Q)−(1−Q)Log(1−Q), where S₀ ≈ ln(10/2.718); 10 in the ln is the degrees of freedom per residue whereas factor 2.718 in the ln takes into account the constraints of the phase space upon collapse. The first term of the entropy is the entropy loss forming a contact whereas the rest of the two terms is responsible for the entropy associated with the possible ways of forming a contact.

We can now formulate the dynamics with the Onsager-Machlup 21 functional path integral as:

(5)

Each path in the path integral contributes a weight, but not every path gives the same contribution. In fact, the contribution from the paths to the weight is on the exponential, so the dominant paths with the largest weight contribute significantly larger than the ones with the subdominant or even smaller weights. We can then approximate the path integrals with a set of dominant paths and ignore the subleading terms. One can easily see to find the paths with the optimal weights, the dominant paths should satisfy the Euler-Lagrangian equation (see Fig. 1):

(6)

(7)

(8)

Display large version of this figure

Figure 1

Folding paths from initial to final configuration.

The equation of motion can be integrated out to obtain:

(9)

The free energy as a function of Q at various temperatures T is plotted in Figure 2AC, (for mixed but mainly four-body, five-body, and six-body interactions). The potential V as a function of Q is plotted in Figure 3AC (for mixed but mainly four-body, mixed but mainly five-body, and mixed but mainly six-body interactions) as a function of Q at folding temperature T_f.

Because around folding temperature, the free energy F as a function of Q often has a double-well shape (Fig. 2, with given parameters specified later in the article) with one well corresponding to the nonnative unfolded states and the other one corresponding to the native folded state. The free-energy barrier is closely linked with the cooperative nature of multibody hydrophobic interactions for protein folding. We have done a careful analysis with different degrees of cooperativity in the inherent interactions in the Hamiltonian with mixed but mainly four-body interactions, mixed but mainly five-body interactions, and mixed but mainly six-body interactions. We see from Figure 2AC, that as the degree of cooperativity increases, the barrier height increases too. In other words, for low cooperativity the barrier is small, but for high cooperativity the barrier is large. One can substitute the shape of F(Q) into the expression of V and obtain the shape of the potential V as a function of Q (Fig. 3). Again, we see that the V has a minimum. The position of the minimum is close to the original minimum in Q in the free-energy profile F. The dominant contribution for the paths are from solving the equation of motion for Q. The effective potential U=−V. In the long time limit, there exists possibilities that the paths go back and forth many times from the hill (maximum) in the effective potential U (the minimum or the valley in V since U=−V) corresponding to the nonnative states to the other bounce-back point near the native state, where the value of U at the bounce is equal to that of the hill. This corresponds to the traversal of multiple times passing through the barrier to reach the native folded state. These oscillating back and forth solutions are called the instanton solutions 35,36,37,38,39. In Figure 4AC, the instanton solutions are shown for dominant four-, five-, and six-body interactions. Each instanton (antiinstanton) corresponds to one transition from nonnative (native) to native (nonnative) states. The dominant path is composed of multiple instantons. The contribution can be summed in the dilute gas approximation by assuming no instanton-instanton interactions to obtain the final contribution to the probability of folding. The one instanton contribution to the weight is given by:

(10)

Display large version of this figure

Figure 2

Free energy F as a function of Q at various temperatures near folding temperature. T_f for dominant four-body (A), dominant five-body (B), and dominant six-body interactions (C).

Display large version of this figure

Figure 3

Potential V as a function of Q at folding temperature. T_f for dominant four-body (A), dominant five-body (B), and dominant six-body interactions (C) when D=D(Q=0).

Display large version of this figure

Figure 4

Instanton path Q as a function of time t of the protein folding dynamics at folding temperature. T_f for dominant four-body (A), dominant five-body (B), and dominant six-body interactions (C) when D=D(Q=0).

The probability is determined from the optimal paths contributed by the sum of the multiinstanton contribution. It can be written as (Figure 3 and Figure 4):

(11)

The above expression can be easily evaluated in the Laplace representation s:

(12)

By inverting the Laplace transform, we obtain:

(13)

(14)

When V(Q_min)=V(Q_max) as is the case in instantons, the expression is simplified as:

(15)

We take number of residues as N=30, roughness or spread of the landscape Δϵ=1, the bias or slope of the landscape toward folded state δϵ=3 and α=0.05, c₁=0.05, c₂=0, c₃=0, and 1−α−c₁=0.90 (for mixed two-, three-, and four-body interactions, but dominant four-body interactions); α=0.05, c₁=0.05, c₂=0.05, c₃=0, and 1−α−c₁−c₂=0.85 (for mixed two-, three-, four-, and five-body interactions, but dominant five-body interactions); α=0.05, c₁=0.05, c₂=0.05, c₃=0.05, and 1−α−c₁−c₂−c₃=0.80 (for mixed two-, three-, four-, five-, and six-body interactions, but dominant six-body interactions). The diffusion coefficients are given as 40: D(Q)=D₀exp[−S₀(Q)] for T<T_g; and D(Q)=D₀exp[−β²ΔE(Q)²] for 2T_g<T; and D(Q)=D₀exp[−S₀(Q)+(β_g(Q)−β)²]ΔE(Q)²] for T_g<T<2T_g. Here,

Figure 4AC, shows the multiinstanton solutions at T=T_f for dominant four-, five-, and six-body interactions, respectively. Figure 5AB, show the temperature dependence of the logarithm (ln) of the folding rate (K)−lnK for dominant four-, five-, and six-body interactions when diffusion coefficient is constant with D=D(Q=0) and when diffusion coefficient is reaction coordinate dependent with D=D(Q).

Display large version of this figure

Figure 5

The logarithm of the kinetic protein folding rate ln K as a function of temperature T (in the units of folding temperature T_f) for dominant four-body, dominant five-body, and dominant six-body interactions when diffusion coefficient is constant, D=D(Q=0) (A) and when diffusion coefficient is Q dependent, D=D(Q) (B).

As we can see folding kinetic rate has a “bell”-like shape dependence with respect to the temperature. At high temperatures, the folding kinetic rate is small. This is due to the instability of proteins at high temperature. On the other hand, at low temperatures, the folding rate drops again. This is due to the possible trapping into the local valleys. Thus, the temperature varying kinetics provides a way of exploring the structure of the underlined folding energy landscape. The maximal rate for folding happens at certain optimal temperature. This is in good agreement with kinetic folding experiments (Chevron rollover) and theory/simulation studies 41,42,43,44,45,46. We also observed that as the cooperativity of the inherent hydrophobic interactions increases (from dominant four-body interactions to dominant six-body interactions), the free-energy barrier increases (as shown in Fig. 2), and the associated kinetic rate decreases. Furthermore, we can see that when the diffusion coefficient depends on the reaction coordinate, the kinetic rate for folding is significantly changed, especially at the low temperature regimes. In the low temperature regimes, the thermodynamic free-energy barrier for folding is less and less compared with the corresponding higher temperature case (see Fig. 2), and the effect of the diffusion on kinetics becomes more and more important.

This indicates that the kinetics is not only controlled by the inherent thermodynamic free energy but also by the diffusion. This is particularly important because the fast folding experiments are now approaching the speed limit where the kinetics of pure diffusion can be measured 47,15,48.

We can simplify the expression of the kinetic rate by assuming that diffusion coefficient is relatively small. In this case, we can substitute the instanton solution to the action of the probability expression of the path integral 49 and obtain analytic form of equilibrium probability as:

(16)

(17)

The effective activation energy for transitions from nonnative unfolded state at Q=Q_min to the transition state Q^# is given by:

(18)

It is very important to realize that the current formalism implies both the diffusion and thermodynamic free-energy barrier control the kinetics of protein folding as mentioned above. When the underlying process is barrier limited, both the thermodynamic barrier and diffusion contribute to the kinetics although free energy contribution might be larger. The role of diffusion is to modify the effective free-energy profile and the corresponding barrier. In the case where there is no inherent free energy barrier, the kinetics is controlled by diffusion. Thus, the formalism in this article provides a route to look for the switching roles from thermodynamic-barrier-driven kinetics to downhill diffusion-driven kinetics, which is quite relevant for the experimental study of fast folding proteins where the speed of folding is determined from the thermodynamic-driven to the essentially diffusion-controlled process 47,15,48.

The current formalism can also be used to discuss the transition state property of protein folding. In the case of constant diffusion, the current formalism reduces to the normal transition state theory and kinetics is controlled by the free energy barrier. As mentioned above, when the diffusion coefficient is not a constant, the kinetics is controlled by both the free energy barrier and diffusion. In the case when the thermodynamic barrier is large, the kinetics is dominated by the free energy profile. On the other hand, when the thermodynamic barrier is moderate, the effect of diffusion will come into play by modifying the original free energy profile. Both the position and value of the resulting effective transition state free energy will be shifted. So the kinetics will be modified by this shift of the original transition state. Details of the study will be given in a future publication.

Let us discuss the possible connections of our approach with another set of experimental observations 50. When the diffusion coefficient is a constant, the kinetics is controlled by the free energy profile as we have derived above. In the barrier limited case, if the barrier is caused mainly by topology instead of heterogeneity of the interactions, then the free energy barrier is mainly from entropy contribution of loop contacts 32,33,51. Thus, the free energy change with respect to the mean sequence length of making contacts can be shown as:

(19)

Thus, the free energy barrier is linked to the average sequence length of making the contacts The effect of increasing the mean loop length is to increase the barrier height. So the kinetics is faster (slower) when the mean contact distance is small (large). When the diffusion coefficient is not constant and interaction heterogeneity are taken into account, the free energy dependence on the mean contact distance might not be as strong. Further detailed investigations on this are needed and will be carried out in a future publication.

We discussed in this article the long-time dynamics of folding. In principle, the short-time dynamics can be revealed by solving the Euler-Lagrangian equation for the optimal paths. Because the time is short, the solutions typically don’t have enough time forming multiple instantons. Finding dominant paths becomes solving ordinary differential equation for fixing two end points. One can expand around the dominant solution up to quadratic order and obtain the contribution to the probability of folding. In general the results are good for short times and the kinetics is usually nonexponential. This is in contrast with the long-time case where the dynamics is usually controlled by the longest timescale as we discussed here.

We obtain in this study the optimal instanton paths that determines the folding rate dynamics in the long time limit. We should mention that the optimal paths are actually a set of paths in the multidimensional configurational space. They represent the dominant flow of paths directed toward the native state. At low temperatures, the folding might be trapped into the local valleys, while the current continuous path approach can give some qualitative features as to approximately when the continuous flow of paths might break down; instead, the more appropriate approach seems to be the discrete version of the path integral we presented here. The formulation is currently under development. With this formulation, one can study and understand the transition from the multiple paths to the discrete path transition in the case of activated folding transition.

The kinetic rate dynamics is often studied by the Fokker-Planck type rate equations (or Brownian dynamics). This approach to the kinetics is mathematical, related to the path integral formulations presented here but emphasizing different aspects. Although the path method concentrates on intermediate processes and the corresponding contributions to the final kinetics, the Fokker-Planck type rate equation approach concentrates more on the end results. Therefore, it is convenient and advantageous to address the kinetic path issues for protein folding in the path integral formulation.

It is worth mentioning that biomolecular recognition (binding) often involves large fluctuations and conformational changes 52,53,54,55,56,57,58,59; sometimes local unfolding 60,61 for induced fit 62 is necessary, so in general folding and binding are dynamically coupled. It is tempting to study the kinetics of the folding-binding process using the current developed path integral methodology (J. Wang, K. Zhang, H. Y. Lu, and E. K. Wang, unpublished data). The crucial question would be what are the dominant kinetic paths for the folding-binding process in nature.