Protein folding from heterogeneous unfolded state revealed by time-resolved X-ray solution scattering.

One of the most challenging tasks in biological science is to understand how a protein folds. In theoretical studies, the hypothesis adopting a funnel-like free-energy landscape has been recognized as a prominent scheme for explaining protein folding in views of both internal energy and conformational heterogeneity of a protein. Despite numerous experimental efforts, however, comprehensively studying protein folding with respect to its global conformational changes in conjunction with the heterogeneity has been elusive. Here we investigate the redox-coupled folding dynamics of equine heart cytochrome c (cyt-c) induced by external electron injection by using time-resolved X-ray solution scattering. A systematic kinetic analysis unveils a kinetic model for its folding with a stretched exponential behavior during the transition toward the folded state. With the aid of the ensemble optimization method combined with molecular dynamics simulations, we found that during the folding the heterogeneously populated ensemble of the unfolded state is converted to a narrowly populated ensemble of folded conformations. These observations obtained from the kinetic and the structural analyses of X-ray scattering data reveal that the folding dynamics of cyt-c accompanies many parallel pathways associated with the heterogeneously populated ensemble of unfolded conformations, resulting in the stretched exponential kinetics at room temperature. This finding provides direct evidence with a view to microscopic protein conformations that the cyt-c folding initiates from a highly heterogeneous unfolded state, passes through still diverse intermediate structures, and reaches structural homogeneity by arriving at the folded state.

Proteins should retain unique three-dimensional (3D) structures to properly express its biological functions in living cells. The linear sequence of amino acids in a protein released from the ribosome encodes not only its native structure but also the folding mechanism by which this structure is attained (1). According to Levinthal’s argument (2), an ensemble of unfolded proteins heterogeneously populated with various conformations will not likely find its native state through a random search in the vast conformational space but rather recovers its native form through well-defined folding pathways.

Numerous theoretical studies have been conducted to describe the complex folding dynamics in view of free-energy landscapes (3–10). Among various conceptual frameworks, it has been suggested that the funnel-like free-energy landscape with many microstates is a scheme pertinent to describing the protein folding (3, 4, 7, 10–12). In this picture, a key factor determining dynamics of protein folding is the native bias. This refers to the dominance of interactions stabilizing the specific native conformation over nonnative interactions that may guide to topologically distinct traps or frustrated conformations, so-called misfolded states (3–5, 13). Another factor is a conformational heterogeneity in the unfolded state and on-route intermediates, which can potentially alter the kinetic behavior in the protein folding due to the possibility of broadly distributed pathways. Thus, this folding model has an advantage for the explanation of conformational heterogeneity and its energetics along the reaction pathway of protein folding. To verify the scenario of the funnel-like free-energy landscape in protein folding, numerous experiments with diverse techniques such as laser-induced temperature jump, rapid mixing, and single-molecule spectroscopy have been conducted (14–24). Nevertheless, comprehensively understanding protein folding, especially in regard to its global conformational changes paired with the structural heterogeneity, has been elusive. This is simply due to the complex nature of any protein. Observing temporal evolutions of protein molecules with such intrinsic diversity can be ideally achieved with computational means such as molecular dynamics (MD) simulations. At present, however, reliably accomplishing the goal is not a trivial task at all. This is partly because simulating a long enough time scale is still a large burden. Perhaps more importantly, commonly available force field parameters are usually optimized for folded structures and may not equally behave well during the entire course of folding events from diverse unfolded configurations.

The elusivity can be well exemplified with cytochrome c (cyt-c). It is a small metalloprotein with 104 amino acids and a heme group that is covalently bound to the protein backbone. Since the pioneering study of cyt-c folding initiated by the dissociation of carbon monoxide from the heme group (25), cyt-c has been recognized as an appropriate model system to probe the folding of single-domain protein in a time-resolved manner. There have been many experimental studies that tracked the cyt-c folding in the wide time region from submicroseconds to late milliseconds (18, 26–36). These experiments revealed that the self-assembly of protein backbones concomitant with the collapse of extended conformation begins on submillisecond regime. Further, a maximum rate of protein folding was carefully scrutinized by using the optical triggering to refold cyt-c, and the experiment showed the diffusion rate of (∼35 μs)−1 for the local contact of two sites separated by ∼50 residues (37). If we assume that the local contacts in the early process of protein folding correctly reflect the native interactions, we can expect a rapid transition to the native state through the nucleation-collapse mechanism (38). However, the time-resolved experiments showed that the actual folding time in cyt-c takes several tens to hundreds of milliseconds (18, 26–34), which are much slower than the diffusion time of protein backbone reported from an earlier study (37). To explain the kinetic difference between the fast intrachain motion of the protein backbone and the relatively slow folding, the slow phase event on the millisecond time region has been obscurely presumed as either the sequential transitions through successive intermediates (26–29, 31–33) or the detours via frustrated conformations on the route to the native conformation (18, 38). This kinetic complexity during cyt-c folding has been presumed to originate from the complication of folding pathways associated with the ligations by natively or nonnatively contacting residues to the heme group. By employing imidazole binding to the heme group or by protonating misligating histidine residues at low pH toward suppressing the nonnative contact formation, some studies showed that the two-state folding mechanism was operational (39, 40), whereas the absence of imidazole accompanied the additional ligand exchange by methionine (26, 28, 29, 33). Meanwhile, other studies showed deviation from the two-state kinetics irrespective of the suppression of nonnative coordination (41, 42). In the end, even for this rather simple cyt-c with a single domain, there is no consensus on the interpretation of the folding dynamics in view of its conformational heterogeneity.

Time-resolved X-ray solution scattering (TRXSS) may serve as an effective way of obtaining a comprehensive view, mainly owing to the superb structural sensitivity from X-ray scattering. Because the redox state of heme acting as a cofactor is coupled with the folding free energy, it is possible to find a denaturing condition where the protein is fully unfolded in one oxidation state while it is fully folded in the other state. In a certain denaturing condition, the change in the redox state of the cofactor resculpts the free-energy landscape, and this brings about a strong native bias toward the folded state as shown in Fig. 1. In this work, we combined TRXSS, also known as time-resolved X-ray liquidography (43–51), with the external electron injection method in order to investigate the redox-coupled folding dynamics of cyt-c. We focus on examining 1) how the structurally heterogeneous nature of the unfolded state changes during the course of folding and 2) what type of kinetic behaviors are accompanied and linked with the inherent heterogeneous nature. The ensemble optimization analysis for the X-ray scattering data based on MD simulations shows that the initial oxidized state which contributes to the folding dynamics has a heterogeneously populated ensemble with diverse unfolded conformations. During folding, this converts to a narrowly populated ensemble around a well-organized folded conformation through an intermediate state, which still contains conformationally heterogeneous characteristics. Surprisingly, the slow folding phase of cyt-c accompanying the overall structural change manifests the stretched-exponential behavior even at room temperature. The finding from the analysis of X-ray scattering data provides direct evidence for the protein folding from the heterogeneously populated unfolded state through the multiple kinetic pathways and insights into the fundamental principles governing protein folding in view of global protein conformations and the related kinetic behaviors.

Fig. 1.

Free-energy landscapes induced by external electron injection. Conventional free-energy surfaces for two different conditions as a function of representative reaction coordinates ( and ). The reaction coordinates correspond to local structural changes in protein such as vibrational modes of specific residues or distances between residues. The free-energy surface can be resculpted by the change in protein environment, indicated as the yellow arrow, such as oxidation state of cofactor and solvent condition. From this, the basin of the unfolded state in the free-energy surface (Bottom) changes as a hillside in the resculpted energy surface (Top). The change of energetics activates the spontaneous folding pathways, indicated as the magenta arrow, along the reaction coordinates. As indicated with the gray dotted arrow, the pathway is inhibited in the initial condition.

Results

Structural Analysis of X-Ray Scattering Curves.

Prior to analyzing the TRXSS data, we first investigate the static X-ray scattering data for both the ferric state (oxidized) and the ferrous state (reduced) in order to confirm that ferric cyt-c is unfolded at the denaturation condition of 3.5 M GdnHCl and to investigate the heterogeneity of unfolded state. To overcome the limitation of adopting a single conformation as commonly employed in analyzing the structure of a well-structured protein (45–47, 52), we applied the ensemble optimization method (EOM) combined with MD simulations (53–55). Briefly, protein structures describing the conformational landscape are first generated from MD simulation snapshots and are used to calculate theoretical X-ray scattering curves corresponding to the individual conformational states. In this analysis, the experimental X-ray scattering curve is assumed to arise from an ensemble including a certain number of coexisting conformational states. This is well suited for describing structural features of a potentially flexible system (53, 55). The schematic of this analysis and its details are described in (see , for more information). By using a genetic algorithm (56), the ensemble that best fits the experimental data is finally chosen for extracting the optimal structural parameters such as the radius of gyration (Rg) and the 3D protein conformation. To guarantee wide enough conformational sampling in MD simulations, we employed three different simulation times of 100 ns, 500 ns, and 1 μs and checked the dependence of the final EOM results on the time. Specifically, we repeated the entire EOM analyses using the pools generated from six independent MD trajectories with each simulation time.

Fig. 2 shows the distributions of Rg derived from the EOM analysis for the static X-ray scattering curves of the oxidized and the reduced states based on different simulation times. The primitive Rg distributions directly from MD simulations before EOM are also shown in Fig. 2 (). Even though the maximum Rg values from different pools are rather different, all of the pools cover the Rg values (24–29 Å) of unfolded cyt-c reported from previous experiments (30, 57). The results from the EOM analyses show that while the Rg distributions of the oxidized state optimized with six independent MD simulations quite disagree with each other with 100-ns simulation time, they eventually well overlap after 500 ns. This simply demonstrates that the conformational distribution in view of Rg from the MD simulation is converging with the increase in the simulation time (). Even though there are some marginal differences in the optimized ensembles from the MD simulations of 1-μs duration, it is obvious that the various unfolded conformations are heterogeneously populated in the oxidized state under 3.5 M GdnHCl. The average Rg value from the EOM based on the 1-μs-long simulation changes from ∼26 Å for the oxidized state to ∼13 Å for the reduced state (). Based on the representative protein conformations determined from the EOM analysis (), we observe that at least four unfolded conformations with both extended and entangled structures were needed for describing the ensemble of oxidized cyt-c. According to a small-angle X-ray scattering (SAXS) study of oxidized cyt-c at differing GdnHCl concentrations (57), an equilibrium between native, partially unfolded, and fully unfolded states exists in the denaturing condition, and at 3.5 M GdnHCl the molar ratio of the three is ∼0:4:6. In addition, based on the existence of an additional singular vector from the singular value decomposition (SVD) analysis, it was vaguely presumed that one more denatured state may exist. Our model of oxidized cyt-c determined from EOM with at least four representative states has a higher heterogeneity than was proposed in this previous SAXS study (57). For the ensemble of reduced cyt-c, only one folded conformation was enough indeed. We note that the theoretical scattering curve from one individual representative structure in the ensemble of oxidized cyt-c cannot reproduce the experimental data (). This aspect implies the inherent conformational heterogeneity in the oxidized cyt-c under 3.5 M GdnHCl.

Fig. 2.

Structural analysis of X-ray scattering curves based on EOM. (A–C) Distributions of Rg determined from the EOM analysis for the static X-ray scattering curves of the oxidized and the reduced states based on unfolding MD simulation times of (A) 100 ns, (B) 500 ns, and (C) 1 μs. For each simulation time, a total of six independent trajectories were sampled by starting from one common structure. Each panel shows the Rg distributions of finally optimized ensembles (squares) as well as the conformational space coverage by the candidate protein structures from the MD simulations (dotted lines). The theoretical scattering curves of the final ensembles well fit the experimental curves within the error bars (). As the simulation time for sampling the candidate structures increases, the discrepancy among the six unfolded distributions optimized with the six trajectories becomes inconspicuous. (C, Bottom) The representative structures of folded and unfolded state with those relative populations. It is clear that various unfolded conformations are heterogeneously populated in the optimized unfolded ensemble (open squares) and single conformation is homogeneously populated in the optimized folded ensemble (closed squares).

TRXSS.

The reaction scheme for the folding of equine heart cyt-c triggered by external electron transfer is shown in Fig. 3. Upon 355 nm excitation of nicotinamide adenine dinucleotide (NADH) acting as the external electron source, a reactive species with solvated electron is generated that rapidly reduces the oxidized cyt-c with the unfolded protein conformation (26–29, 31, 33, 58). Based on the equilibrium unfolding curves for oxidized and reduced cyt-c (Fig. 3 and ), at neutral pH the concentration of guanidine hydrochloride (GdnHCl) in the range from 2.8 to 4.5 M induces the reduced cyt-c to be fully folded, whereas the same condition drives the oxidized form to the unfolded state with a fraction higher than 90%. This is in good agreement with the results reported from earlier studies (29, 33–35). Therefore, the concentration of 3.5 M GdnHCl is a proper condition to initiate spontaneous folding of unfolded cyt-c through the change in the redox state of the cofactor. Using the oxidized cyt-c dissolved in the phosphate buffer with 20 mM NADH and 3.5 M GdnHCl, we performed TRXSS experiments in order to investigate the folding dynamics triggered by the external electron transfer.

Fig. 3.

Schematic of folding dynamics of cyt-c induced by external electron injection. (A) Unfolded cyt-c with the oxidized heme group undergoes folding dynamics triggered by the electron injection from NADH acting as an external electron donor to the heme group in cyt-c. (B) Equilibrium unfolding curves for oxidized (Fe3+) and reduced (Fe2+) cyt-c as a function of molar concentration of GdnHCl. These curves were reconstructed from the circular dichroism spectroscopic data. The TRXSS experiment was performed by employing the photoinduced electron transfer under the concentration of 3.5 M GdnHCl.

Time-resolved difference scattering curves, q∆S(q, t), were measured in a wide time region covering from 31.6 μs to 316 ms (). The difference scattering curves show not only a significant amplitude in the small-angle region (q < 0.4 Å−1), which is related to the global conformational change, but also oscillatory features in the wide-angle region (0.4 Å−1 < q < 2.2 Å−1), which can be attributed to the subtle structural change such as the rearrangement of the secondary structure (45, 47, 50) or the thermal response of solvent with the relaxation of the photoexcited molecule (45, 51). Also, from previous studies (27, 58), it is known that photoexcited NADH forms photoproducts such as (NAD)2 as well as it recovers into its ground state. To discriminate the contributions by solvent heating and (NAD)2 formation from the TRXSS data, we implemented an additional experiment using the buffer solution with the same composition but without the protein (see , for details). In principle, the formation of (NAD)2 and the accompanying increase in electron density may alter the scattering signal in the small-angle region (q < 0.4 Å−1). The difference scattering signal of the protein-free buffer solution in the small-angle region is small and featureless, whereas the signal in the wide-angle region (0.8 Å−1 < q < 2.2 Å−1) is dominant (). From the comparison of the difference scattering curves of the two solutions (), we can verify that the effect of the formation of photoproduct on the entire q region is negligible and that only the response to solvent heating is prominent in the wide-angle region. Therefore, to extract the photoinduced signal of pure protein from the time-resolved X-ray data, the contribution of solvent heating was removed by subtracting the signal that was obtained from a separate experiment on the protein-free solution by following the established protocols (45, 47, 51).

Kinetics of Protein Folding.

Fig. 4 shows the time-resolved X-ray difference scattering curves, q∆S(q, t), which are the heating-free experimental data. Quantitative kinetic analysis was conducted by using SVD to determine the rate constants that describe the time-dependent change of the heating-free data (45–47, 59). SVD decomposed the original data into time-independent components (left singular vectors), their associated time dependence (right singular vectors [RSVs]), and the relative contribution of each component (singular values). From this analysis (), we identified two significant singular vectors and the corresponding autocorrelation factors, suggesting that two species discernible by X-ray scattering are formed in the measured time range (Fig. 4). The simplest two-state kinetics involving two species will be a single-exponential conversion of the first species to the other one, and globally fitting with this assumption failed to give an acceptable fit (the case with β = 1 in Fig. 4). This deviation points toward the existence of parallel pathways connecting the two species. We considered two candidates for two-state kinetics with parallel pathways: 1) bifurcated kinetics from the first to the second species with two exponentials, which can be considered as the simplest parallel kinetics, and 2) stretched exponential kinetics with , as an extreme case with parallel pathways, in which the β value stands for describing the extent of superposition of many exponential decays (60). Both the biexponential with timescales of 16.4 (±5.1) ms and 193 (±34) ms and the stretched exponential with τ = 185 (±75) ms and β = 0.70 (±0.08) show good agreements with RSVs (). We note that the bifurcated kinetics is different from the consecutive sequential kinetics suggested in earlier optical spectroscopic studies (26, 28, 29), where three species were required to explain two time constants. Because the SVD analysis suggests that only two species formed in the measured time region, we ruled out the possibility of this sequential kinetics ().

Fig. 4.

TRXSS data for the folding dynamics of cyt-c. (A) Experimental (black) and theoretical (red) difference scattering curves. (B) SADSs extracted from the global kinetic analysis. The second SADS (red) is compared with the difference for the static X-ray scattering curves of reduced and oxidized forms (gray). The high similarity between the second SADS and the static difference, q∆Sstatic, confirms that the observed dynamics is linked with the protein folding. (C) Time-dependent population change of the intermediates (dots) is compared with the theoretical population dynamics (lines) based on the stretched exponential kinetics. The values shown on the right side of graph represent the β values in the stretched exponential, , ranging from 0.15 to 1.0. The experimental population dynamics is in good agreement with the stretched exponential kinetics with the time constant (τ) of 185 (±75) ms and the β value of 0.70 (±0.08).

On the basis of both stretched exponential and bifurcated kinetics models, we performed principal component analysis using the time constants determined from the SVD analysis. The details on this kinetic analysis are described in (see , for more information). From this analysis, we extracted the species-associated difference scattering curves (SADSs) and time-dependent contributions of SADSs. Those determined from the stretched exponential kinetics are displayed in Fig. 4. Each SADS contains information on the structure of a protein intermediate. The first SADS formed in the microsecond regime has the monotonous feature with a small negative amplitude around 0.1 Å−1 in the difference scattering curve, whereas the second SADS formed in the late millisecond regime has a prominent oscillatory feature in the q space from 0.05 to 0.8 Å−1. Even though the first SADS has a much smaller contribution than the second one, the omission of the first species in the folding kinetics results in a worse agreement between the experimental and the theoretical data (), indicating that the cyt-c folding should involve the first species. The small decrease of the scattering intensity in the small-angle region of the first SADS implies that the protein conformation of the first intermediate is slightly more unfolded relative to that of the initial state. To identify the origin of the conformational change in the transition from the first SADS to the second one, we compared the difference of the static scattering curves of the reduced and the oxidized forms [∆Sstatic(q) = Sreduced(q) – Soxidized(q)] with the appropriately scaled second SADS as shown in Fig. 4. The good agreement between them indicates the formation of the folded state through the transition from the first to the second species. When the scattering curves were reconstructed from the linear combination of SADSs obeying the stretched exponential kinetics, they exhibited excellent agreements with the experimental curves (Fig. 4), confirming that cyt-c folding faithfully follows nonexponential kinetics. In fact, almost the same SADSs are also obtained from the bifurcated kinetics model because both cases are indistinguishable in the trace (). Therefore, the choice between the stretched exponential and the bifurcated kinetics in explaining the cyt-c folding dynamics cannot be made solely based on the kinetic trace. For this reason, we implemented an additional analysis based on the maximum entropy method (MEM) whereby the time-dependent signal can be explained by a distribution of time constants (61). From this analysis, the optimal distribution of log-time constants was extracted from the experimental data (see , for details). The broad distribution of log-time constants covering from the early to the late millisecond regime supports that the experimental trace likely follows the stretched exponential kinetics rather than the bifurcated one. The details about the folding kinetics in this regard are given in .

Discussion

The kinetic analysis shows that the TRXSS data can be explained by the following scheme for the folding dynamics induced by the redox state change of the cofactor:Here U is the initial unfolded state prior to the photoreduction of cyt-c, and U′ and F are the intermediate and the final photoproduct involved in the folding dynamics, respectively. The results from the EOM analysis for the static scattering curve of oxidized cyt-c shown in Fig. 2 suggest that the initial U state with an oxidized form of the iron-containing heme group has the conformational heterogeneity including various unfolded structures. The U′ intermediate is formed within 31.6 μs, which is the first time delay of our experiment. Considering the time of ∼33 μs taken to reduce the oxidized cyt-c with solvated electron and NAD radical produced from NADH photoexcitation (27, 58), the U′ state likely possesses a reduced form of the iron-containing heme group. In addition, the decrease in scattering amplitude in the small-angle region of the first SADS indicates that its volume and size slightly increases in comparison with the U state. From the comparison of the second SADS and the difference between static X-ray scattering curves of the oxidized and the reduced cyt-c, we find that the final F state has the completely folded conformation with well-organized secondary structures and reduced heme.

As mentioned already, based on the global kinetic analysis of TRXSS data, the main folding event in the transition from U′ to F can be explained by both stretched exponential and bifurcated kinetic models. This is from a typical uncertainty problem in the experimental data with a limited signal-to-noise ratio. While the MEM results favor the stretched exponential kinetics over the bifurcated one, having additional logical arguments other than the fitting quality in the kinetic analysis will be desirable. To address this, we consider the structural information available from the EOM analysis for the static scattering data. The ensemble of oxidized cyt-c linked to the initial U state shows the prominent structural heterogeneity including at least four representative conformations with both extended and entangled structures (Fig. 2 and ). If we apply the stretched exponential kinetics, the different unfolded conformations with the Rg values of 20 to 30 Å in the U state individually undergo reduction to the structurally heterogeneous U′ state, which is kinetically expressed as one representative species. After this, various unfolded conformations in the U′ state will likely adopt different folding pathways with different folding times toward the homogeneous folded state F with Rg of 13.2 Å. The superposition of many parallel pathways with the different timescales is compatible with the stretched exponential behavior. When we assume bifurcated kinetics, two parallel pathways may be attributed to the folding event from the heterogeneously populated U′ state to the folded state. As an example, some portion of the unfolded structures may possibly rapidly fold without a detour to a misfolded state, and the other remaining portion may progress through a misfolding pathway toward the native state. Namely, the two different timescales can originate from the presence or absence of the misfolded state in the folding pathway. Although this hypothetical situation may support the bifurcated kinetics, it is unclear how the extremely heterogeneous conformations in the U′ state will favor only two pathways toward the final folded state. It will be more persuading that the heterogeneous unfolded conformations adopt multiple pathways that potentially pass through various misfolded conformations on the free-energy surface. The superposition of such various folding paths actually merges to a single stretched exponential kinetics, and therefore, we suggest that the stretched exponential will be the operational kinetics in the cyt-c folding.

Previous studies of cyt-c folding suggested the presence of kinetic phases in the microsecond time regime called as burst phase (). Using time-resolved circular dichroism (TRCD) spectroscopy on electron-transferred cyt-c folding, Kliger and coworkers reported the time constant of 5 μs in the burst phase and attributed it to the rapid formation of partial secondary structures linked with the collapse of extended random coil (33). This earliest event was similarly observed with time-resolved optical rotatory dispersion (TRORD), where the fast phase process was assigned to the misligation of His33 (26, 28, 29). These studies suggested that the early-phase intermediate generated with the sudden collapse of random coil converts directly to the native state during the slow phase in hundreds of milliseconds. In terms of timescales, the U′ state observed in our study can be compared to the misligated state accompanying the partial formation of the secondary structure observed by TRCD and TRORD (26, 28, 29, 33). The first SADS, however, exhibits a decrease of the scattering amplitude in the small-angle region, which implies a slight increase of the protein volume and size in comparison with the U state. In other words, this slightly expanded U′ state, occurring before the rate-determining folding, is structurally distinguished from the molten globule state that possesses more partially folded conformations formed by the sudden collapse of the random coil (30, 41, 62). This difference is likely a consequence of adopting a strong denaturing condition with 3.5 M GdnHCl. According to an earlier fluorescence study (41), the relative amplitude of the early-phase kinetic trace of fluorescence significantly decreases with an increase in the concentration of GdnHCl. This is natural as GdnHCl tends to destabilize compact states that may form in the initial stage of folding (30, 57). In addition, destabilization of higher charged heme after oxidation in the low dielectric protein environment may be mitigated at high GdnHCl concentrations (63), and the reduction of the heme may not energetically prefer the early-phase protein motion necessary for forming the molten globule state such as the contraction of the backbone chain.

The kinetic analysis of TRXSS data shows that the U′ state undergoes transition to the final state with the time constant of 185 ms. In the TRCD study with multiexponential fitting, the kinetics in this time regime showed a decay with a time constant of 110 ms, which was attributed to the displacement of misligated residues at the heme cofactor (33), although the uncertainty is rather large due to the limited signal-to-noise ratio in the TRCD data. Using TRORD spectroscopy, a single decay of 160 ms in the slow phase was observed at a slightly lower GdnHCl concentration of 3.3 M and was assigned to the formation of the native structure with the ligation process between Met-80 and the heme (29). A similar folding time of 185 ms was reported for tuna heart cyt-c, containing a tryptophan residue at the site occupied by His33 in equine heart cyt-c (28). With transient absorption spectroscopic experiments, it was shown that the slow phase folding of reduced cyt-c follows a single exponential decay with a lifetime of hundreds of milliseconds that coincides with the change in tryptophan fluorescence (31, 32). Our time constant of 185 ms in the stretched exponential kinetics identified by TRXSS relatively well matches those from the previous studies (). However, we note that the heterogeneous transition of U′ → F expressed as the stretched exponential is in stark contrast to the normal exponential behavior interpreted as sequential transitions of intermediates with the conformational homogeneity in earlier studies (26, 28, 29, 31, 32). On the other hand, a stopped-flow mixing study for the reduced cyt-c combined with spectroscopic measurements observed two phases in the refolding kinetics when the GdnHCl concentration was smaller than ∼1.5 M, although at 3.5 M GdnHCl, which is the same concentration used in our work, only a single phase was observed (42). The two phases could be explained by either the conformational heterogeneity of unfolded cyt-c or the kinetic partitioning of the homogenous unfolded state (64). Our observation on the global conformational changes supports the former scenario even at the high GdnHCl concentration condition. This distinction is possible as TRXSS directly detects the overall structure unlike typical spectroscopic tools.

Indeed, the heterogeneous transition to the final folded state in slow phase folding of cyt-c may represent the experimental epitome of the funnel-like free-energy landscape in protein folding, which assumes multiple transitions from different microstates in the unfolded state to a unique global minimum of the folded state (12). A plausible scenario consistent with our observation is that the different unfolded structures in the initial oxidized state adopt conformation-specific pathways and stop by on-route intermediate configurations located at different positions of the free-energy landscape. Specifically, at an early stage of folding, some fractions of unfolded protein structures will be forming native interactions and are led to the correctly folded state, whereas others will face trap sites at local minima on the folding energy landscape. The differences in folding pathways originating from the differences in the unfolded conformations will likely bring about different folding times for each folding pathway. The stretched exponential kinetics observed with cyt-c in the slow phase can be interpreted as the sum of many single exponentials, and this aspect supports the kinetic partitioning of the heterogeneous unfolded state in folding, which may arise from the rugged features around the high-energy region of folding free-energy landscape. Our scenario for the folding dynamics of cyt-c can also be connected to the highly stretched exponential behavior in the early phase of folding, reported from a low-temperature stopped-flow study for the refolding of ferric cyt-c (18). As noted earlier, this behavior will correspond to the many parallel pathways initiating from as many different unfolded conformations. Thus, our results provide clear evidence for the heterogeneous folding of cyt-c involving structurally distinct multiple unfolded conformations and are in contrast to the view with unidirectional sequential pathways from a homogeneously defined unfolded state discussed in earlier studies (26–29, 31, 33, 34). A previous SAXS study reported that the high concentration of GdnHCl increased the population of fully unfolded state involving a high content of random coils (57), and another spectroscopic study showed that as the concentration of GdnHCl increased, folding events on the timescale of hundreds of milliseconds slowed down at room temperature (31). If these results are applied to our observation, we expect that the folding time may change rather sensitively with the GdnHCl concentration as it will affect the population ratio of microscopic unfolded states and their thermodynamic energetics. For example, at a higher concentration of GdnHCl the population of unfolded conformation with larger Rg values and extended structure will increase, which may bring about the increase of folding time correlated with the unfolded conformation of cyt-c.

We note that the chain diffusion timescale in cyt-c is several tens of microseconds (65). Considering this relatively fast backbone diffusion, the slow interconversion among many conformers in the U′ state may seem incompatible. This is resolved when they cannot be interconverted by chain diffusion alone. There are at least two factors that contribute to increasing the barriers between the unfolded conformers in the U′ state: 1) heme–residue misligation and 2) nonnative configurations of proline residues. In the case of misligated conformers, the interconversion toward a native-like ligation requires the displacement of the misligated residues, which is not likely to occur by the chain diffusion. In addition, we cannot rule out the existence of nonnative configurations of proline residues in the unfolded cyt-c. According to previous studies (42, 66), the timescale for transcis isomerization of proline is much larger than the chain diffusion time. Overall, we speculate that the origin of the heterogeneous kinetics from U′ to F is the heterogeneity of unfolded conformers in the U′ state due to misligating residue–heme interactions and the proline isomerization which cannot be interconverted only by fast chain diffusion. This aspect eventually results in a broad distribution of folding times.

We also note that we have utilized MD simulations for generating decoy structures for the EOM process based on selections with the genetic algorithm. Thus, how well the structural ensemble from the MD trajectories represents the folded and the unfolded state is not of concern, as long as the trajectories can cover the thermally accessible configuration space well within the limit of accessible simulation times. Of course, it would be ideal if the simulated ensemble could directly match the scattering data without any ensemble adjustment through EOM. In fact, when we compare the Rg distribution from any structure pool presented in Fig. 2 against the distribution after its EOM treatment for the unfolded state, we find that the structure pool data display higher probabilities on the smaller Rg side. Namely, the MD structures are generally more compressed than the actual experimental structures, which was also reported by the theoretical study of disordered proteins (67). Thus, we can infer that the adopted force field parameters unphysically prefer more compact structures even at the high temperature (800 K) employed for our simulations. It is likely that reproducing folded state characters were more importantly considered while optimizing the force field. Improving it such that the characters of less structured states are better reproduced will be important for studying proteins in general, especially with intrinsically disordered proteins. Our approach may provide a useful tool toward this end as it generates a straightforward target for benchmarking. Namely, by comparing the theoretical scattering data based on the simulated ensemble against the actually measured experimental figure, one can evaluate the fidelity of the adopted force field model for both folded and unfolded states. Although our EOM analysis scheme with MD simulations successfully extracted the structural information about the heterogeneous ensemble of unfolded cyt-c, there remains room for further improving on how to generate the candidate protein conformations. For example, the molecular form factor (MFF) analysis can be an alternative way for denatured and intrinsically disordered proteins as SAXS data can be analyzed by predetermined dimensionless MFFs without any additional MD simulations (68).

Conclusions

In summary, using TRXSS, we were able to establish the folding dynamics of cyt-c with the nonexponential kinetics in the transition between two states (U′ and F). The combination of the results obtained from the kinetic and structural analyses of X-ray scattering data revealed that the various protein conformations in the unfolded state take the multiple relaxation pathways toward the folded conformation, resulting in stretched exponential behaviors. Fig. 5 summarizes the detailed scheme for the spontaneous folding dynamics of cyt-c, especially highlighting the multiple relaxation pathways linked with the conformational heterogeneity in the unfolded state. This observation is direct evidence of protein folding from the heterogeneous unfolded states in view of microscopic protein conformations. The populations that can form nonnative ligation on the heme group are likely to induce various misfolded structures located in local minima in the folding free energy. Fast escaping from these minima will result in fast folding kinetics, while getting trapped in deep energy minima will require much longer time and detours toward the native state. The existence of all those complex pathways, associated with various states of heme ligations and the related conformers, will surely add a level of heterogeneity to its folding and unfolding dynamics. When such complicating factors combine together, they will contribute to rendering the kinetics to appear diverse with a stretched exponential behavior for cyt-c. Thus, how much our conclusion regarding the folding path diversity can be generalized should be considered as an open question, and answering that question in the future will definitely be interesting and important. We believe that our approach can surely contribute toward addressing this by further studying other folding systems. If the TRXSS experiment is implemented at various temperatures as employed in previous studies (14, 18), it may allow us to estimate the ruggedness of the folding free-energy landscape and to unravel the conformational heterogeneity in the early phase. This addition will ultimately facilitate topologically mapping the overall folding pathways in view of the microstate ensemble and its conformations. Such a capability will also be important for developing new computational models. In fact, constructing protein interaction models that can properly describe unstructured states is not a trivial task, primarily because experimental data that can validate the models are severely lacking. Unstructured protein states are more affected by interfacial behaviors such as solvent exposure and related polarization effects on protein surfaces, for which modern models are still being heavily developed. Our method can provide a convenient tool for evaluating their reliability. Thus, experimental developments in this area may contribute well toward advancing prediction powers of newly designed theoretical tools.

Fig. 5.

Schematics of the folding dynamics in a funnel-like free-energy landscape. Upon the photoreduction of cyt-c, the heterogeneously populated initial unfolded state (red dots in the projected plane) undergoes spontaneous folding along the multiple parallel pathways (white arrows). Based on the EOM analysis with one MD trajectory over 1-μs duration, four conformations were obtained to represent the unfolded state. The number of representative unfolded conformations somewhat depends on the trajectory identity. The rugged features around the top of the free-energy funnel may be associated with the existence of intermediates (yellow dots in the projected plane) that are visited in the early phase of folding. The observed stretched exponential folding kinetics is the result of the superposition of the multiple pathways from the unfolded substates (Ui, i = 1 to N) to the folded state (F) via intermediate substates (U′i, i = 1 to N), indicating the conformational heterogeneity in the cyt-c folding.

Methods

Sample Preparation.

Equine heart cyt-c purchased from Sigma Aldrich was dissolved in 100 mM of Na-phosphate 20 mM, 20 mM of NADH, and 3.5 M of GdnHCl at pH 7.0. A centrifuge with the cellulose acetate membrane filter was used to remove the aggregated particles in the protein solution. After that, the protein solution was bubbled using nitrogen gas for 30 min in order to make an oxygen-free buffer. The protein samples were prepared just before the X-ray scattering experiments.

Steady-State Spectroscopic Measurements.

Prior to starting the time-resolved experiment, we performed the steady-state spectroscopic experiments using UV-visible absorption and UV circular dichroism (UV-CD) spectroscopies (JASCO J-815) to check the state of the protein sample.

Measurement of X-Ray Scattering Curve.

Steady-state X-ray scattering data for the oxidized and the reduced form of cyt-c were collected at the 14-ID-B beamline at Advanced Photon Source. For the steady-state measurement for the reduced cyt-c, the oxidized cyt-c was reduced by the addition of sodium dithionite (4 mM) into the oxygen-free buffer. During the measurement, a 1-mM cyt-c solution dissolved in 100 mM Na-phosphate, 20 mM NADH, and 3.5 M GdnHCl at pH 7.0 was flowed to the capillary flow cell in order to avoid the signal from the damaged sample. The X-ray probe pulse at 12 keV, generated from the 324-bunch mode operation at APS, was delivered to the protein sample, and X-ray scattering patterns were recorded using a Mar165 CCD detector (Rayonix). The energy profile of X-ray pulse used for the collection of scattering pattern was measured by using a photodiode and was used for the structural analysis of X-ray scattering data in order to correct the polychromaticity in theoretical scattering curve (see , for details).

TRXSS data were collected at the ID09B beamline at the European Synchrotron Radiation Facility (ESRF) by following an established experimental protocol based on the pump-and-probe scheme (45–47, 50, 51). For the measurement, a 2-mM cyt-c solution dissolved in the buffer, which is the same as that used in the steady-state measurement, was flowed into the capillary system by using a syringe pump. To prevent reoxidation during the time-resolved measurement, we used an airtight circulation system where the bulk solution was fully degassed. Furthermore, to fundamentally remove the possibility that the folded form (reduced form) generated as photoproducts can undergo reoxidation, the sample portion exposed to the excitation laser and the X-ray pulses were flown through and discarded. Namely, the protein sample was exposed for laser and X-ray pulses only once. The cyt-c solution was irradiated by circularly polarized nanosecond (ns) laser pulses (Opotek Vibrant) with the fluence of 1 mJ/mm2 at 355 nm and was probed by X-ray pulses incident at well-defined time delays. Polychromatic X-ray probe pulse with the energy of 18 keV was generated from the 16-bunch mode operation at ESRF. Time-dependent X-ray scattering patterns were recorded using a Mar165 CCD detector (Rayonix), covered up to the wide-angle region. The scattering pattern at a negative time delay (−50 μs) was also collected and used as a reference for calculating time-resolved difference scattering curves. The scattering pattern at −50-μs time delay contains structural information of the initial state, while the data at positive time delays contain the contributions from a mixture of the initial state and photoproducts. The time-resolved difference scattering curves, ∆S(q, t) = S(q, t) – S(q, −50 μs), were generated by subtracting the scattering curve at −50 μs from the curves at positive time delays. The ∆S(q, t) curves provide the information on the change caused by the laser photoexcitation. We measured the TRXSS curves in a wide time range from 31.6 μs to 316 ms. At each time delay, more than 300 scattering curves were averaged to achieve a high signal-to-noise ratio.

Data Availability.

All data are included in the manuscript and .