Modeling time series with dynamic disorder

Markov chain approach is ubiquitously used in modeling biological systems. For example, reversible conformational transitions of biomolecules probed by single molecule fluorescence resonance energy transfer (smFRET) or force spectroscopy are often modeled as a homogeneous Markov process in which the transition rates between experimentally discernible conformational states are uniquely decided. To decipher time series with dynamic disorder that change their dynamic pattern from one time interval to another we assume that there are hidden “internal states”, each of which determines the rate of conformational transitions. A signature of the transition between internal states, which gives rise to dynamic disorder in time series, are difficult to detect using the value of FRET efficiency or end-to-end distance alone when the values observed along the time series are indiscernible even if the internal state is altered. By assuming that the transition between internal states is described by a homogeneous Markov process, and that transition between observable (in this study, FRET) follows non-homogeneous Markov process, whose transition rates are slaved to the internal state at each time, we model time trajectories made of these two layers of Markov chains. This algorithm corresponds to the Double Chain Markov Model (DCMM) [39–43] (Fig 2C and 2D).

DCMM is characterized by the following model parameters: (i) Transition matrix A for homogeneous Markov chain, which describes the transition probability between the K-distinct internal states along the time series (x = (x(1), x(2), …, x(t), ⋯, x(T − 1))). Here K is a total number of internal states in the model, and x(t), specifying internal state at time t, takes one of the values between 1 and K. T is the total observation time. The internal state at time t+1 (x(t + 1)) is determined by the previous internal state at time t (x(t)), whose transition to x(t + 1) is determined by a K × K Markov transition matrix A as where P(x(t) = ν) denotes the probability of x(t) being in the ν-th internal state; (ii) K transition matrices B^(μ) with μ ∈ {1, 2, …, K} for non-homogeneous Markov chain describes the transition probability between the observable states along the time series (o = (o(1), o(2), …, o(T))). o(t) specifies the state of the observable among N possible states {1, 2, …, N} at time t. Transition from o(t) to o(t + 1) is determined by an N × N transition matrix B^x(t)(t), the matrix elements of which are slave to the value of x(t)(= μ ∈ {1, 2, …, K}).

For example, if there are two (K = 2) internal states, and each internal state has three (N = 3) observables in a given time trace recorded with time resolution Δt, then two transition matrices for o with μ = 1, 2 can be considered (i.e., B⁽¹⁾ and B⁽²⁾): Next, the transition matrix A for the interconversion between two internal states is: In the above matrices, the matrix elements must satisfy, for each i = 1, 2, 3 in B^(μ), and γ^(1)→(1)Δt + γ^(1)→(2)Δt = γ^(2)→(1)Δt + γ^(2)→(2)Δt = 1 in A. More detailed descriptions about DCMM are available in the original papers [39–43] particularly in ref. [39] (see also SI). A similar but more general version of DCMM, which can accommodate inputs variables as well as multiple number of internal state sequences, has been suggested by extending the factorial hidden Markov model [46, 47].

Determining the number of internal states

DCMM can estimate the transition matrices A and B^(μ) quantitatively, and hence determine the most probable sequence of internal state and associated kinetic rates, and {γ^(μ)→(ν)}. However, the likelihood (the probability of observing data for given model parameters), maximized by DCMM, P(o|π, A, B) where π ≡ (π₁, π₂, …, π_K) with π_μ = P(x(1) = μ|o, A, B), is prone to increase when more number of parameters are used in the model. DCMM can select the best set of parameters for a given model, but not suited to select the best model (i.e., cannot determine the optimal number of internal states K for a given time trace). To overcome this limitation, often used is the maximum evidence method, where the evidence (P(o|K), also called marginal likelihood) is defined as the conditional probability of observing data (o) for a given model (K), so that (1) where λ ≡ (π, A, B) represents the parameter space. In this method, the penalty against model complexity is naturally incorporated during the calculation, allowing to select the best model (see SI). By calculating the evidence for each different model (different K, the number of internal states in data), one can select the best model with an optimal number of internal states that maximizes the evidence. The calculation of the evidence, however, involves a massive computational cost to explore the entire parameter space for a given model.

Variational Bayes double chain Markov model

To alleviate the computational cost in employing the maximum evidence method in Eq 1, we employ the Variational Bayes [48], a method that effectively uses a mean-field approximation. The method has previously been used to determine the number of observable states (FRET states) from smFRET data [49–51], the number of diffusive states from single molecule tracking data [52], and the number of DNA-protein conformations from tethered particle motion data [53]. It has also been used inside the empirical Bayes method which can analyze several smFRET time series simultaneously [54, 55]. In our study, the variational Bayes method combined with DCMM (VB-DCMM) was used to analyze single molecule time traces with dynamic disorder. The analytical expression of the lower bound of the evidence (F), offered by VB-DCMM, makes clear where the model penalty comes from, thus providing guidelines to choose the prior parameters to incorporate a prior knowledge of data (see SI). Once prior parameters are selected, VB-DCMM iteratively increases the lower bound of log(evidence)(= log P(o|K)) by identifying a better approximation to the true probability distribution. (2) where q(Z) is an arbitrary probability distribution of a set of variables, Z(≡ (x, λ)) consisting of parameters and hidden variables of model, and where D_KL(q||p) is the Kullback-Leibler divergence of q(Z) from P(Z|o, K), which we want to minimize. Once the solution from the algorithm converges, the approximate value of log P(o|K*)(≃ F[q*]) and the (locally) best model parameters (a set of the best kinetic rates), π*, A* and B*, which determines all the rate constants to describe the given time traces ( and {γ^(μ)→(ν)}), are acquired from an approximated probability distribution (See SI for the mathematical details). The performance of VB-DCMM is quite robust over a wide variation of prior parameters (S21 and S22 Figs).

Implementation of the algorithm

The observable sequence o is obtained by filtering the noise in the experimental data (o_n) using Hidden Markov Model (HMM) following a similar procedure as the previous studies [49, 56] using a custom code written based on the code from Sagemath software [57]. Next, the o is analyzed using VB-DCMM to select the best model and to estimate the best model parameters. The optimal sequence of internal states x is determined by using Viterbi algorithm [39]. All the implementations and data analysis are done by using our custom code. VB-DCMM is freely available at “https://github.com/TBiophysG/VBDCMM”

Validation of VB-DCMM

To first validate the efficacy of VB-DCMM in identifying internal states in a given SM time trace, we applied VB-DCMM algorithm on synthetic data that mimic a SM time trajectory with dynamic disorder (see Methods). To generate a synthetic SM time trajectory, we first produce a time trajectory specifying the value of internal state from t = 1 to t = T − 1. The time trajectory of internal state is represented with a symbol x ≡ (x(1), x(2), ⋯, x(t), ⋯, x(T − 1)). When the total number of distinct internal states in the model is K, one of the values in {1, 2, ⋯, K} is assigned to x(t). Thus, for K = 2 a typical time trajectory of internal state x looks like (1, 1, 1, ⋯, 1, 1, 2, 2, 2, ⋯, 2, 2, 1, 1, ⋯, 1, 1), (1, 1, 1, ⋯, 1, 2, 2, ⋯, 2, 2), (2, 2, ⋯, 2, 1, 1, ⋯, 1, 1, 2, 2, ⋯, 2, 2, 2), etc. The time trajectory given in Fig 3A(i) is an example generated with K = 2 and T = 8801. Next, similar to the structure of x, the time trajectory of noiseless observables is represented using o ≡ (o(1), o(2), ⋯, o(t), ⋯, o(T)). In Fig 3A(ii), a trajectory of o is shown, also demonstrating the influence of x on o. Finally, Gaussian noise was added on o and the range of signal was adjusted to produce the final trajectory o_n ≡ (o_n(1), o_n(2), ⋯, o_n(T)) which now resembles a time trajectory of SM FRET signal (Fig 3A(iii)).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Validation of VB-DCMM on synthetic data.

(A) (i) A time trace of internal state generated with γ^(1)→(2)Δt = γ^(2)→(1)Δt = 0.001. (ii) An observable time trace generated based on the trace of internal state in (i) by using internal state-dependent parameters , , . (iii) An synthetic FRET data with Gaussian noise overlaid on the trace in (ii). (iv) Noised filtered FRET state by HMM (blue line). (v) Traces of internal state with different K, estimated using VB-DCMM on the noise-filtered FRET trace from (iv) (black line is the true internal state trace while red, orange, and blue are internal state estimated from the model with K = 1, 2, and 3, respectively. The indices of internal state were determined by comparing B^(μ) estimated for each internal state with B^(μ),true which is used to generate the synthetic data). (B) Estimated lower bound of the evidence function F(K) of DCMM models with K = 1, 2, and 3. (C) Accuracy of detecting internal states. The overlap function χ calculated for 100 synthetic FRET traces generated under the identical condition used for generating the trace of internal state shown in (A).

https://doi.org/10.1371/journal.pcbi.1005286.g003

Deciphering the information of internal states from an observed time trace involves solving an inverse problem, i.e., decoding o_n to obtain x. To decode the trace of internal states from the synthetic data, we follow a 3-step procedure: (1) Filter the noise from o_n to obtain o using Hidden Markov Model (HMM) [56] (Fig 3A(iv), blue line); (2) Analyze o by applying VB-DCMM algorithm with different models 1, 2, …, K (again, K is the total number of internal states assumed in each model); (3) To select the best model we calculated the conditional probability of observing data for a given model parameter K, P(o|K), which is often called evidence or marginal likelihood in machine learning community (Eq 2) [48]. Calculation of P(o|K) is conducted using the Variational Bayes (VB) method, which gives the lower bound of log P(o|K) denoted by F(K). Details of the evidence function F(K) and approximation procedure are provided in the Supplementary Information (SI). Finally, we select the best model K* which maximizes F(K), i.e., K* = arg max F(K).

To be specific, in order to identify the best model parameter K for the time trace o(t) given in Fig 3A(iv), we varied K from 1 to 3. The most probable trace of internal states, , was calculated for each model with K = 1 (red), K = 2 (orange), K = 3 (blue) (see Fig 3A(v)). The evidence F(K) calculated using VB method was maximized at K = K* = 2, and the resulting time trace of the internal states, , most closely recovers the trajectory of x (black trace in Fig 3A(v)) except at the time interval where the transitions of x(t) between 1 and 2 occur only transiently or at the boundaries of transitions (red arrows on Fig 3A(v)). This result shows that VB-DCMM can avoid the over-fitting problem that other methods based on maximum likelihood are often fraught with [48].

Conditions required for an accurate recovery of internal states

VB-DCMM detects a signature of change in internal state (x) from a given observable time trace (o) by evaluating the statistical difference in transition rates. Thus, in the absence of an enough number of transitions in the trace o, the algorithm becomes less reliable. For example, we obtained F(K = 2) ≈ F(K = 3) although F(2) ≫ F(3) is more desirable (Fig 3B. See another example in S1 Fig). This is due to the lack of statistics in transition events in this particular test trace given in Fig 3A. For example, when only a part of the time trace is selected and analyzed using HMM, the estimated rates of transition from high (H) to low (L) FRET value are in 1500 ≲ t ≲ 4000, and in 5700 ≲ t ≲ 8700. Thus, in (K = 3)-model the two time intervals, originally generated by using the same kinetic parameter (), are determined to be distinct from each other (blue trace in Fig 3A(v)). By contrast, in (K = 2)-model, was estimated over these two time intervals. This type of statistical error is unavoidable for a small T_obs. A more systematic evaluation on the accuracy of the algorithm as a function of T_obs and transition rate between distinct internal states will be discussed in the next section.

To assess the accuracy of the best model predicted by VB-DCMM against the solution x, the following overlap function can be used. (3) where δ_i,j is the Kronecker delta and T = T_obs/Δt is the total number of data in the traces (Δt denotes the temporal resolution of the data). For 100 synthetic time traces, generated under the identical parameters used for producing the time trace in Fig 3A, we found that χ ≈ 0.9 on average (Fig 3C). Note, however, that x(t), only available for the case of “synthetic data”. Thus, to assess the accuracy of our method against a real time trace from SM experiments, we devised other metrics.

For a given time trace with dynamic disorder, our algorithm quantifies the kinetic features of the time trace in terms of the transition rate between the observable states a and b within the μ-th internal state and the transition rate from the μ-th internal state to ν-th internal state γ^(μ)→(ν) (1 ≤ μ, ν ≤ K, 1 ≤ a, b ≤ N. Here, μ is the index for internal state whereas a and b are indices for observable (In FRET displaying low/high two state transitions, these states correspond to the low and high FRET values). K is the total number of hidden internal states, and N denotes the total number of observables). To be able to extract the information of multiple internal states reliably from a time trace using VB-DCMM, two general conditions are required for the time trace being analyzed.

1. A large time scale separation should be present in the kinetics within each internal state, i.e., and (μ ≠ ν) should be disparate.
2. There should be a clear time scale separation between intra-basins and inter-basin transitions (i.e., τ_conf and τ_int). More precisely, the intra-basin transition probability should be much greater than the transition probability from the μ-th to any other internal state ∑_ν≠μ γ^(μ)→(ν) Δt (= 1 − γ^(μ)→(μ) Δt).

To substantiate the above-mentioned conditions 1 and 2, we define two metrics D_conf and D_int, which compute the average Hamming-like distances between the distinct rate constants extracted from a given time trace using VB-DCMM analysis: (4) and (5) D_conf measures the dissimilarity between distinct internal states in terms of the intra-basin transition rates. Two distinct internal states (μ, ν (μ ≠ ν)) can be better discerned if the intra-basin transition rate of one internal state (say, ) differs greatly from that of other internal state (), so that is maximized. D_int measures the average number of intra-basin transitions in each internal state using the ratio between the transition probabilities, and ∑_ν≠μ γ^(μ)→(ν)Δt (= 1 − γ^(μ)→(μ) Δt). A greater D_int ensures a large time scale separation in dynamics between intra-basin and inter-basin transitions, which improves the reliability of our method to decode the internal state from a given time trace. In general, D_int or D_conf shows a good correlation with 〈χ〉 (see below); thus, one can use (D_int, D_conf) to assess the accuracy of predicted internal states. Note that the metrics D_int and D_conf can be estimated for real data, while 〈χ〉 can be calculated only against the synthetic data. Since there is a good correlation between (D_int, D_conf) and χ, one can evaluate (D_int, D_conf), alternative to χ, to assess the reliability of a predicted result of .

To be more concrete, we applied VB-DCMM algorithm to analyze synthetic data generated with N = 2 (transitioning between high and low FRET values) and K = 2 (two internal states; μ = 1 and 2) under various scenarios.

• We fixed the transition rates in the state μ = 1 as , and varied the rates associated with the state μ = 2 over the range of (Fig 4A, left). For the interconversion probability between the two internal states we set γ^(1)→(2) Δt = γ^(2)→(1) Δt = 0.001. The accuracy of the model prediction (〈χ〉, Eq (3)) is on average greater than 0.9 as long as the transition rates and (μ ≠ ν) differ more than the factor of 4. Note that in Fig 4A (left), the value of 〈χ〉 is greater for , than for , ; this is because a statistically sufficient number of transitions make the detection of internal states more reliable. In contrast, when , , i.e. when the kinetics inside the two internal states are essentially identical, it is difficult to discern the two internal states. In this case, K = 1 instead of K = 2 is effectively the correct number of internal states. Indeed, when K = 1 is assumed (i.e., assuming true internal state x(t) = 1 for all t in Eq (3)), the re-calculated 〈χ〉 is close to 1 (see S2 Fig).
• To explore the effect of interconversion between distinct internal states on the performance of algorithm, we generated synthetic data with , , and by, this time, varying γ^(1)→(2) Δt and γ^(2)→(1) Δt = 0.00025 ∼ 0.005 (Fig 4B, left). The results clearly show that the case with smaller γ^(μ)→(ν) results in a higher 〈χ〉, which is expected because each internal state can have more number of transitions in the traces o when the interconversion is slower (Fig 4B). Re-plotting 〈χ〉 as a function of D_conf and D_int reveals clear dependence of the accuracy on D_int (Fig 4B, right). Similar trends are observed for other conditions of and (S3 Fig).
• Analyses on synthetic data generated using the same input parameters with those in Fig 4, but with a different number of data points in each trace, T_obs/Δt = 4400, and 2200 (S4 Fig) show a similar trend as observed in Fig 4 with T_obs/Δt = 8800 but with slightly smaller 〈χ〉 values.
• Extension of VB-DCMM algorithm to a more complicated case for K > 2 (S5 and S14 Figs) or N > 2 (S6 and S14 Figs) is straightforward. Application of VB-DCMM to a trajectory in which each internal state trajectory has different N is also straightforward (S7 Fig). In the latter case, the data is analyzed by assuming that all internal states have the same number of possible observables, N; but the analysis would indicate that transition associated with a small transition rate is essentially disallowed. In all situations considered for various K and N, VB-DCMM can be used for the reliable recovery of the sequence of true internal states.
• Analyses of synthetic traces show that the accuracy of the algorithm improves with both D_conf and D_int (Fig 4 right panels and Fig 5). Thus, these two metrics allow one to judge the reliability of the information on internal states extracted from a given time trace. Alternatively, a single parameter D_tot(= D_conf + αD_int) with an empirically acquired coefficient α ≈ 0.8) can be used to judge the reliability of the extracted information. Note that 〈χ〉 remains similar as long as D_tot remains constant (Fig 5B). Hence, when 〈χ〉 is plotted against D_tot, all synthetic data generated using different parameters approximately collapse onto a single universal curve (S8 Fig).
• There are multiple ways of assessing the efficacy of VB-DCMM in decoding the internal states. In addition to 〈χ〉, D_conf, D_int, and D_tot as the possible measures for the assessment, one can also use the statistical property that the dwell times of homogeneous Markov process satisfies (see SI for details).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Accuracy of VB-DCMM in detecting internal states under various conditions of and γ^(1)↔(2) with T_obs/Δt = 8800.

(A) The color bar denotes the accuracy of analysis in terms of 〈χ〉 under varying with K = 2, , , and γ^(1)→(2)Δt = γ^(2)→(1)Δt = 0.001. (B) 〈χ〉 under varying γ^(1)→(2) and γ^(2)→(1) with K = 2, . 〈χ〉 was calculated by averaging over the results from analysis of 100 traces in each condition. The panels on the right show the relation between the value of 〈χ〉 and pairs of D_int and D_conf values which are evaluated at varying kinetic parameters. Results from the analysis over the data with the same parameters but different length of time trace (or different number of data points T_obs/Δt = 2200, 4400) are provided in S4 Fig.

https://doi.org/10.1371/journal.pcbi.1005286.g004

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Average accuracy of internal state detection as a function of D_conf, and D_int.

To construct this diagram, we employed various synthetic data in Fig 4 (circle, two internal states (K = 2), two FRET states (N = 2)), S5 Fig (left triangle, K = 3, N = 2), and S6 Fig (hexagon, K = 2, N = 3). The right triangle symbol denotes the result from the similar analysis shown in S5 Fig with K = 3 but with smaller relative differences in the transition rates, k’s. Pentagon represents the result obtained with K = 4 and N = 2. (A) Color code denotes the accuracy of internal states predictions in terms of 〈χ〉, averaged over 100 traces for each condition. (B) The dashed lines corresponding to Δ = D_conf + 0.8D_int = 4, 5, … 9 are overlaid on the 2-D scatter plot of 〈χ〉(D_conf, D_int) calculated in (A).

https://doi.org/10.1371/journal.pcbi.1005286.g005

Application of VB-DCMM on H-DNA data

Now, to analyze the duplex-triplex transitions of H-DNA (Fig 6), we obtain o by filtering the noise from FRET signal (Fig 6(ii), blue line) and apply the VB-DCMM algorithm to decode the hidden internal state in the signals. Fig 6(iii) shows time series of internal state, , calculated from the VB-DCMM by varying K from 1 to 5. It is of note that the number of actually observed internal states in the for a given input parameter K does not change after some K_obs(≤ K) (K_obs = 2 (Fig 6A), 2 (Fig 6B), 2 (Fig 6C), and 1 (Fig 6D)). (See also other time traces of synthetic data and H-DNA analyzed in SI: S1A Fig (K_obs = 1), S9A Fig (K_obs = 3), S9B Fig (K_obs = 2), S9C Fig (K_obs = 3), S9D Fig (K_obs = 3), S10A Fig (K_obs = 4), S10B Fig (K_obs = 2), S10C Fig (K_obs = 2), S10D Fig (K_obs = 3), and S12 Fig (K_obs = 3)). A similar behavior is also observed when analyzing data using the variational Bayes Gaussian mixture model [48].

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Representative time traces of H-DNA at [Na⁺] = 100 mM and their analysis.

(A) (i) Fluorescence signal and (ii) their FRET state. (iii) Internal states estimated for K = 1, 2, …, 5. Right panel shows G(K) (blue circle) where K_obs specifies the number of detected internal states in individual traces (blue). (B, C, D) Other representative time traces and their G(K) obtained under the same experimental condition.

https://doi.org/10.1371/journal.pcbi.1005286.g006

To account for the contribution due to degeneracy in labeling the internal states, log K! term is conventionally considered in formulating the evidence function F(K) (See SI for the details); however, in our problems, the actual number of degeneracy in labeling internal states should be instead of K!. Therefore, we replace the log K! term in F(K) with log [K!/(K − K_obs)!], and considered a modified evidence function, G(K), to identify an optimal K for a given time trace: (6) G(K) shows a clear peak, allowing us to identify the optimal K(= K*) with ease (blue circles on the right side of Figs 6, S9 and S10). Use of G(K) instead of F(K) in analyzing synthetic data does not alter K* (S12 and S14B Figs).

Among the time traces of H-DNA, traces with more than 3 interconversions between distinct internal states, which enables us to estimate γ^(μ)→(ν), are rare, especially when [Na⁺] = 100 mM; thus it is not feasible to get a statistically meaningful scatter plot of (D_conf, D_int) (see S23D, S23E and S23F Fig); however, for those displayed in S23D, S23E and S23F Fig, 〈D_tot〉 ≈ 7 suggests that χ ≳ 0.9 (from Fig 5). Therefore, at least the intra-basin rate constants extracted from H-DNA data using VB-DCMM are reliable.

Time traces that have τ_int comparable to experimental observation time () would exhibit on average no or only a single transition event between distinct internal states. Indeed, we find that only a subset of total number of internal states is sampled by individual time traces due to the limited observation time. For instance, at [NaCl] = 100 mM, our analysis identified K* ≤ 2 in 265 out of 269 traces, and that only 4 time traces display K* > 2 (S11 Fig). Therefore, in order to identify the internal states present in the transition dynamics of H-DNA, clustering analysis is required against the whole ensemble of time trajectories. We provide the procedure of clustering analysis and results in details in the following section.

Clustering H-DNA data

VB-DCMM algorithm allows us to decompose individual H-DNA time traces with dynamic disorder into multiple “components”, each of which should satisfies the property of homogeneous Markov chain. In order to understand the structure of conformational space of H-DNA, the ensemble of components acquired from the VB-DCMM analysis should be clustered into the same kind. To this end, we produce scatter plots of (k_L→H, k_H→L), representing the kinetic property of the ensemble of time traces, using the transition rates estimated for individual time traces. The scatter plots of (k_L→H, k_H→L) were calculated for the ensemble of H-DNA time traces (i) before (Fig 7A, left) and (ii) after decomposing the individual heterogeneous time traces retaining multiple components into the homogeneous ones (Fig 7A, right). The scatter plot of (k_L→H, k_H→L) after the decomposition has a greater dispersion, which is expected since a data point (k_L→H, k_H→L) for a time trace with dynamic disorder is a mixture of with μ = 1, 2, … K. In the presence of clear distinction between internal states (μ ≠ ν), the clustering of would be straightforward, which is indeed the case for the synthetic data (S13A Fig). However, for the H-DNA data, even after the decomposition, the clustering of data on (k_L→H, k_H→L) plane (Fig 7A) is not that clear.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Clustering H-DNA data at [Na⁺] = 100 mM condition.

(A) The scatter plots of (k_L→H, k_H→L) before (left) and after (right) applying VB-DCMM from [Na⁺] = 100 mM data. (B) The “average pairing distance” as a function of the number of clusters ()(see Methods). The minimum value of is found at . (C) Left: the scatter plot of clustered data projected on plane. Right: scatter plot of clustered data projected on plane. The data belonging to different clusters are depicted in different colors, and the centroid of each cluster is marked with the × symbol. Total 98 data points were used for analysis. (D) The result of the above clustering is represented using 6 kinetic arrows which represent the centroids of each cluster represented in (k_L→H, k_H→L) plane. The starting point of the arrow is (, ) whereas the ending point of the arrow is (, ), where the superscript μ represents the index of each cluster μ = 1, 2, … 6. The colors used for depicting kinetic arrows are consistent with the data points in (C). (E) A schematic of the conformational landscape of H-DNA.

https://doi.org/10.1371/journal.pcbi.1005286.g007

To improve the quality of clustering, we extended the clustering of the kinetic data to a higher dimension by considering the kinetic information of internal states that are contiguous (kinetically linked) along time traces. To be specific, for a time trace exhibiting a transition from the μ-th to ν-th internal state (μ ≠ ν), one can consider that the inter-basin transition has occurred from the time interval represented by its pair of kinetic rate to the next time interval represented by , where the superscripts, ‘bf’ and ‘af’ denotes ‘before’ and ‘after’ the transition, respectively. Thus, instead of (k_L→H, k_H→L), a clustering at a higher dimension can be carried out by measuring the Euclidean distance between a pair of the four-dimensional (4-dim) arrays, .

In order to cluster the 4-dim arrays we used the k-means clustering algorithm. Application of the algorithm to the H-DNA data at [Na⁺] = 100 mM reveals that the average pairing distance, (see Methods), is minimized when the number of clusters is 6 (), namely, the model with 6 clusters provides the best interpretation of the data (Fig 7B). Although the model with 14 clusters shows a smaller , we selected as the best solution, since for each of 12 clusters out of 14 has less than 10 data points, which makes the result of clustering statistically less significant (S16 Fig). This results remain qualitatively identical when L1 distance (so called “city block” distance) was used instead of “square-euclidean” distance (S19 Fig). Furthermore, the clustering algorithm using “affinity propagation” [58], which considers all the data points as possible exemplars (analogous to centroids in k-means clustering method) and iteratively exchanges messages between them, also gives qualitatively identical results, confirming the robustness of the conclusion on H-DNA dynamics obtained from VB-DCMM and k-mean clustering (see S20 Fig).

We present the result of clustering either (i) by projecting it on the two separate kinetic planes, and , which visualize the inter-basin transitions in terms of variable L → H and H → L transition rates (see Fig 7C), or (ii) by using “interconversion arrows” linking the kinetic rates of two internal states, on the (k_L→H, k_H→L) plane (Fig 7D). Note that in the scatter plot visualized with , the distinction between different clusters is clear (the right panel of Fig 7C). Furthermore, for a system in equilibrium or at least near equilibrium, the interconversion between two internal states, say μ and ν, should occur in both directions, i.e., μ → ν and ν → μ. In the representation (i), a symmetry of is expected in the both panels of Fig 7C; and in the representation (ii), the “arrows”, amounting to the kinetic connectivity between distinct internal states, should be bi-directional. The symmetry of the data plotted in Fig 7C or the bidirectionality of the kinetic arrows confirms the condition of detailed balance being satisfied in the system in equilibrium. Fig 7D depicts 6 kinetic arrows (3 pairs of reversible kinetic arrows) connecting the centroids of or data.

Application of the above clustering method to synthetic data with K = 3, N = 2 (S13 Fig) is straightforward. To check the efficacy of clustering method for a more complicated case, we have tested with synthetic data generated with K = 4, N = 4, i.e. when there are as many as 4 observable states in each internal state (S14 and S15 Figs). In the case with 4 observable states, total 12 possible intrabasin transitions are conceivable. Thus, the dimension of the array associated with interbasin transition is 24. As long as there is a clear time scale separation, it is expected that the pairing distance shows minimum as there are 12 connection paths between 4 internal states. Indeed, is minimized at (S15A Fig).

Lastly, it is noteworthy that the clustering method presented here is not limited to data analysis for systems in equilibrium, but can be extended to systems in nonequilibrium steady state [59] where the individual state-to-state kinetic transition rate is well defined using the reversible Markov process although the condition of detailed balance is no longer anticipated [60, 61]. The symmetry of data point and bidirectionality of kinetic arrows as in Fig 7C and 7D are still of use to cluster the kinetic information generated from a system in nonequilibrium steady states.

Folding energy landscape of H-DNA

We classified the “components” of a similar kinetic pattern (k_{L → H}, k_{H → L}) obtained from VB-DCMM into a single cluster which represents a kinetic path linking two independent basins of attraction (or internal states). For example, the kinetic paths in Fig 7D can be best understood by hypothesizing 4 internal states (four basins) linked by 6 kinetic paths. Thus, the conformational transition landscape of H-DNA at [Na⁺] = 100 mM condition consists of 4 internal states with 3 reversible kinetic paths being established as illustrated in Fig 7E. At lower salt concentrations ([Na⁺] = 50 mM (S17 Fig) and [Na⁺] = 26 mM (S18 Fig)), H-DNA transitions slow down and the dispersion of data also increases; however, the overall structure of conformational landscape of H-DNA remains unchanged from the picture suggested in Fig 7E; thus, there is a central superbasin to which three other superbasins are kinetically connected (S17 and S18 Figs).

Contributions of our work

In comparison to other pre-existing methods, the advantage of our VB-DCMM in decoding dynamic disorder from a given trajectory is highlighted as follows:

(1) Dynamic disorders in single molecule time trajectories are modeled using DCMM by assuming the presence of hidden internal states. While Aggregated Markov Model (AMM), which has been adopted in ion-channel community for time trace analysis of varying current [62–76], can be employed to analyze our data with dynamic disorder, DCMM is better in correctly decoding dynamic disorder than AMM. We found that AMM is prone to overpredict the transition between kinetic patterns (S25 Fig). Our method is more suitable to the data showing persistent dynamic patterns by suppressing unwanted frequent transition between kinetic patterns. Detailed explanations of connection and quantitative comparison between DCMM and AMM are provided in SI and S25 Fig.

(2) In this paper, Bayesian version of DCMM was developed by using variational Bayes (VB) method, which enabled us to determine the number of internal states straightforwardly. Although Bayesian version of DCMM using Markov chain Monte Carlo (MCMC) method has previously been developed for the credit portfolio modeling [77], the idea of Bayesian inference in ref. [77] was used only for the purpose of calculating a posterior distribution of model parameters. To determine the number of hidden states corresponding to the internal states in this study, the authors in ref. [77] used the economic cycle fluctuation model, instead. Our study combining VB with DCMM (i) can determine the number of internal states in a more objective fashion, (ii) offers intuitive way to incorporate prior knowledge, and (iii) is computationally more efficient than MCMC (See SI for details).

(3) We tested VB-DCMM under various conditions, by varying the kinetic rates, the number of observables, the number of hidden states, and prior parameters. New metrics were also devised to quantify the performance of algorithm systematically.

(4) Finally the connection paths (kinetic arrows) between internal states of H-DNA are clustered by using the kinetic components extracted from VB-DCMM and by applying k-means clustering algorithm to high dimensional arrays.

To recapitulate, our entire process of analyzing single molecule data is composed of three stages: (i) noise-filtering using HMM; (ii) decomposition of heterogeneous time traces into the homogeneous components using VB-DCMM; (iii) clustering the decomposed components into the same cluster.

In principle, this three-stage analysis can be made more systematic by combining the noise-filtering and clustering procedure with VB-DCMM. To be more specific, (1) The noise-filtering of observable trace (o_n) is processed, independently from the main VB-DCMM algorithm, by using HMM, which has been proved to be reliable in noise-filtering [56], and the maximum number of observables (N) are predetermined as an input parameter. Current version of algorithm can be further automated by combining with the Bayesian version of HMM [49], which can determine the number of observables while filtering the noise in data (See Fig 2D). The resulting model will have a similar structure with the modified factorial HMM [46, 47]. (2) The heterogeneous components identified from individual time traces are clustered separately from our main algorithm. It would be also desirable to unify the post-processing step (clustering) with VB-DCMM using empirical Bayes method which has been applied recently to analyze single molecule data [54, 55].

However, it should also be noted that a blind integration of noise-filtering and clustering steps inevitably complicates the implementation of VB-DCMM, as more number of prior parameters are ought to be decided by users. For example, Bayesian implementation of HMM for noise filtering demands manual determination of additional N(N + 5) prior-parameters [49]. Compared to this, currently VB-DCMM requires users to pre-determine only one prior parameter which characterizes the final transition rate matrix, A (see the subsection: Selection of prior parameters in SI). Moreover, the integration of other methods will obscure the flow of analysis, making it difficult to identify an error-causing step. Keeping each step in the algorithm separate makes the integration of VB-DCMM to other applications more transparent (for example, if noise-filtering by HMM is unsuccessful, other advanced method can be employed [49]). We leave it as our future work to develop an algorithm that integrates the above-mentioned three procedures (noise-filtering, VB-DCMM, and clustering) without increasing complexity or obscuring the flow of analysis.

In decoding SM FRET data, the most notable difference of our VB-DCMM from the previous studies employing the probabilistic models such as maximum likelihood and Bayesian statistics is that VB-DCMM explicitly considers the situation that transition rates can change from one time interval to another within individual time traces. The previous studies [49–51, 56, 78–80] assumed that the transition rates were constant within individual time traces. Also, currently, VB-DCMM is applicable to window-averaging FRET trajectories. It will be of great interest to extend VB-DCMM to analyzing time trajectories in which arrival times for individual photons are available. VB-DCMM is particularly powerful when there is a separation in time scales between τ_int and τ_conf.

Concluding remarks

While the notion of dynamical heterogeneity or broken ergodicity seems better recognized in the research field of nucleic acids [81] than in proteins, which likely arises from more homopolymer-like nature of building block of nucleotides [82], biomolecules in general can have a rugged folding landscape with many local basins of attraction and kinetic barriers with varying heights [83]. Conformational dynamics of biomolecules on rugged landscapes can be heterogeneous, which gives rise to static or dynamic disorder depending on the time scale of observation or the height distribution of kinetic barrier. The presence of heterogeneity or disorder among individual molecules, unveiled by in vitro SM experiments could be surprising at first sight; however, it is also important to note that the general hypotheses in the conventional molecular biology towards a single native state have been put forward based on the observations from ensemble experiments where the heterogeneity, if any, is usually masked by the process of ensemble averaging. Given that the complexity of a molecular system increases with the system size (N_sys) as [29], it should not be too surprising to find such disorder in biomolecules in itself. Cells are equipped with molecular chaperones that can tame misfolding-prone biomolecules with rugged landscapes [84–87]; thus the principle of optimization in biology, if it fails at the level of a molecule in isolation, can be extended further to the molecular system including its environmental factors.

It is not easy to elucidate the molecular origin of disorder in a conclusive manner; yet, it has recently been suspected that interactions of biomolecules with cofactor such as ATP and multivalent metal-ions could be the microscopic causes for those molecules exhibiting dynamical heterogeneity [12, 13, 16, 88, 89]. Modulating the concentration of Mg²⁺ ions from high to low and again to high induced inter-conversions of dynamic patterns in equilibrium conformational fluctuations of T. ribozyme [12] and Holliday junctions [13]. Distinct velocities of ATP-empowered individual RecBED helicase motors, which can move progressively along dsDNA by unwinding it into two separate strands, can be reset by introducing a long pause by halting the supply of ATP. For the time trajectories of biomolecules displaying quenched disorder, a method to analyze such data was proposed using a concept from glass physics [13]. Here, to deal with more general scenarios, we have developed a method to analyze single molecule time traces with dynamic disorder.

As demonstrated by testing the VB-DCMM algorithm on synthetic data, the algorithm is quite accurate in decoding dynamic disorder as long as a time trajectory of interest contains multiple time intervals, each of which display kinetic pattern distinct from others. When a clear separation in timescale is present between two distinct kinetic patterns, large value of D_conf, D_int, and D_tot would be acquired.

While we developed the VB-DCMM algorithm primarily to analyze dynamic disorder in duplex-triplex transitions of H-DNA, the method is applicable to any data in the form of one-dimensional time series with multiple transitions. Together with a further technical advance in SM, which eliminates experimental artifacts as well as extends the measurement time, our algorithm developed here will contribute to better understanding of biomolecules that display heterogeneous dynamics.

Generation of synthetic data

Internal state sequence x was generated by using Monte Carlo method with a constant transition matrix (homogeneous Markov chain model). The observable sequence o was generated by using the same method but with the transition matrix that was defined at each time t based on the internal state x(t). Finally, Gaussian noise was added on o to produce o_n.

Single-molecule FRET measurements to monitor duplex-triplex transitions of H-DNA

We purchased triplex forming oligonucleotides from Integrated DNA Technologies (Coralville, IA, USA). The oligonucleotides were dissolved in T50 buffer solution (10 mM Tris-HCl, 50 mM NaCl, pH = 7.5) and were heated beyond the melting temperature of DNA duplex (∼ 90°C), and slowly cooled down on a heat block to room temperature over 8 hour to properly hybridize them. The DNA prepared as such is called “H-DNA” here. The sequences of the triplex forming strands (purine-rich and pyrimidine-rich) are: Purine-rich strand: 5’ AAG AAG AAG AAG AAG (Cy5) TGG CGA CGG CAG CGA (Biotin) 3’, Pyrimidine-rich strand: 5’ TCG CTG CCG TCG CCA CTT CTT CTT CTT CTT TTT TCT TCT TCT TCT TCT TC (Cy3) 3’. In the purine-rich strand, the biotin at 3’ terminus is used to attach the H-DNA molecule to a neutravidin-coated cover-glass. The Cy3 and Cy5 dyes in the H-DNA molecule correspond to a donor and an acceptor for FRET measurements, respectively. In order to observe the transition between folded triplex and unfolded DNA, we used the reaction buffer containing 50 mM HEPES(Sigma-Aldrich) and various concentrations of Na⁺ (26, 50, 100 mM). These buffer solutions also contained 2 mM trolox, 10% glucose and gloxy for single-molecule fluorescence experiments. We utilized a home-made TIRF (Total Internal Reflection Fluorescence) microscope to measure the FRET efficiency between donor and acceptor dyes, which reveals the conformational state of the H-DNA molecule. A 532-nm laser (CrystaLaser DPSS, 10 mW) was used to excite donor molecules and fluorescence intensities of both dyes were measured by an EMCCD (Andor iXon DV887, Andor technology). To observe the change of FRET efficiency in real time, we measured the time-lapse FRET traces with the repetition rate of 10 Hz. To study kinetic features of the conformational transition with dynamic disorder, we acquired the FRET time traces for a long period (>100 sec).

Clustering at a higher dimension

For given N and K, total N(N − 1) intra-basin transition rates (a, b ∈ {1, 2, ⋯, N}, a ≠ b) are defined in the μ-th basin (or μ-th internal state) and total K(K − 1) inter-basin transitions are conceivable. To cluster the kinetic information of H-DNA data obtained from VB-DCMM, we consider the kinetic arrow, 2N(N − 1)-dimensional array of data, which has the structure of where the subscript i denotes an index referring to one of K(K − 1) possible inter-basin transitions linking two internal states (μ ≠ ν). For a kinetic scheme made of a network of reversible transitions between K internal states, the transition between two internal states should be bidirectional; thus for a given inter-basin transition path i, there should be a kinetic path j antiparallel to the path i, satisfying , where . In our problem, the set of all the data generated as an outcome of VB-DCMM can in principle be clustered into the disjoint subsets of size 2 partitioning the transition paths, , and one realization of such disjoint subsets will minimize the pairwise sum of Euclidean distances for all α and β; however, the method suffers from high computational cost as the possible number of clusters increases rapidly with N and K.

To alleviate the computational cost for large N and K, we modified the original method. We first searched the the best partitioning set of data for a given that minimizes the Euclidean distance between all the pairs of centroids, (7) where with , , and 〈…〉 denotes the centroid of clustered data. To obtain the best clustering result for a given , we conducted k-means clustering using k_means function from scikit-learn libraries [90] with 20,000 different random initial conditions in each analysis. It is expected that . The summation, ∑_(i,j), signifies that the sum is taken over the disjoint subsets of size 2 partitioning a set with being an even number) and is the best partitioning set that minimizes the value of for a given . For example, provided that there are 4 kinetic arrows made of centroids (i = 1, 2, 3, 4), which minimizes at when i = 1 is paired with i = 3 and i = 2 with i = 4, then S*(2) = {{1, 3}, {2, 4}} and .

Next, in order to decide the optimal , we calculated pairing distance between paired clusters in again, but this time using all the elements in each cluster. The total pairing score (8) where the average pairing distance between two clusters i and j is defined as where , , and n refers to an index for the element in the i-th cluster and m to an index for the elements in the j-th cluster. M_i is the total number of the elements in the i-th cluster. Finally, the optimal , minimizing , is selected, i.e., , and the interpretation of data is conducted for the best partitioning set .

For H-DNA data at three different Na⁺ concentrations, the optimal are determined at for [Na⁺] = 100 mM (Fig 7A), for [Na⁺] = 50 mM (S17B Fig), for [Na⁺] = 26 mM (S18B Fig). This implies that the complexity of conformational space of H-DNA increases at low salt condition (also see the scatter plot of (k_{L → H}, k_{H → L}) in Figs 7A, S17A and S18A).

The clustering results presented in this study remain robust regardless of the choice of distance metric. K-means clustering using L1 distance (“city block”) measure with 20,000 different random initial conditions also was led to qualitatively similar results (S19 Fig). Furthermore, as an alternative clustering algorithm, we also tested “affinity propagation” [58] on our data, and the results remain qualitatively identical (see S20 Fig). In the affinity propagation method, negative square-euclidean distance was employed as a similarity metric (s(i, j) = −||x_i − x_j||²) where x_i denotes the coordinate of the i-th data point. The objective of the algorithm is to optimize the factorized probability distribution which approximates the net similarity , defined as . Here, c_i is the index of the exemplar of i-th data point x_i. For example, if c_i = k, x_k is an exemplar of x_i and x_i belongs to the cluster represented by x_k. Multiple iterations of message passing are carried out until convergence is achieved in the result and the best result of clustering is acquired. For implementation, we used AffinityPropagation class from scikit [90] library with varying “preference” as an input parameter, where the preference denotes the logarithm of probability that i-th data point x_i selects itself as an exemplar. Further details of the algorithm are available in Ref. [58].