Modeling protein thermodynamics and fluctuations at the mesoscale

We use an extended Gō model, in unfrustrated and frustrated variants, to study the energy landscape and the fluctuations of a model protein. The model exhibits two transitions, folding and dynamical transitions, when changing the temperature. The inherent structures corresponding to the minima of the landscape are analyzed and we show how their energy density can be obtained from simulations around the folding temperature. The scaling of this energy density is found to reflect the folding transition. Moreover, this approach allows us to build a reduced thermodynamics in the inherent structure landscape. Equilibrium studies, from full molecular dynamics (MD) simulations and from the reduced thermodynamics, detect the features of a dynamical transition at low temperature and we analyze the location and time scale of the fluctuations of the protein, showing the need of some frustration in the model to get realistic results. The frustrated model also shows the presence of a kinetic trap which strongly affects the dynamics of folding.

Figure 1

(Color online) Representation of the main structural elements of the $B_1$ domain of protein $G$ (drawn with VMD [30]). The purple part corresponds to the $\alpha$ helix, while the yellow ribbons are $\beta$ sheets. Residue number 1 is the end of the chain at the top right of the figure. The main structural elements of this protein are two $\beta$ strands comprising residues 1–7 and 14–19, which form a $\beta$ sheet denoted by ${\beta }_N$, one $\alpha$ helix made of residues 23–35, and two $\beta$ strands comprising residues 42–46 and 51–55, which form the $\beta$ sheet called ${\beta }_C$. Other residues belong to loops.

Figure 2

(a) Temperature evolution of the dissimilarity between the average structure of the protein model and the native structure. (b) Temperature evolution of the specific heat of the protein model.

Figure 3

Comparison of the density of states ${\Omega }_{\mathit{IS}}\left(e_{\alpha }\right)$ computed from the probability densities $P_{\mathit{IS}}(e_{\alpha },T)$ deduced from molecular dynamics trajectories at three different temperatures $T=T_f$, $T=0.83T_f$, and $T=0.69T_f$ for the unfrustrated model $U$ (top part) and the frustrated model $F$ (bottom part) $(d{e}_{\alpha }=0.02)$.

Figure 4

Top figures: density of state ${\Omega }_{\mathit{IS}}\left(e_{\alpha }\right)$ for the inherent structures, which is obtained from $P_{\mathit{IS}}(e_{\alpha },T_f)$ according to Eq. (9). The inset shows a magnification of the curve for the energies close to the ground state. The dotted lines are two distinct exponential functions with exponent corresponding to Eq. (12). Bottom figures: probability distribution (in logarithmic scale) of the inherent state energies at temperature $T_f$.

Figure 5

Thermodynamic quantities defined from the inherent structures. Top figure: conformational specific heat $C_{\mathit{IS}}∕k_B$. The arrows on the top axis point to the values of the temperatures $T_u$ and $T_s$ which enter in the expression of the density of inherent states [Eq. (12)]. Middle figure: magnification of the variation of the specific heat in the low temperature range, showing also the contribution $C_W=C_v-C_{\mathit{IS}}$, which is expected to reflect the properties of thermalization inside each well. Notice that the figure shows the value $C_W∕N_r$, i.e., the vibrational specific heat per residue, $C_{\mathit{IS}}$, which is expected to reflect a global property of the protein, is shown for the full molecule. The dashed line indicates the harmonic limit $C_W∕N_r=3k_B∕2$. Bottom figure: conformation entropy $S_{\mathit{IS}}$ and conformation energy $U_{\mathit{IS}}$. Inset: conformational entropy in logarithmic scale.

Figure 6

Temperature variation of the average fluctuations of the positions of residues for models $U$ and $F$, measured by $\Delta u^2$ [Eq. (19)].

Figure 7

Fluctuations of the position of the residues as a function of the index of the residue along the protein chain, at different temperatures, for model $U$ and model $F$. As indicated in the figures, residues with the lower indices correspond to the ${\beta }_N$ terminal part of the protein, residues with intermediate indices belong to the $\alpha$ helix, while the residues with the higher indices are part of the ${\beta }_C$ terminal region. The different curves show $\Delta u^2∕(T∕T_f)$ at temperatures $T=0.20–0.60$ by steps of 0.05, the curves being in order of increasing temperatures from bottom to top.

Figure 8

Temperature evolution of the fluctuations of the distances between different structural elements of the protein [defined by Eq. (20)] for model $U$ and model $F$. The residues belonging to each structural element are listed in the caption of Fig. 1.

Figure 9

(a) Variation of the fluctuations of the dissimilarity $D$ between the model conformations and the native structure averaged over a time $\tau$ vs $\tau$ for model $U$ at $T=0.7T_f$. The solid line shows the decay given by the law of large numbers. (b) Arrhenius plot of the relation time ${\tau }_r$ deduced from the fluctuations of the dissimilarity factor for models $U$ and $F$. The solid lines are fits to ${\tau }_r\propto \exp (E_B∕k_{B}T)$ with activation energies $E_B=6.2k_B{T}_f$ for model $F$ and $E_B=5k_B{T}_f$ for model $U$.

Figure 10

Example of the time evolution of the dissimilarity factor $D$ at $T=0.5T_f$. The presence of two types of fluctuations is visible: fast small amplitude variations and large changes occurring from time to time.

Figure 11

Density of inherent states of models $U$ and $F$ obtained at different temperatures, showing the difficulty to get a complete sampling of the phase space for model $F$ at low temperature.

Figure 12

(a) Histogram of the folding time of model $U$ and model $F$ at two temperatures. The large peak observed at low $T$ for model $F$ occurs because most of the initial states do not reach the native state within the simulation time of ${10}^8$ time units. (b) Time evolution of the folding of model $F$ at different temperatures. The time points $t_i$ selected to calculate $U_{\mathit{IS}}\left(t_i\right)$ evolve with a logarithmic scale according to $t_{i+1}=\sqrt{2}{t}_i$ and for each calculation the time interval $\Delta t_i$ for the integration is $\Delta t_i=(1-\sqrt{2})t_i$ so that the time domains which are analyzed are adjacent to each other.

Figure 13

Effective free energy $F\left(D\right)$ (dotted line) and inherent structure free energy $F_{\mathit{IS}}\left(D\right)$ (full line) of model $F$ at three temperatures in the vicinity of the folding temperature. The graphs have been shifted along the vertical axis to get the same values for the left minimum.

Figure 14

Inherent structure free energy $F_{\mathit{IS}}(D,T)$ as a function of the dissimilarity factor $D$ of the inherent structures, obtained from the analysis of folding pathways at different temperatures $T_2$ for model $F$.

Figure 15

Probability distribution of the inherent structures $P_{\mathit{IS}}^t(e_{\alpha },T_2)$ derived from the sampling of folding trajectories at temperature $T_2$ in a given time range for various folding temperatures $T_2$. In each case the top figure shows $P_{\mathit{IS}}^t(e_{\alpha },T_2)$ averaged for the whole duration of the simulation ${10}^8$ time units, except for the first ${10}^5$ time units which correspond to a short transient time during which the results might be influenced by the initial condition of each simulation. The lower figure shows $P_{\mathit{IS}}^t(e_{\alpha },T_2)$ in the last part of the folding trajectory ($0.7\times {10}^8<t<1.0\times {10}^8$ time units). (a) $T_2∕T_f=0.34$, (b) $T_2∕T_f=0.48$, (c) $T_2∕T_f=0.69$. A logarithmic scale is used for the probability distributions in all these figures.