Calculation of the Fraction of Native Contacts, Q, and Its Usage to Determine Folded Trajectories

Two residues are considered to form a native contact if their α carbons are less than 8 Å apart in the crystal structure. To account for thermal fluctuations in contact distances during simulation, a flexibility parameter Δ = 1.2 was used: a native contact between two residues is classified to be formed in a current frame of the simulated trajectory if their distance is shorter than 1.2 times the distance in the crystal structure. The fraction of native contacts, Q, was calculated for each protein during their posttranslational folding or refolding simulations. Only contacts between pair of residues i and j both within secondary structural elements as identified by STRIDE³⁷ and satisfying the criterion |i – j| > 3, where i and j are the residue indices, were considered; we excluded any secondary segment that is shorter than four residues from the analysis. To determine when a given trajectory of a protein is folded, we first characterized the fraction of native contact, Q, of each protein’s native state by performing ten 1.5 μs coarse-grained simulations at 310 K initialized from the native-state coordinates. The threshold for protein folding during refolding or posttranslational folding simulations, Q_threshold, was determined as Q_threshold = Q_mode^NS – 3σ, where Q_mode^NS is the average Q_mode over all 15 ns windows of the ten 1.5 μs native-state simulations and σ is the standard deviation of Q_mode^NS. To determine when folding occurred during refolding or posttranslational folding simulations, the mode of the Q values over a sliding 15 ns window was compared to the Q_threshold. A given trajectory is defined as folded if during its time evolution, Q_mode^15-ns ≥ Q_threshold, the folding time is the first time that the above condition is met.^35,38 The threshold value of Q for each protein is presented in Table 1.

Table 1

Threshold Value of Q of Three Proteins Computed from 10 Native-State Simulations Used to Determine if a Given Trajectory of Protein Folds

protein	Qthreshold = QmodeNS – 3σ
DHFR	0.9221
CAT-III	0.9269
DDLB	0.9521

Open in a separate window

Estimating the Folding Time of Slow-Folding Proteins with a Large Proportion of Unfolded Trajectories

When the portions of folded trajectories are less than 50% of total trajectories, it is not possible to estimate the folding time as the median first passage time.

We consider three-state folding kinetics with parallel pathways. State A folds rapidly to the native state N at the rate k₁, and state B folds slowly to the native state with a much smaller rate k₂ (k₁ k₂), and there is no interconversion between A and B. We have a set of ordinary differential equations respecting the rate of changing portion of states A and B

2

where [A] and [B] are the portion of non-native states A and B. The portion (survival probability) of non-native states at time t: S_U(t) = [A](t) + [B](t) = c₁ exp(−k₁t) + c₂ exp(−k₂t), where c₁ and c₂ are arbitrary constants. The initial condition that at time t = 0, the survival probability of non-native state = 1, we have S_U (t = 0) = c₁ + c₂ = 1, this yields: c₂ = 1 – c₁.

Hence, we computed the survival probability of the unfolded state as a function of time from simulations, and the resulting time series were then fit to the double-exponential equation

3

c₁, k₁, and k₂ are the fitting parameters. The time constants of the two kinetic phases are , with the larger of these two times determining the overall time scale of the folding process, τ₂ τ₁. To estimate the uncertainty of the folding time when fitting to double-exponential folding kinetics, we apply bootstrap resampling by randomly selecting trajectories from the list of simulations. We only consider the random sample with the coefficient of determination R² > 0.9. This procedure was applied to estimate the folding time of CAT-III.

Definition of the Progress Variable ς and Use to Monitor the Sequence of Pairs of Native Secondary Structure Elements Formed during the Folding Process

To account for the significant variation in folding times among different trajectories, we monitored folding pathways as a function of a progressive variable,³⁹ ς, defined as

4

where <···> indicates the average over all folded trajectories, and t_pair,i and t_fold,i are the folding time of pair and the whole protein folding time of the folded trajectory i, respectively. With this definition, we have 0 ≤ ς ≤ 1, ς = 0, which means that the pair under studied folds at the start of the simulation, and ς = 1 indicates the pair folds as the last step in the folding process. To determine the sequence of pairs of the secondary structure formation (defined in Figure S1 and Table S2), we consider a pair between two secondary structure elements that have more than one native contact. A pair is considered to be folded if its fraction of native contacts is larger than the threshold determined from native simulations. In our analysis of folding pathways, trajectories that did not fold within the 6 μs simulation duration were excluded.

Identifying Entanglement and the Changes in Entanglement

We use the approximation to the partial Gaussian double integration method proposed by Baiesi and co-workers⁴⁰ to calculate these partial linking numbers for a closed (loop) and open curve (termini). To identify lasso-like entanglements, we used the numerically invariant linking numbers,⁴¹ which describe the linking between a closed loop and an open segment in a three-dimensional space. This procedure is a modified version of the original protocol proposed by Baiesi to detect entanglement in coarse-grain protein structures. The original protocol is not computationally efficient to analyze trajectories since for each pair of contact, we have to calculate the linking number for all possible combinations of loop and threading segments. In our modified protocol, we only have to calculate the linking number between the closed loop (closes by native contact) and two tails. The closed loop is composed of the peptide backbone connecting residues i and j that form a native contact. Outside this loop is an N-terminal segment composed of residues 5 through i – 4 and a C-terminal segment composed of residues j + 4 through N – 5, where we exclude the first five residues of the N-terminal curve, the last five residues of the C-terminal curve, and four residues before and after the native contact to eliminate the error introduced by both the high flexibility and contiguity of the termini and trivial entanglements in the local structure; this metric is similar to whGLN.⁴² We can characterize the entanglement of each tail with the loop formed by the native contacts with two partial linking numbers denoted g_N and g_C. For a given structure of an N-residue protein, with a native contact present at residues (i, j), the coordinates R_l and the gradient dR_l of the point l on the curves were calculated as

5

where r_l is the coordinates of the C_α atom in residue l. The linking numbers g_N(i, j) and g_C(i, j) were calculated as

6

The total linking number for a native contact (i, j) is therefore estimated as

7

Comparing the absolute value of the total linking number for a native contact (i, j) to that of a reference state allows us to detect a gain or loss of linking between the backbone trace loop and the terminal open curves as well as any switches in chirality. Therefore, there are six changes in linking cases we should consider when using this approach to quantify entanglement (see Supplementary Figure S1 and Table 1 of ref (43)).

The degree of entanglement G is defined as

8

where (i, j) is the native contact in the crystal structure; NC is the set of native contacts formed in the current structure at time t; and g(i,j,t) and g^native (i,j) are, respectively, the total linking number of the contact (i, j) at time t and native structures estimated using eq 7. M is the total number of native contacts in the native structure and Θ is a Heaviside step function, equals 1 if the condition is true and equals 0 if the condition is false.

The difference between g(i,j,t) and G(t) is g(i,j,t), which is characterized by the entanglement in a given structure of contact (i, j) at time t, while G(t) provided information about the total number of contacts that changed the entanglement at time t.

Clustering and Coarse-Graining Conformational Space (Q, G)

The projection of conformation space onto (Q,G) reveals intermediate states that may be hidden when projected onto Q alone, as two states can have the same value of Q but one may be entangled while the other is not. Entanglement can prevent a protein from reaching its native state, as the loop-threading segment is improperly organized. Entangled states thus can form kinetic traps with large energy barriers preventing progression to the folded state, as large sections of the protein must unfold to allow disentanglement. To derive the log probability surface as a function of (Q,G), we first combined (Q, G) data from refolding and posttranslational folding for each protein and applied the Min–Max algorithm⁴⁴ for normalization. K-mean++ clustering⁴⁵ was then utilized to identify microstates, with 200, 400, and 400 clusters (microstates) being used for DHFR, CAT-III, and DDLB, respectively. As k-mean++ is a distance-based clustering algorithm, the normalization of data was necessary to prevent one-order parameter from dominating the distance measure. The resulting clusters were further coarsened into a small number of metastable states using the PCCA+ algorithms⁴⁶ to facilitate the interpretation of the folding pathways. The number of metastable states was determined based on the presence of a gap in the eigenvalue spectrum of the transition probability matrix; 11, 14, and 13 metastable states were used for DHFR, CAT-III, and DDLB, respectively. Both the clustering and coarse-graining processes were performed by using the PyEmma⁴⁷ and Deeptime⁴⁸ packages.

Identify Folding Pathways along the Order Parameters (Q, G)

To identify folding pathways from the simulated trajectories, the following procedure was followed:

• (1)

  For each discrete trajectory, the starting state of the first frame is added to the pathway.

• (2)

  The trajectory is then advanced, and the next state that differed from the last recorded state in the pathway was identified. If this state had not yet been recorded in the pathway, it was added to the pathway. If the state is already been recorded in the pathway, the pathway was truncated at the first instance of the recorded state and the trajectory was advanced from that point.

• (3)

  Repeat Step (2) until the end of the trajectory is reached.

This process resulted in pathways that contained no loops, and only recorded the on-pathway states for each discrete trajectory. The distribution of distinct pathways and the probabilities of transitioning from one state to another was then estimated based on the pathways of all of the discrete trajectories. The initial, folded, and misfolded states (in the folding/misfolding pathways plots) are colored yellow, blue, and red, respectively. A state is considered misfolded if there is a trajectory that becomes trapped in that state, and there is no direct transition to the native state. The size of the nodes is proportional to the probability of the state appearing in the coarse-grained trajectories. The size of the edges connecting the nodes is proportional to the number of transitions between states, and the red number beside the edge is the total number of transitions observed in the coarse-grained trajectories.

Protein Synthesis Does Not Increase the Folding Efficiency of CAT-III and DDLB

Using the same simulation protocol as DHFR, we performed refolding and posttranslational folding for CAT-III and DDLB proteins. In contrast to DHFR, the folding dynamics and population of folded trajectories for CAT-III and DDLB are relatively insensitive to posttranslational folding versus refolding. Specifically, for CAT-III, the log probability landscape of CAT-III is almost identical between posttranslational folding and refolding (Figure ​Figure22c). The progress of normalized Q of the folded trajectories is similar (Figure ​Figure22d), and the difference in the number of folded trajectories is insignificant (p-value = 0.14; Table 2). There are a large number of misfolded trajectories (Q < Q_threshold; Table 2) within the simulation time of 6 μs. The proportion of folded trajectories for CAT-III is less than 50%; we, therefore, estimated its folding time by fitting the survival probability of the unfolded state as a function of time to a three-state kinetic model (eq 3; see Materials and Methods section). There is no statistical difference in folding times for CAT-III between refolding (2.3 × 10⁵ ns, 95% CI [6.5 × 10⁴ ns, 1.7 × 10¹² ns]) and posttranslational folding (2.05 × 10⁵ ns, 95% CI [7.8 × 10⁴ ns, 1.6 × 10¹² ns]), p-value = 0.96 (Table 2).

In the case of DDLB, more than 50% of trajectories are folded; hence, the median folding time could be estimated. We find that the median folding time in refolding is 522.5 ns (95% CI [412.1 ns, 712.2 ns]), compared to the folding time in posttranslational folding, which is 426.3 ns ([264.7 ns, 690.9 ns]). We find that there is no difference in the median folding times or the number of folded trajectories between the refolding and posttranslational folding simulations (p-value = 0.87 for the number of folded trajectories and p-value = 0.18 for the median folding time comparisons; Figure ​Figure22f and Table 2). However, there are some observed differences: the log probability landscape in the posttranslational folding of DDLB sampled a smaller region along the Q coordinate, and the local minima were deeper compared to refolding (Figure ​Figure22e). This suggests that the cotranslational formation of native contacts may have occurred after translation.

To test the influence of simulation time on the results, the misfolded trajectories for CAT-III were extended to 30 μs and the misfolded trajectories for DDLB were extended to 15 μs. We find that only one additional trajectory each from the refolding and posttranslational folding simulations of CAT-III folds during this extended duration, at 15 and 29.2 μs, respectively. No misfolded trajectories of DDLB folded in either the refolding or posttranslational folding simulations. This suggests that these misfolded trajectories are kinetically trapped and unlikely to convert to the folded state at longer time scales—consistent with previously published results.¹¹

Measuring the Folding Mechanisms of Proteins Using Progress Variable ς Reveals the Differences for DHFR and Remains Robust for CAT-III and DDLB

Protein folding is typically thought to occur in a hierarchical fashion, with secondary structural elements first forming individually and then cooperatively coalescing into tertiary structures. With this in mind, we characterize the folding process of DHFR, CAT-III, and DDLB as the temporal sequence of formation of their stable pairs of native secondary structural elements with the aid of a progress variable, ς (see the Materials and Methods section, eq 4). The value of ς is relative to the time of complete folding of the protein, with ς = 0 indicating that the pair folds at the start of the simulation and ς = 1 indicating the pair folds as the last step in the folding process. To simplify the analysis, we restrict ourselves to pairs of secondary structures that have more than one native contact, as described in the Materials and Methods section and Table S2.

Based on this analysis, we observe a significant difference in DHFR. In posttranslational folding, all pairs of the native secondary structural elements belonging to the ABD domain fold cotranslationally (ς ~ 0), while in refolding, most of the pairs fold at the end of the folding process (ς ~ 1) (Figure ​Figure33a and Table 3). This suggests that the vectorial synthesis from the N-terminus to the C-terminus prevents the spontaneous cotranslational folding of some β-sheets in the C-terminal (C1, C2) and that the complete folding of DHFR occurs immediately upon release of the C-terminal from the ribosomal exit tunnel. These observations are consistent with previous experimental studies that have found that the central domain (ABD) acts as an independent folding unit during translation, while the DLD domain folds posttranslationally.³⁰ For CAT-III, the sequence of secondary structure pair folding is similar in both refolding and posttranslational folding, with all pairs folding late during the folding process (ς ~ 1; Figure ​Figure33b and Table 3). For DDLB, the overall folding order is similar, but some differences were observed, such as in posttranslational folding, four pairs in the center domain (C13, C19, C22, and C24) fold cotranslationally (ς = 0), two pairs in the N-terminal domain (C7, C8) fold posttranslationally but before the complete folding occurs (ς ~ 0.65; Figure ​Figure33c and Table 3), while these pairs fold at the last event in refolding. Thus, protein synthesis and posttranslational folding do not significantly perturb the folding mechanisms of CAT-III and DDLB.

Open in a separate window

Figure 3

Comparisons of folding processes of DHFR, CAT-III, and DDLB in posttranslational folding versus refolding are shown as temporal sequences of secondary structure pairs formed over time with the aid of a progress variable ς. (a) Folding mechanism of DHFR is significantly different: all pairs in the ABD domain fold cotranslationally in posttranslational folding simulations, (b) CAT-III: there is no difference between posttranslational folding versus refolding, and (c) DDLB protein exhibits a small difference in four pairs of the center domain (C13, C19, C22, and C24) and two pairs (C7, C8) in the N-terminal domain. Refolding and posttranslational folding data are represented by black circles and red stars, respectively.

Native Entanglements Exist in the Crystal Structure of CAT-III and DDLB Proteins

We hypothesized that there is something distinct about the native topologies of CAT-III and DDLB that leads to a large proportion of misfolding. Several recent papers have predicted a link between misfolding involving a change in the entanglement status and long-lived misfolded states,^11,35 including the failure to form native entanglements. Indeed, this is the molecular hypothesis explaining the observation that experimental folding rates of proteins decrease as the number of times the threading segment pierces the loop increases.⁴⁰ To further understand this phenomenon, we investigate whether entanglement may play a role here by calculating the degree of entanglement for these proteins using eq 7.

We find that the crystal structure of DHFR does not contain any entanglements. In contrast, CAT-III has 16 native entanglements, with 14 of them consisting of a loop located near the N-terminus and a threading segment at the C-terminus. The remaining two native entanglements have a loop located near the C-terminus and a threading segment at the N-terminus. Similarly, DDLB has 36 native entanglements, half of which consist of a loop located closer to the N-terminus and a threading segment at the C-terminus, while the other half has a loop located closer to the C-terminus and a threading segment at the N-terminus. Representative examples of these entanglements are shown in Figure ​Figure44a,b for CAT-III and DDLB, respectively. Furthermore, proteins with an entanglement loop closer to the N-terminus were found to be folded more difficultly than the proteins with a loop closer to the C-terminus.⁵⁰ This explains why DHFR (without native entanglement) can fold easily and small portions of CAT-III trajectories (most entanglement loops are located near the N-terminus) folds in our simulation. This observation suggests that entanglement plays an important role in the proper folding of proteins.

Open in a separate window

Figure 4

Example of native entanglements in the crystal structures of CAT-III and DDLB. The closed loop and crossing section of the threading segment of their entangled regions are colored red and blue, respectively. The loops are closed by noncovalent contacts between two residues (colored yellow), and the rest part of the protein is colored gray. (a) Representative native entanglement in CAT-III: the loop (colored red) is closed by a native contact between residues 8 and 77, and the threading segment consists of residues 177–208. (b) Representative native entanglement in DDLB: the loop (colored red) consists of residues 98–146, and the threading segment consists of residues 160–184.

Protein Synthesis Assists the Folding of DHFR by Avoiding Misfolded States with Non-Native Entanglements

Entanglement plays an important role in the proper folding of proteins. To further characterize the folding pathways of DHFR, we clustered the conformational space based along the order parameters Q and G and then assigned them to metastable states (see the Materials and Methods section). In posttranslational folding, DHFR can spontaneously fold to its native state once the C-terminus is released from the ribosome. The two-dimensional log probability surface is concentrated in the region around the folded state (small G, high Q; Figure ​Figure55b), which is consistent with our 1D log probability landscape from the previous section. Specifically, the posttranslational folding simulations of DHFR only sample two states, 5 and 10 (the folded state). There are no misfolded trajectories in posttranslational folding, as all trajectories reach the folded state at the end of the simulation. The protein cotranslationally folds to the ensemble state 5, which has about 60% of native contacts formed (the fraction of native contacts with non-native entanglement is negligible, around 0.16%), and the folding process simply involves diffusion to the folded state (state 10). Folding network analyses reveal that 100% of folding pathways go straight from the initial state 5 to the folded state 10 (Figure ​Figure55d). There is no off-pathway state in the posttranslational folding of DHFR.

Open in a separate window

Figure 5

Ribosome helps DHFR fold more efficiently. (a, b) −ln(P) surface in refolding and posttranslational folding, respectively, where P is the probability of sampling particular Q and G values. The centers of metastable states and their corresponding indices are shown on top of the surface (black points). (c, d) Transition network from discrete trajectories of refolding and posttranslational folding simulations. The yellow, red, and sky blue nodes correspond to the initial, misfolded, and folded states, respectively. The black numbers on the nodes match the indices of metastable states in panels (a) and (b). The red numbers beside the edges indicate the number of direct transitions between states observed in discrete trajectories.

Refolding from the thermally unfolded ensemble is more complicated, compared to posttranslational folding. The −ln(P) surface has sampled a broad region in the non-native (low Q) or near-native (high Q) regions. We found that the population of DHFR refolding samples had a large number of entangled states, indicated by high values of G (Figures ​Figures55a and S3). The protein follows two parallel pathways to reach the native state: we find that the dominant pathway (*→ 5 → 10), which is the only pathway observed in posttranslational folding, accounts for 87% of the total trajectories in refolding simulation and a small portion (four trajectories, accounts for 4% of total trajectories) folds via intermediate state 7 (*→ 7 → 10). In addition, we find that 9% of trajectories become trapped in misfolded states (states 7–9). The broader −ln(P) surface in refolding is caused by a small number of misfolded trajectories. Five trajectories become trapped in state 7, three trajectories become trapped in state 8, and one trajectory becomes trapped in state 9. States 8 and 9 are off-pathway misfolded states, as we do not observe any folding events (conversion to the folded state 10) if the protein visits these states. When the protein samples the near-native state 7, only 40% of trajectories can fold successfully (* → 7 → 10/folded), while the remaining 60% fold to misfolded states (Figure ​Figure55c).

M.S.L. acknowledges that this work was supported by the National Science Centre, Poland (grant 2019/35/B/ST4/02086). E.P.O. acknowledges support from the National Science Foundation (MCB-1553291) as well as the National Institutes of Health (R35-GM124818). This research was supported in part by the TASK Supercomputer Center in Gdansk and PLGrid Infrastructure in Poland (Prometheus and Ares supercomputers).