| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
Department of Physics and Howard Hughes Medical Institute, University of Illinois at Urbana-Champaign, Urbana, Illinois
Correspondence: Address reprint requests to Taekjip Ha, E-mail: tjha{at}uiuc.edu.
| ABSTRACT |
|---|
|
|
|---|
| INTRODUCTION |
|---|
|
|
|---|
To accomplish this, historically one needed to determine, based on the FRET values, which state the system was in at a given time and for how long it remained there. Typically, histograms were built out of the resulting dwell times and then fit to an exponential decay (1
7
). The method of deciding when the molecule is in what state can vary from simplistic "by eye" analysis to the slightly more sophisticated thresholding algorithm. In the "by eye" case (6
), one determines, based on one's own experience, what changes in FRET are legitimate state-to-state transitions as opposed to noise or photophysical effects. In this approach, shorter transitions are likely to be missed and the results may vary from person to person. The other common approach is the thresholding algorithm (8
) that requires a user to set cutoffs between FRET states. Still needed is a way of distinguishing a noise-induced change in FRET from a legitimate short-lived transition, and one needs to account for the fact that different molecules, even from an entirely homogenous population, can show different state FRET values due to instrumental noise. In both cases, analysis is performed with a preconceived number of states in mind and is subject to the user's prejudice for systems of great complexity. Likewise in both cases the ability to reliably and reproducibly analyze data becomes almost impossible when the number of states being analyzed increases beyond three. To overcome the shortcomings of these approaches we have devised an algorithm based on hidden Markov modeling that seeks to determine in a probabilistic and less user-dependent way the number of conformational states of a system and the rates of exchange between them.
| BASIS OF THE ALGORITHM |
|---|
|
|
|---|
A Markov process consists of any combination of state-to-state transitions with kinetics governed by single-exponential decay. A Markov process becomes hidden when an observer's ability to detect states is obscured by noise. Let us consider a system that can adopt one of two conformational states (A or B), where the probability of observing state A transit to state B in one time step is governed by single exponential decay and independent of how long the molecule has been in state A. If a FRET pair has been placed on the molecule such that the ideal, noiseless FRET values obtained from each state are FRETA and FRETB (Fig. 1 A), in principle, one could obtain a molecule's true A
B trajectory by following the FRET value in time as it alternates between FRETA and FRETB (Fig. 1 B). In reality, experimental noise broadens the FRET distribution of each state (Fig. 1 C) yielding significantly more complicated data (Fig. 1 D). This makes the process a hidden Markov process as the underlying reality (the sequence of true A
B exchanges) is hidden in the data because of noise. To reconstruct the underlying reality using HMM, first we need to define model parameters.
|
(either A or B) transiting to state
(either A or B) in the next time step, in general called tp

. Of the total of four transition probabilities (tpA
B, tpA
A, tpB
A, and tpB
B) forming a transition probability matrix (Fig. 2 A), only two are independent since by definition tpA
A + tpA
B = 1 and tpB
B + tpB
A = 1. It is important to note that in using a transition probability matrix of this form one is implicitly assuming that the underlying process is a Markov process. This means that transitions are assumed to be governed by single exponential decay kinetics.
|
(with idealized FRET value FRET
). We have modeled the noise-induced FRET distribution using a Gaussian function with a characteristic width
. Then, the two emission probability functions (epA and epB) are given as follows:
![]() | (1) |
Deciding between different event sequences
Once we have a transition probability matrix and emission probability functions for a system we can evaluate the relative likelihood of any proposed sequence of A and B states given the observed data FRETdata(i) (where i represents time step). We term the time-indexed proposed sequence of states
so that
(i) = A or B. To evaluate how likely a proposed sequence
is we calculate, for each time point, the probability that the proposed state yields our observed data (using the emission probability function), and then determine the point's transition probability to the next time step's proposed state. Multiplying the two terms we obtain the probability of the proposed state yielding the observed data at that time point:
![]() | (2) |
![]() | (3) |
As a specific example, we provide the seven-time point-trajectory in Fig. 2. The system follows the parameters in Fig. 1 but with a narrower width (
= 0.16 instead of
= 0.35). Fig. 2, C and D, attempts to fit the same set of data with two different proposed trajectories. In Fig. 2 C, the proposed trajectory
1 has a transition from B to A in time step 3 and the reverse in time step 4. In the trajectory
2 proposed in Fig. 2 D no transition takes place. To determine which of the two is the most likely given the model parameters in Fig. 2, A and B, we use Eqs. 2 and 3. Although it is much more likely for a molecule in state B to remain in state B as in the trajectory
2 than to transit from state B to state A and then back again as in
1, it is even more likely for a molecule in state A to emit a FRET value at 0.39 as in
1 than a molecule in state B to as in
2. Thus, the more likely trajectory among the two is
1.
To try every possible sequence of states and find the total probability for a FRET trajectory of N time points would require O(2N) computational power for a two-state system. Fortunately this can be cut down to O(N) using the Viterbi algorithm (19
), which is guaranteed to find the most likely sequence of underlying states given a set of data and corresponding transition probability matrix and emission probability functions.
Determining emission probability functions and transition probability matrices
Thus far we have assumed that the emission probability functions and transition probability matrix are known beforehand. But in an experiment, one does not usually know the idealized FRETA, FRETB values or their distributions. In addition, the transition probabilities are what one is after in the first place. As the Viterbi algorithm also provides the total probability of the most likely sequence of states, we can vary the model parameters (FRETA, FRETB,
, tpA
B and tpB
A) until we achieve a maximized total probability. This is a well known multi-dimensional optimization problem which we solved using Brent's algorithm (20
). We have found that as long as reasonable first guesses are made with respect to the parameters the algorithm converges to the true values rather than to a local maxima. As initial guesses, we use uniform transition probabilities of 0.005 and uniformly distributed FRET states between 0 and 1, or rough estimates based on the data. Similar results could in principle be obtained with other HMM algorithms such as the faster but more complicated Baum-Welch forward-backward algorithm (21
). All algorithms were encoded in C++ and analysis performed on the UIUC Turing CPU cluster.
Utilizing multiple trajectories
In the hopes of obtaining a more accurate tpA
B and tpB
A, one would generally like to determine them for a number of different molecules and then find some way to average them. We have found that transition probabilities are distributed asymmetrically as in Fig. 3 A. Therefore to obtain a representative average value and corresponding uncertainty we first take the transition probabilities returned from the HMM algorithm and find their logarithm (as in Fig. 3 B). The logarithm is chosen partially for because it yields symmetry, but also because of the relationship between free energy and kinetic rates:
. That is, if the free barrier has heterogeneous broadening that is symmetric, the distribution of the logarithm of the transition rate would be symmetric. Once the logarithm has been taken, a mean is determined together with its standard error. These are then converted back to transition probabilities by exponentiation. Explicitly, for a collection of N transition probabilities labeled by i:
![]() | (4) |
![]() | (5) |
B and
TPA
B denote the representative mean and the representative error, respectively.
|
. We chose a single value of
for all states to reduce computational time. It is often the case that the value of S itself is unknown. In such cases the solution is simple: assume some large number of states (say M = 10) and perform analysis. If there are in reality only five states (S = 5), five legitimate (between 0 and 1) FRET values and their associated transition probabilities will be returned together with five extraneous states. The extraneous states are never populated.
Another potential difficulty in systems with a large number of states is one of how to categorize FRET states. Let us consider an example of a five-state system (states labeled A, B, C, D, and E). In some cases individual FRET trajectories may not show all five states (for instance, if the trace is not long enough). Or the FRET values found for each individual state may not agree completely between trajectories (one molecule's state A may have an apparent FRET = 0.25, whereas for another molecule it may have an apparent FRET = 0.2 due to experimental variability). Occasionally one might even see an additional state, C', which was found by the algorithm but differs only slightly from one of the existing states C.
To visualize the number of actual FRET states (S) in a way that avoids these pitfalls, we developed a simple two-dimensional pseudo-histogram we call the transition density plot (TDP). For each FRET trajectory analyzed (where the algorithm attempted to fit M different states), the algorithm finds M idealized FRETm levels (m = 1, 2, ..., M), the transition probability matrix, and the number of transitions found for each of the FRETm
FRETm* pairs of FRET levels. A two-dimensional Gaussian function is constructed for each pair of start and stop FRET values (startm, stopm), with an amplitude (am,m*) equal to the number of transitions found that started at startm and ended at stopm* and width
(we chose
2 = 0.0005 empirically to yield clear plots). Summing each Gaussian over all state combinations (startm, stopm*) in all traces yields the TDP, with the horizontal axis corresponding to the starting FRET values and the vertical axis corresponding to the final FRET values:
![]() | (6) |
The advantage of this method can be seen in the example of Fig. 4. The simulated FRET trajectory in Fig. 4 A shows no peaks when a histogram is constructed out of its FRET values (Fig. 4 B). But if there are distinct, reproducible FRET values from trace to trace, then peaks should develop in the TDP (Fig. 4 C). For a general S-state system with exchange possible between every state, S x (S 1) peaks should appear.
|
![]() | (7) |
hetero and negative diagonal width
homo.
|
|
![]() | (8) |
![]() | (9) |
i
S) FRET values of the S different states, the
hetero peak width in the positive diagonal direction, the
homo peak width in the negative diagonal direction, and ampi,j (1
i
S, 1
j
S, i
j), the normalized peak amplitudes of the S x (S 1) different peaks. Using Brent's algorithm again we vary these parameters until the probS is maximized for different choices of S, the number of underlying states. To determine which S is best we use the Bayesian information criterion (BIC) (22
![]() | (10) |
We can then determine BIC(S) for different S values until a minimum is reached; at that point the value that minimizes BIC(S) is the most likely number of underlying FRET states. For the data in Fig. 5 A, such a minimum was reached with the five-state prob5(x,y) overlaid in Fig. 5 B. With probS(x,y) known, one can easily determine which of the underlying (si,sj) peaks each detected (startm, stopm*) transition belongs to (if any). Once this list is compiled, each (si,sj) peak's mean transition probability and associated error are determined from Eqs. 4 and 5 above.
| ROBUSTNESS OF THE ALGORITHM |
|---|
|
|
|---|
|
FRET
|FRETB FRETA| = 0.4, FRET noise width
of 0.144, a typical value found in real data, and transition probabilities of tpA
B = 0.05 and tpB
A = 0.02. Each parameter was varied, whereas the others remained constant so as to determine the impact of that parameter, and 100 traces were analyzed for each choice of parameters. The success of the algorithm was measured in two ways. First, we determined what fraction of the 100 traces returned the true FRET values (0.3 (±0.05) or 0.7 (±0.05)) (solid squares in Fig. 6). Second, the deduced transition probability k* was compared to the true value k and |ln(k/k*)| plotted (open squares in Fig. 6).
Fig. 6 A shows the effect of changing
FRET, the spacing between the two idealized FRET values. The algorithm responds well with <40% error in transition probability and close to 100% yield of returning correct FRET values down to a
FRET of 0.1. Fig. 6 B shows the effect of FRET noise, parameterized by
. The algorithm responds well to increased noise, breaking down only when
> 0.4. Fig. 6 C shows the effect of data integration time relative to the state lifetimes. Here we vary the transition probabilities tpA
B and tpB
A to change their absolute values while preserving their ratio. Since tpA
B > tpB
A, we plot the algorithm response against the ratio of the data integration time and the dwell time of state A. Adequate state detection and transition probability determination (11% error) can be made even when the data integration time is only 60% longer than the state dwell time.
Effect of trace length
So far we have considered the cases wherein the observation time window was practically infinite, but in real experiments, the observation time is limited by photobleaching. Fig. 6 D shows the algorithm performance versus the average number of transitions per trace. For the standard transition probabilities used here,
100 data points are necessary to determine their values within 20%, which corresponds to only 1.5 transitions per trace on average.
Number of traces needed
Up to this point all discussions of error have related to systematic error, i.e., unavoidable error that could not be reduced by analyzing more traces. We found that even with as few as 20 traces, a number routinely obtained in single molecule experiments, the statistical contribution to error is minimal, at only 8% for the standard parameters (data not shown).
More complex systems
In moving from simple two- to three- or more-state systems, results will follow those described for the three-state system as long as (A), the FRET states are separated by at least the same distance as the noise width of FRET states; (B), the sampling rate is greater than the transition rate for all states; and (C), the system transits between each of the states more than once per trace. Of these criteria, by far the hardest to satisfy as the system becomes more complex is the last one. Since for every state added there are 2 (S 1) new entries in the transition probability matrix it becomes difficult to obtain traces that show transitions from every state to every other state. Fortunately in many systems not all transitions are possible or some transitions are very rare, allowing one to neglect several entries in the transition probability matrix and significantly reduce the requirements. For example, Fig. 4 shows the transition density plot (TDP) analysis of a simulated five-state data set of 400 traces, 500 data points/trace, FRET noise width of 0.144, and FRET values of 0.18, 0.35, 0.48, 0.63, and 0.78. Transition probabilities ranged from 0.002 to 0.1 and were chosen to favor transitions between neighboring states. The lowest peaks showing a small number of transitions have poor statistics and generate transition probabilities that deviate significantly from the true values. For the 10 highest peaks showing a significant number of transitions, the calculated transition probabilities agreed with the true values within 23%.
Heterogeneous broadening
In real single molecule FRET trajectories different molecules often exhibit different FRET values, even when observed over a sufficiently long time to remove statistical noise. This is particularly noticeable in wide field measurements when looking at trajectories obtained from different image files. This heterogeneous broadening is usually due to instrumental sources such as imperfect alignment of the donor and acceptor detection channels, heterogeneous background on the slide surface, and changes in microscope focusing, and is not necessarily indicative of an underlying molecular heterogeneity. Since this broadening is usually uniform, the result is to shift all FRET values by a nearly constant amount, largely preserving spacing. To test if our algorithm works well even in the presence of heterogeneous broadening, we modified the five-state system described above to exhibit heterogeneous broadening in FRET with a width of 0.15 (typical value obtained experimentally). The resulting TDP is found in Fig. 7. Using the thresholds plotted, the transition rates were calculated, and for the largest 10 peaks (those showing an appreciable number of transitions) the calculated values agreed with the underlying values with a typical error of 15%.
|
| APPLICATION TO EXPERIMENTAL DATA |
|---|
|
|
|---|
RecA binding and dissociation
Finally we applied our algorithm to experimental data that could not otherwise be analyzed reliably: the tracking of individual RecA protein binding and dissociation events on a single DNA molecule. Detailed experimental procedures and analysis are published elsewhere (23
), but in brief a construct was designed which would exhibit a change in FRET upon binding of a RecA monomer (Fig. 8 A). As more RecA proteins bind the FRET continues to decrease. There was only an initial guess as to the maximum number of RecA monomers that can bind to the DNA, and for that reason the analysis made no assumption about the number of states. Rather, the algorithm was asked to attempt to find a large enough number of states, 10, in the data.
Typical data with fit can be seen in Fig. 8 B. The resulting TDP (Fig. 8 C) constructed from raw, uncorrected intensity data shows clearly resolved peaks at 0.15, 0.3, 0.45, 0.6, and 0.8 along each axis. The first, centered at approximately FRET = 0.15 is the result of acceptor blinking, and is not indicative of a true physical state. The remaining four were presumed to be different legitimate FRET states of the system. The 0.8 peak coincides with the value obtained in the absence of RecA and the 0.3 peak coincides with the value obtained with RecA in the presence of ATP
S (where RecA binds stably), leading to the conclusion that the two remaining peaks (0.45 and 0.6) are the result of an intermediate number of RecA molecules bound. Since there are three FRET peaks other than the DNA-only peak indicating that up to three RecA monomers can bind, we label the FRET peaks accordingly (0.8
M0, 0.6
M1, 0.45
M2, and 0.3
M3). It is worth noting that the TDP shows transitions occurring primarily between neighboring FRET values, suggesting that RecA binds and dissociates as a monomer. With states defined we proceeded to determine the state-to-state transition rates at various RecA concentrations and graphed the results in Fig. 9. Note that the transition probability of going from M0 to M1, thus for binding, increases significantly as higher concentrations of RecA are added (Fig. 9 A). No corresponding change was seen for dissociation for example for the transition probability of going from M3 to M2 (Fig. 9 B). We conclude that our algorithm can analyze a complex four-state system without any preconceived model and can return nontrivial conclusions such as transitions being predominantly between nearest neighbors.
|
| SUMMARY |
|---|
|
|
|---|
Both the HMM software and the transition density plotting software are publicly available and can be downloaded from our website at http://bio.physics.uiuc.edu/HaMMy.html.
| ACKNOWLEDGEMENTS |
|---|
|
|
|---|
Submitted on February 1, 2006; accepted for publication May 31, 2006.
| REFERENCES |
|---|
|
|
|---|
2. McKinney, S. A., A. D. Freeman, D. M. Lilley, and T. Ha. 2005. Observing spontaneous branch migration of Holliday junctions one step at a time. Proc. Natl. Acad. Sci. USA. 102:57155720.
3. Nahas, M., T. J. Wilson, S. C. Hohng, K. Jarvie, D. M. J. Lilley, and T. Ha. 2004. Observation of internal cleavage and ligation reactions of a ribozyme. Nat. Struct. Mol. Biol. 11:11071113.[CrossRef][Medline]
4. Okumus, B., T. J. Wilson, D. M. J. Lilley, and T. Ha. 2004. Vesicle encapsulation studies reveal that single molecule ribozyme heterogeneities are intrinsic. Biophys. J. 87:27982806.
5. Rueda, D., G. Bokinsky, M. M. Rhodes, M. J. Rust, X. W. Zhuang, and N. G. Walter. 2004. Single-molecule enzymology of RNA: essential functional groups impact catalysis from a distance. Proc. Natl. Acad. Sci. USA. 101:1006610071.
6. Tan, E., T. J. Wilson, M. K. Nahas, R. M. Clegg, D. M. Lilley, and T. Ha. 2003. A four-way junction accelerates hairpin ribozyme folding via a discrete intermediate. Proc. Natl. Acad. Sci. USA. 100:93089313.
7. Zhuang, X., H. Kim, M. J. Pereira, H. P. Babcock, N. G. Walter, and S. Chu. 2002. Correlating structural dynamics and function in single ribozyme molecules. Science. 296:14731476.
8. McKinney, S. A., A. C. Declais, D. M. J. Lilley, and T. Ha. 2003. Structural dynamics of individual Holliday junctions. Nat. Struct. Biol. 10:9397.[CrossRef][Medline]
9. Eddy, S. R. 2004. What is a hidden Markov model? Nat. Biotechnol. 22:13151316.[CrossRef][Medline]
10. Qin, F., A. Auerbach, and F. Sachs. 2000. A direct optimization approach to hidden Markov modeling for single channel kinetics. Biophys. J. 79:19151927.
11. Qin, F., A. Auerbach, and F. Sachs. 2000. Hidden Markov modeling for single channel kinetics with filtering and correlated noise. Biophys. J. 79:19281944.
12. Yang, H., and X. S. Xie. 2002. Statistical approaches for probing single-molecule dynamics photon-by-photon. Chem. Phys. 284:423437.[CrossRef]
13. Yang, H., G. B. Luo, P. Karnchanaphanurach, T. M. Louie, I. Rech, S. Cova, L. Y. Xun, and X. S. Xie. 2003. Protein conformational dynamics probed by single-molecule electron transfer. Science. 302:262266.
14. Andrec, M., R. M. Levy, and D. S. Talaga. 2003. Direct determination of kinetic rates from single-molecule photon arrival trajectories using hidden Markov models. J. Phys. Chem. A. 107:74547464.[CrossRef]
15. Schroder, G. F., and H. Grubmuller. 2003. Maximum likelihood trajectories from single molecule fluorescence resonance energy transfer experiments. J. Chem. Phys. 119:99209924.[CrossRef]
16. Ha, T., I. Rasnik, W. Cheng, H. P. Babcock, G. H. Gauss, T. M. Lohman, and S. Chu. 2002. Initiation and re-initiation of DNA unwinding by the Escherichia coli Rep helicase. Nature. 419:638641.[CrossRef][Medline]
17. Zhuang, X., L. E. Bartley, H. P. Babcock, R. Russell, T. Ha, D. Herschlag, and S. Chu. 2000. A single-molecule study of RNA catalysis and folding. Science. 288:20482051.
18. Dahan, M., A. A. Deniz, T. J. Ha, D. S. Chemla, P. G. Schultz, and S. Weiss. 1999. Ratiometric measurement and identification of single diffusing molecules. Chem. Phys. 247:85106.[CrossRef]
19. Viterbi, A. J. 1967. Error bounds for convolutional codes and an asymptotically optimum decoding algorithm. IEEE Trans. Inf. Theory. 13:260267.[CrossRef]
20. Press, W. H., S. A. Teukolsky, W. T. Vetterling, and B. P. Flannery. 1992. Numerical recipes in C. Cambridge University Press, Cambridge.
21. Baum, L. E., T. Petrie, G. Soules, and N. Weiss. 1970. A maximization technique occurring in the statistical analysis of probabilistic functions of Markov chains. Ann. Math. Stat. 41:164171.[CrossRef]
22. Hastie, T., R. Tibshirani, and J. Friedman. 2001. The Elements of Statistical Learning: Data Mining, Inference, and Prediction. Springer-Verlag, New York.
23. Joo, C., S. A. McKinney, M. Nakamura, I. Rasnik, S. Myong, and T. Ha. 2005. Real time observation of RecA filament dynamics with single monomer resolution. Cell. In press.
24. Eddy, S. R. 1998. Profile hidden Markov models. Bioinformatics. 14:755763.
This article has been cited by other articles:
![]() |
F. Vanzi, L. Sacconi, and F. S. Pavone Analysis of Kinetics in Noisy Systems: Application to Single Molecule Tethered Particle Motion Biophys. J., July 1, 2007; 93(1): 21 - 36. [Abstract] [Full Text] [PDF] |
||||
![]() |
M. Linden and M. Wallin Dwell Time Symmetry in Random Walks and Molecular Motors Biophys. J., June 1, 2007; 92(11): 3804 - 3816. [Abstract] [Full Text] [PDF] |
||||
![]() |
J.-D. Wen, M. Manosas, P. T. X. Li, S. B. Smith, C. Bustamante, F. Ritort, and I. Tinoco Jr. Force Unfolding Kinetics of RNA Using Optical Tweezers. I. Effects of Experimental Variables on Measured Results Biophys. J., May 1, 2007; 92(9): 2996 - 3009. [Abstract] [Full Text] [PDF] |
||||
![]() |
J. F. Beausang, C. Zurla, C. Manzo, D. Dunlap, L. Finzi, and P. C. Nelson DNA Looping Kinetics Analyzed Using Diffusive Hidden Markov Model Biophys. J., April 15, 2007; 92(8): L64 - L66. [Abstract] [Full Text] [PDF] |
||||
![]() |
Y. Zhou and X. Zhuang Robust Reconstruction of the Rate Constant Distribution Using the Phase Function Method Biophys. J., December 1, 2006; 91(11): 4045 - 4053. [Abstract] [Full Text] [PDF] |
||||
| ||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
| HOME | HELP | FEEDBACK | SUBSCRIPTIONS | ARCHIVE | SEARCH | TABLE OF CONTENTS |