## 1.

## Introduction

The Förster resonance energy transfer (FRET) is a radiationless process between two fluorophores, donor and acceptor, whose intensities ratio defines transfer effeciency ($E$) and reports the in-between distance $\{E={[1+{(R/{R}_{0})}^{6}]}^{-1}\}$, where ${R}_{0}$ is the Förster distance between the donor and the acceptor. As a single biological molecule is labeled by a FRET pair and immobilized,^{1} its conformation dynamics associated with the function can be recorded by using a two-channel fluorescence microscope to track FRET changes for an extended period of time.^{2} In practice, this type of single-molecule FRET (smFRET) microscope can be realized through a confocal configuration,^{3} by which one molecule is imaged at a time, or a wide-field configuration that allows for hundreds of molecules to be simultaneously monitored by a pixelated detector. To achieve the detection of single-molecule fluorescence in a wide-field configuration, total internal reflection (TIR) has been employed to generate an evanescent wave to excite a thin layer near the interface such that the fluorescence background from the bulk can be greatly reduced.^{1}^{,}^{3}4.5.^{–}^{6}

The wide-field smFRET data are embedded in movies of $N$ time-frame of CCD images. Each image is divided into two half-images: one of the donor (usually Cy3, TMR, Alexa, or Atto 555) and the other of the acceptor channel (Cy5, Alexa, or Atto 647). To maximize the temporal resolution and achieve the longest single-molecular time traces before photo-bleaching occurs, one should increase the CCD frame rate as fast as possible and keep the laser illumination power as low as possible, yielding the time trace data with poor signal-to-noise ($\mathrm{S}/\mathrm{N}$) ratio. A wide-field smFRET microscopy can be realized by a prism type TIR (PTIR)^{1} or a high-numerical aperture (NA) objective type (OTIR) (Fig. 1).^{4}^{,}^{6} We chose the OTIR as it is easy to build and offers an empty space on top of the sample slide that would allow for sample manipulation from above. Compared to PTIR, OTIR requires an additional dichroic mirror underneath the objective to direct the excitation beam into the objective lens and yet to prevent the reflected excitation beam from entering the pathway of the fluorescence collection optics. The photon collection is thus comprised in the OTIR configuration despite having a high NA. We have experienced higher fluorescence background that might come from the objective lens or its associated elements as reported.^{7}^{,}^{8} All those conditions have together yielded very noisy data so that it becomes very challenging to studying hundreds of thousands of noisy time traces to select meaningful ones by human efforts. To facilitate the data reduction from smFRET movies to smFRET traces and address the noise-limited issues, we developed a denoising recipe that utilizes novel statistical approaches based on the spatial and/or temporal correlation at different steps in the work flow. Due to these new algorithms, the work flow can be automated.

As to the estimation and removal of the background, local subtraction^{5}^{,}^{9}^{,}^{10} and profile fitting^{11} are two commonly used methods, both of which utilize local information around the fluorophores. Interestingly, as we introduced a global method based on higher-order singular value decomposition (HOSVD) to extract spatial and temporal features in a movie for denoising the system errors, we found that it could estimate the non-uniform TIR background quite well. As to mapping and registering the fluorophore coordinates in the two channels, we used temporal correlation to locate those spots representing the donor signals bleeding through the acceptor channel and then used their coordinates to determine the best linear transformation matrix. We found the transformation matrix was not adequate because non-uniform deviations on the coordinate mapping were observed across the whole image. To identify the pairs, our algorithm performs a spatially exhaustive search in the neighboring pixels of the spot coordinates predicted by the matrix. Once the trace pair was found, we checked whether it was a FRET trace by studying if there were donor and acceptor anti-correlated patterns along the pair traces. For those FRET traces, we used change point detection (CPD) to detect the status change points and reported the dwelling time. To do so, we derived a likelihood ratio statistic given a change point position under a multivariate Gaussian model and used it to search all time points in the interval between adjacent change points to find the new change point candidates. Conventionally, the status change points are solved by hidden Markov model (HMM),^{14} multivariate Gaussian HMM (MGHMM),^{15} or time order clustering (TOC).^{16} In contrast to those methods, CPD is a deterministic approach that does not require assigning initial values such as state means and variance. As such, our work flow employing the CPD algorithm could be executed without human intervention.

This paper is organized as follows: Section 2 introduces the analysis methods, and these methods are summarized in five algorithms. In Sec. 3 we present our real data analysis and simulation studies. This paper ends with a brief discussion.

## 2.

## Methods

The main goal of the FRET data analysis is to target the FRET traces and detect the status change points along the trace. Our strategy includes three steps, which are summarized in a flow chart in Fig. 2. The first step is to denoise the system error. We applied HOSVD to estimate the non-homogenous system error. Traditionally, the corresponding approach is to perform local subtraction of the intensities in the surrounding pixels of a spot. These two approaches are shown to match quite well on both simulated data and real data. One extra benefit from the HOSVD approach is that the estimated system error could give feedback to the experimenters so as to optimize their instrument setup.

The second step is to find the paired FRET traces. In order to make it automated, three sub-steps are executed by applying three algorithms. They involve 1. locating the fluorophores from the acceptor images; 2. mapping the coordinate system of the donor images and acceptor images; 3. finding all possible fluorophore paired traces.

The third step is to detect the change points along the FRET trace pairs. Instead of detecting the change points along the FRET trace resulting from the ratio of the acceptor to the combination of the acceptor and donor, we treated the donor and the acceptor traces as multivariate variables and detected the change points along them. To detect the change points, we adopted the multivariate Gaussian model and used the log-likelihood ratio as our statistics.^{12}^{,}^{13} We also checked whether those traces are locally anti-correlated across time. The detection criterion is based on the $p$-value. Once the threshold (given a $p$-value) is determined, the computation is deterministic, involving no random initials, which includes the state mean and variance values or pre-specified number of states as in commonly used methods such as HMM,^{14} MGHMM,^{15} and TOC.^{16}

To test our algorithm, both simulation and real data were used. The sample used to generate real data shown in Figs. 3 to 6 was a 16 bp GC-rich double strand DNA with Cy3 and Cy5 attached at the 5’ ends modified from a design previously described^{17} to allow for specific binding of a protein, of which the details will be published in a separate paper. To compare the performance of CPD with MGHMM (Fig. 7), the movie data were provided by Ha’s laboratory website ( http://www.cplc.illinois.edu), and the sample is a DNA plus its binding protein.

## 2.1.

### System Error Denoising

One particular feature of the smFRET experiment is that the fluorescent molecules are immobilized. Thus, we proposed a local model: let $({i}_{0},\text{\hspace{0.17em}}{j}_{0})$ be the center of a fluorophore, then the signal can be modeled as

## (1)

$${I}_{{i}_{0},{j}_{0}}^{L}(i,j,k)={{F}^{\mathrm{ps}}}_{{i}_{0},{j}_{0}}(i,j)\xb7{g}_{{i}_{0},{j}_{0}}(k),\phantom{\rule[-0.0ex]{1.0em}{0.0ex}}\phantom{\rule{0ex}{0ex}}\mathrm{for}\text{\hspace{0.17em}}\text{\hspace{0.17em}}1\le i\le I,1\le j\le J,1\le k\le K,$$^{18}${I}^{G}(i,j,k)={\sum}_{m,n,l}{a}_{m,n,l}{f}_{x,m}(i){f}_{y,n}(j){f}_{z,l}(k)$, where each ${f}_{\xb7,\xb7}$ represents one linear mixture on one mode of the array and each product term of ${f}_{x,\xb7}{f}_{y,\xb7}{f}_{z,\xb7}$ represents one component. ${a}_{\xb7,\xb7,\xb7}$ refers to the projected size of the data on the corresponding component. The solution to this decomposition problem can be obtained by the alternating least square algorithm,

^{19}

^{,}

^{20}described in Algorithm 1. It may be helpful to link this model with the commonly used SVD (singular value decomposition), which can be viewed as a two-order version of this model.

Most data analysis would model signal components with larger variance (eigenvalue) and noise with smaller variance; however, this is not so in our case. In real data analysis, we found the first component with the largest eigenvalue describes the system error $B(\xb7,\xb7,\xb7)$ very well. The data structure indicates that this result would not be a surprise. Each signal ${I}_{\xb7,\xb7}^{L}(\xb7,\xb7,\xb7)$, localized within 5 by 5 pixels, would unlikely produce large variance. On the other hand, the system error, representing global information, can contribute a large portion of the total variance. The phenomenon that the first component captures the global information has also been reported on an integrative analysis of DNA microarray set combining different studies.^{21} To be specific, we used ${a}_{1,1,1}{f}_{x,1}(i){f}_{y,1}(j){f}_{z,1}(k)$ to estimate $B(i,j,k)$ and delete it from the data as a denoising process. We further developed Algorithm 2 to smooth this term in order to avoid the possible mixture of the signals. Fig. 3 provides a typical example.

We also designed simulation experiments to check the performance of this algorithm. The result is very supportive. We compared this algorithm with the concurrent local background subtraction method. This part is reported in section 3.1.

## 2.2.

### Extract FRET Traces

## 2.2.1.

#### Locate fluorophores

The fluorophores we are interested in have two characters: they are immobilized, and they last for a period of time. Therefore, the corresponding pixel values are supposed to be comparatively high over a range of spatial index and over a period of time. The image with the system error denoised is transformed to a binary image, depending on whether the pixel value exceeds a given threshold. The threshold is set to be the overall mean plus only one standard deviation to avoid missing the candidate fluorophore cluster. Subsequently, the clusters with the value of one are searched throughout the spatial index and the time frame index as well. The detailed procedure is stated in Algorithm 3.

## 2.2.2.

#### Image registration (map the coordinate systems)

Traditionally, the mapping of coordinates relies on imaging a calibration slide made of fluorescence beads of broad spectrum before collecting the data from the experimental slides as a daily practice.^{5}^{,}^{6} The fluorescence spots in the donor channel and their leakage pairs in the acceptor channel are used for the registration of the two coordinate systems by solving a transformation matrix. Alternatively, we took an approach that used the leakage spots found in the experimental images to determine the mapping between the two channels so that the work on the fluorescence beads could be omitted. Leakage spots refer to those donor-only fluorophores carrying abnormal intensities so that they bleed through into the acceptor channel. As a result, the values of those leakage pixels increase or decrease simultaneously in both channels. Having a spot located in the acceptor channel by the previous sub-step, we looked for its leakage pair in the donor channel by searching in the neighborhood of the corresponding spot coordinates. (Note that a simple translation relationship between the acceptor and the donor coordinate systems represents a good approximation to begin with when the two channels are aligned roughly parallel to each other.) When a set of leakage pairs was found, we applied least squared error method to determine the transformation matrix. This algorithm is briefly described in Algorithm 4.

## 2.2.3.

#### Search fluorophore paired traces

Given the acceptor fluorophores and the coordinate transformation matrix, we applied Algorithm 3 again to search all possible donor fluorophores on the corresponding locations. In this part, the traditional method searches the candidate spots with good $S/N$ at the corresponding locations in the donor channel by averaging the initial, say, 10 time frames. With this regard, we adopt a temporally exhaustive search along the time trajectory to include all possible fluorophores. The necessity comes from considering the following scenario in which the traditional approach will likely miss a donor spot while the greedy search will not: the FRET occurs for a period longer than 10 time frames in the beginning, such that no fluorophore would be detected in the donor channel based on the initial average.

## 2.3.

### Change Point Detection

## 2.3.1.

#### Check if FRET occurs

Given a fluorophore trace pair, we check if there exists any real FRET event by anti-correlations. One means to do so is to calculate the correlation over the whole trace pair, which can be easily automated. However, this approach would sometimes generate an artifact when a time trace carries a dark state period and/or a long range of random noise. Instead, the human inspection approach usually focuses on the local anti-correlation around the candidate change point. To mimic human inspection, we applied a change point method to locate the possible change points and calculated the local correlations around them for checking if FRET occurs. Our approach at this step can greatly shrink the size of the candidate set.

## 2.3.2.

#### Detecting change points

Change point detection is a statistical method to detect different states in a sequential data set. A very brief introduction on CPD is presented in the Appendix. For the theoretical background, please refer to Siegmund^{12} and Chen and Gupta.^{13} Previously, Watkins and Yang^{22} applied a likelihood ratio test to locate the intensity jumps as the change point based on individual photon arrival times in a single molecule trace extracted from one channel. In contrast, the CPD used here is based on multivariate Gaussian model that allows us to simultaneously analyze the paired intensity traces.^{23} Fig. 6 is a typical example. The tuning parameter for this approach is the significant size (the threshold for $p$-value). Reasonable $p$-value is in the range between 0.1 and $0.01:.1$ (more flexible) or 0.05 (marginal) or 0.01 (for the sake of multiple testing). The computation is deterministic, such that it involves neither random initials nor a pre-specified number of states as those commonly used methods like HMM,^{14} MGHMM,^{15} or TOC.^{16} which not only saves the computation time but also minimizes human intervention.

## Algorithm 1

### Alternating least square algorithm^{19}.

Given initial vectors ${f}_{x}^{0}$, ${f}_{y}^{0}$, ${f}_{t}^{0}$. Denote ${f}_{x,i}$ as $i$’th element of vector ${f}_{x}$.

**input** $\mathrm{I}(\mathrm{i},\mathrm{j},\mathrm{k})$ {# input movies of smFRET data}

**for** $\mathrm{p}=0,1,\dots \mathrm{P}$ **do** {# repeat $\mathrm{P}$ iterations}

**for** $\mathrm{i}=1,\dots \mathrm{I}$ **do**

${f}_{x,i}^{p+1}={\sum}_{j,k}I(i,j,k){f}_{y,j}^{p}{f}_{t,k}^{p}$

**for** $\mathrm{j}=1,\dots \mathrm{J}$ **do**

${f}_{y,j}^{p+1}={\sum}_{i,k}I(i,j,k){f}_{x,i}^{p}{f}_{t,k}^{p}$

**for** $\mathrm{k}=1,\dots \mathrm{K}$ **do**

${f}_{t,k}^{p+1}={\sum}_{i,j}I(i,j,k){f}_{x,i}^{p}{f}_{y,j}^{p}$

normalize so that ${\Vert {f}_{x}^{p+1}\Vert}_{2}={\Vert {f}_{y}^{p+1}\Vert}_{2}={\Vert {f}_{t}^{p+1}\Vert}_{2}=1$

**output** ${f}_{x}^{P}$, ${f}_{y}^{p}$ and ${f}_{t}^{p}$

## Algorithm 2

### Background estimation

**input** $\mathrm{I}(\mathrm{i},\mathrm{j},\mathrm{k})$, ${f}_{x}$, ${f}_{y}$.

${f}_{x}$ and ${f}_{y}$ are outputs of Alogrithm 1.

${\tilde{f}}_{x}$, ${\tilde{f}}_{y}$ are degree two local polynomial regressions of ${f}_{x}$, ${f}_{y}$ with span 0.2.

${\tilde{f}}_{t,k}={\sum}_{i,j}I(i,j,k){\tilde{f}}_{x,i}{\tilde{f}}_{y,j}$.

**output** ${({B}_{\mathit{HOSVD}})}_{i,j,k}={\tilde{f}}_{x,i}{\tilde{f}}_{y,j}$.

## Algorithm 3

### Localization of fluorophore

Let ${D}_{l}({i}_{0},{j}_{0})=\{(i,j)\}|i-{i}_{0}|\le l,|j-{j}_{0}|\le l\}$ be the square index centered at $({l}_{0},{j}_{0})$ with edge size $2l+1$ and $|A|$ be the size of the set $A$.

**input** $\mathrm{I}(\mathrm{i},\mathrm{j},\mathrm{k})$. Denote the $k$’th frame of $I$ by ${I}_{k}$.

${\mu}_{k}=\mathit{mean}({I}_{k}),{\sigma}_{i}=\mathit{std}({I}_{k})$.

**for** $\mathrm{k}=1,\dots ,\mathrm{K}$ **do**

get ${E}_{k}=[(i,j)|I(i,j,k)>{\mu}_{k}+{\sigma}_{k}]$. {# select high-intensity pixels}

**for** ${i}_{0},{j}_{0}$ in ${E}_{k}$ **do**

**if** $|{D}_{1}({i}_{0},{j}_{0})\cap {E}_{k}|<5$, eliminate $({i}_{0},{j}_{0})$ from ${E}_{k}$. {# keep aggregated high-intensity pixels}

$E={\cup}_{r=1}^{K-W}{\cap}_{0\le l<w}{E}_{r+l}$ for some chosen $W$ {# keep sustained and aggregated high-intensity pixels}

**for** ${i}_{0}$, ${j}_{0}$ in $E$ **do**

${A}_{{i}_{0},{j}_{0}}(k)=\mathit{argmax}[I(i,j,k)|(i,j)\in {D}_{2}({i}_{0},{j}_{0})]$

$C({i}_{0},{j}_{0})={\mathit{mode}}_{k}({A}_{{i}_{0},{j}_{0}}(k))$

**if** $C({i}_{0},{j}_{0})\ne ({i}_{0},{j}_{0})$, eliminate $({i}_{0},{j}_{0})$ from $E$.

**output** $E$. {# select a local maximum as a representative}

## Algorithm 4

### Coordinate transformation for two channels

Denote the time trace in donator and acceptor channels at $(i,j)$ point by ${I}^{D}(i,j,k)$ and ${I}^{A}(i,j,k)$, respectively.

**input** ${I}_{k}^{A}$ into Algorithm 3 with $E$ as its output. {# get the coordinate indices for acceptor fluorophores}

Let ${E}_{A}$ and ${E}_{D}$ be two empty sets.

**for** ${i}_{0}$, ${j}_{0}$ in $E$ **do**

**for** $(r,s)\in {D}_{4}({i}_{0},{j}_{0})$ **do**

${\rho}_{r,s}=\text{correlation coefficient}$ of ${I}^{A}({i}_{0},{j}_{0},\xb7)$ and ${I}^{D}(r,s,\xb7)$.

$({r}^{*},{s}^{*})={\mathit{argmax}}_{r,s}{\rho}_{r,s}$.

**if** ${\rho}_{{r}^{*},{s}^{*}}>0.8$, {# search the possible leakage pairs around $({i}_{0},{j}_{0})$}

$(i,j)\to {E}_{A}$, $({r}^{*},{s}^{*})\to {E}_{D}$ {# register the donor and acceptor coordinates for leakage pairs}

use least square method to find $A$ which map the index $({r}^{*},{s}^{*})\in {E}_{D}$ to its corresponding index $({i}_{0},{j}_{0})\in {E}_{A}$.

**output** the transformation matrix $A$.

## Algorithm 5

### Change point detection.

Initially, let the set of change points $C$ be $\{1,N\}$. Denote $({X}_{D},{X}_{A})$ as bi-variate traces of donor, and acceptor ${S}_{r}$ is defined in the Appendix.

**repeat**

Sort $C=\{{c}_{1},{c}_{2},\dots ,{c}_{n}\}$ such that ${c}_{1}<{c}_{2}<\dots <{c}_{n}$.

**for** $\mathrm{i}=1,\dots ,\mathrm{n}$ **do**

${r}^{*}={\mathit{argmax}}_{r\in ({c}_{i},{c}_{i+1})}{S}_{r}$

**if** ${S}_{{r}^{*}}>$ threshold $({c}_{i+1},{c}_{i})$ and $\mathit{min}({r}^{*}-{c}_{i},{c}_{i+1}-{r}^{*})>5$

${r}^{*}\to C$ {# find the possible change point in interval $({c}_{i},{c}_{i+1})$ and update the change point set $C$.

**until** no more ${r}^{*}$ is significant in each interval

**for** $\mathrm{i}=1,\dots ,\mathrm{n}$ **do**

${\rho}_{i}=\text{correlation coefficient}$ of $({X}_{D,j},{X}_{A,j})$,${c}_{i}-20\le j\le {c}_{i}+20$.

**if** there exists an $i$ such that ${\rho}_{i}<-0.5$,}

**then** label $({X}_{D},{X}_{A})$ as a FRET candidate. {# check an anti-correlated pattern around change points}

Collect mean values and dwelling time in each interval $[{c}_{i},{c}_{i}+1]$ for all FRET candidate to plot a histogram.

## 3.

## Real Data Analysis and Simulations

## 3.1.

### Results of Background Subtraction

To investigate the performance of HOSVD in estimating the system error, we simulated the data according to model 2. Poisson random variables ${X}_{i,j,k}$ are generated with conditional mean: $B(i,j)+\sum _{l=1}^{100}1000\xb7{{F}^{ps}}_{{i}_{l},{j}_{l}}(i,j){g}_{{i}_{l},{j}_{l}}(k)$, $1\le i\le 256$, $1\le j\phantom{\rule{0ex}{0ex}}\le 512$ and $1\le k\le 40$. We generated the system error $B(i,j)$, simulating real data; thus, we fitted a bi-variate function $B(i,j)$ on a real image as

^{5}($l\times l$ LS, $2l+1$ refers to the size of the background) at $({i}_{0},{j}_{0})$, which can be specified as

## Table 1

Errors in background estimation. The simulation setup is in Section 3.1. Note that the sever maximum error in 5×5 LS occurs when averaging over a region containing other fluorophores. Since fluorophores locate at least 6 pixels apart, the 3×3 LS is free from this concern.

Maximum error | Root mean square error | |
---|---|---|

Algorithm 1 only | 36.0±4.9 | 15.6±1.1 |

Algorithms 1 and 2 | 13.6±1.6 | 5.1±0.5 |

3×3 LS | 19.5±1.6 | 5.7±0.1 |

5×5 LS | 57.4±24.3 | 8.3±3.3 |

We further compared the performance between HOSVD and LS on real data sets. In this comparison, the size of system error is unknown. We applied both methods to estimate system error underlying 134 manually selected time traces. Fig. 4 shows a linear relationship between estimated values of HOSVD and LS, indicating compatibility between the two methods. Particularly, the estimated system error by HOSVD is most consistent with those from $7\times 7$ LS whereas $3\times 3$ LS usually overestimates the system error. Choosing the size $l$ in local subtraction method is an art: larger $l$ may reduce the point spread function impact but would have the risk of entering another fluorophore’s influence circle.

## 3.2.

### Results of Searching FRET Trace Pairs

Fig. 5 shows an image example of the leakage pairs to demonstrate its feasibility on mapping the coordinates. Fig. 6 gives an example of the searched FRET trace pair. Generally, all the single-molecule FRET analysis methods have the first three sub-steps in common: target fluorophores, perform image registration, and search for the trace pairs. Thus, we set up a task to check the performance of searching the FRET trace pairs from all the fluorophore trace pairs. There were 4421 trace pairs from an experiment data set consisting of 14 movies passing through our algorithm. We let an experimenter who is experienced and rigorous with regards to manual selection screen the FRET traces. The manual selection is based on a more stringent criteria listed as follows:

1. the traces exhibit single-step photo-bleaching

2. the average fluorescence intensity along the trace is constant

3. the traces have a normal signal strength

4. the acceptor channel is fluorescent

Thus, the experimenter selected only 79 FRET trace pairs out of the pool of 4421 (a typical example of such unambiguous FRET traces is shown in Fig. 6). On the other hand, by using the FRET screening step in our algorithm, more than 80% of the 4421 trace pairs were efficiently eliminated, yielding a subset of 583 candidate FRET trace pairs. Seventy-seven among the 79 manually chosen ones were preserved in the subset of 583. We checked those two traces that were selected manually but not by our algorithm to find they could be omitted as their anti-correlations, -0.4 and -0.2, were not significant. This result, summarized in Table 2, demonstrates that our algorithm can shrink the candidate set by a factor of approximately 8.

## Table 2

Results of selecting FRET trace pairs. There are 14 movies in this data set. The upper row shows the ratio of manually selected targets to the candidate pairs for each movie. The lower row shows the ratio of targets to those pairs selected by our algorithms.

1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10 | 11 | 12 | 13 | 14 | |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|

targetstotal traces | 10172 | 5204 | 11361 | 8285 | 11268 | 4323 | 3400 | 2335 | 2344 | 7329 | 5335 | 1331 | 3392 | 7342 |

selected targetsselected traces | 1040 | 529 | 1052 | 836 | 1044 | 428 | 335 | 241 | 252 | 750 | 544 | 128 | 355 | 749 |

## 3.3.

### Results of Change Point Detection

To compare the performance of our CPD method with that of MGHMM, we simulated a three-state system. While CPD needs the threshold for $p$-value, MGHHM requires inputs of random initials and the number of states. To be specific, let $({X}_{D},{X}_{A})$ follow a bi-variate three-state Markov model with parameters as follows:

1. Initial distribution $\pi =(1/3,1/3,1/3).$

2. Transition probability matrix

3. Observation distributions of $({X}_{D},{X}_{A})$ are two-dimensional Poisson distributions with (800,200), (500,500), and (520,480) as the mean values for the three states. Their corresponding standard deviations are (28.3, 14.1), (22.4, 22.4) and (22.8, 21.9).

In other words, we simulated a situation in which two compact conformations exist plus an unfolded state with distinct FRET efficiencies 0.5, 0.48, and 0.2, respectively. Furthermore, this transition matrix $A$ would allow the dwelling time to be distributed as a geometric random variable with mean 100 units. We would like to test whether or not these algorithms can distinguish these two conformers. Figs. 8(a) and 8(b) display some typical histograms based on 100 simulated time traces. Note that two populations of different FRET efficiency are distinguishable by both approaches. Surprisingly, as the difference was made smaller by letting the two FRET efficiencies be 0.5 and 0.49, our algorithm appeared to perform better with respect to the resolution [Figs. 8(c) and 8(d)]. The resolution of CPD is increased because it provides more efficient denoising by taking the average of the intensities over the dwelling time. The reason MGHMM fails in this case is because the likelihood does not gain as much as the paid penalty associated with introducing one more state.

We also studied the performance of our CPD algorithm by employing a real smFRET model data set provided by TJ Ha’s group at UIUC ( http://www.cplc.illinois.edu/). The data contains movies taken from two donor-acceptor labeled DNA samples. One is the DNA-only control, and the other is treated with a protein that can induce the change in FRET efficiency. A total of 15,548 raw time traces were processed by the UIUC’s package to select the FRET traces.^{5} Those traces were further studied by the MGHMM and the CPD method for comparisons. Fig. 7 displays the histogram results generated by the two methods. MGHMM and CPD both enhance the resolution in the histograms while CPD seems to give slightly sharper modes, consistent with the results obtained from the simulation data (Fig. 8). In summary, the comparison of the performance suggests that our CPD approach can be either better than or as good as the MGHMM in terms of removing the random noise. Most importantly, to perform MGHMM usually requires many tries to escape the trap of local minimum whereas there is no such need for the CPD approach.

## 4.

## Discussion

The development of our algorithms has involved implementation of several novel ideas. First, we utilized the knowledge on the smFRET system to model the data and applied HOSVD to clean the system error effectively. Second, we adapted to the leakage spots for registration of the donor-acceptor mapping coordinates. Finally, we employed the characteristic of coincidental change of the donor and acceptor traces to implement the multivariate change point analysis to detect the smFRET events. By using these approaches, automatic computing has replaced human manual efforts, largely if not completely, including the fluorescence bead slide setting,^{5}^{,}^{9}10.^{–}^{11} human inspection on choosing FRET traces, and parameter assignment for the micro-status change detection. Nevertheless, the full automation in the current version of our package falters at the step of selecting meaningful FRET traces because we have implemented only an anti-correlation rule using a less stringent condition in order not to miss any interesting traces. Furthermore, by module segmentation, our algorithm package allows for combination with other software packages. The sharpness of the CPD methods depends on the dwelling time in each state. When the dwelling time is too short, we may not see a clear boundary between states. Thus, our package also includes an option for MGHMM analysis, which applies a model selection criterion to set the boundary. For those interested in the details of various analysis methods for smFRET time trajectories^{14}^{,}^{24}^{,}^{25} and their comparisons, we recommend a comprehensive article by Bianco and Walter,^{26} which appeared while we were preparing this article.

A wide-field TIR equipped with a CCD has made it possible for collecting single-molecule FRET data in a high throughput manner. Recently, the smFRET method has been applied to various biological processes to reveal dynamic behaviors. These processes include catalytic RNA,^{10}^{,}^{27} polymerase-nucleic interactions,^{28}^{,}^{29} ribosome translation,^{30}^{,}^{31} spliceosome assembly,^{32}^{,}^{33} vesicle fusion,^{34} and the intrinsic protein disorder involved in the process.^{35} In addition to revealing the dynamics of biological molecules, smFRET can be employed to study the structure of biological molecules. However, most researchers in this field find it difficult to relate FRET efficiency to a physical distance. Since the anisotropy of the fluorophores’ dipoles can be characterized^{5} and the quantum efficiency of the two channels can be calibrated,^{36} it is now possible to convert smFRET to absolute distance. By obtaining a set of distances between different sites of a biological complex through a large number of smFRET measurements and triangle analysis, the partial structure of the complex might be determined.^{17}^{,}^{37}^{,}^{38} As the TIR solution, the dual-view splitter, and the high-sensitivity CCD are available from the market nowadays, it would be easy for an experimenter to set up a wide-field smFRET and generate a large volume of data. Interestingly, very often the time trace from a wide-field smFRET experiment is compared to that from the single ion-channel technique.^{39} However, unlike the single ion-channel recording, from which a time trace is reported in real-time, the time traces from wide-field smFRET experiment are obtained through post-imaging data processing. Therefore, to turn the wide-field smFRET microscope into a friendly tool with real-time data analysis potential, it is crucial to improve its data processing efficiency dramatically so that an experimenter can see the time traces immediately at the end of an experimental session to have a clue about what to do for the next. In order to facilitate the post-imaging processing for the smFRET movie data in an online manner, we have created a highly efficient analytical package that combines a series of algorithms that allow for full automation of the work flow of the smFRET data processing.

## Appendices

## Appendix:

### Change Point Detection

This section serves as a very brief introduction of change point detection. At first, a hypothesis test is constructed as

Because $r$ is unknown, ${\mathrm{max}}_{r}\text{\hspace{0.17em}}{S}_{r}$ is our test statistic for the hypothesis. The critical value to reject the null hypothesis can be calculated by a $p$-value approximation, once the significance size is determined. The approximation has been obtained by Srivastava and Worsley^{23} as follows:

## Acknowledgments

We thank Academia Sinica (AS-99-TP-AB5) and the National Science Council of Taiwan (NSC 98-2118-M-001-022) for the funding to this work. We are thankful for critical discussions with Chi-Fu Yen and Dr. Yen-Chen Lin for OTIR smFRET instrumentation and data analysis and to Dr. Chin-Yu Chen for helping the smFRET DNA experiments.