title: Maximum likelihood analysis of non-equilibrium solution-based single-molecule FRET data
authors: Marijn de Boer
text: 
  1.5
  
  
  
  Maximum likelihood analysis of non-equilibrium solution-based single-molecule FRET data
      Marijn de Boer 
   Groningen Biomolecular Sciences and Biotechnology Institute 
   University of Groningen 
   Groningen, The Netherlands 
   mailto:marijndeboer4@gmail.commarijndeboer4@gmail.com
      
  ===============================================================================================================================================================================================
  
  
  
  
  
  
  Measuring the Förster resonance energy transfer (FRET) efficiency of freely diffusing single molecules provides information about the sampled conformational states of the molecules. Under equilibrium conditions, the distribution of the conformational states is independent of time, whereas it can vary over time under non-equilibrium conditions. In this work,  we consider the problem of parameter inference on non-equilibrium solution-based single-molecule FRET data. With a non-equilibrium model for the conformational dynamics and a model for the conformation-dependent FRET efficiency distribution, the likelihood function could be constructed. The model parameters, such as the rate constants of the non-equilibrium conformational dynamics model and the average FRET efficiencies of the different conformational states, have been estimated from the data by maximizing the appropriate likelihood function via the Expectation-Maximization algorithm. We illustrate the likelihood method for a few simple non-equilibrium models and validated the method by simulations. The likelihood method could be applied to study protein folding, macromolecular complex formation, protein conformational dynamics and other non-equilibrium processes at the single-molecule level and in solution.
  
  
  
  
  
  
  § INTRODUCTION
  
  Single-molecule methods have emerged as a powerful tool to study biological processes on the molecular level <cit.>. Compared to bulk methods, at which populations of molecules are observed, single-molecule methods allow one to observe individual molecules one after another. The advantage of observing individual molecules is that information about the distribution can be obtained, rather than only providing information about the population average. However, the disadvantage of single-molecule methods is that only a subsample of the population of molecules is seen, and thus, to draw any meaningful conclusion about the entire population, statistically correct inference is needed. In this work, we look at the problem of obtaining parameter estimates from single-molecule Förster resonance energy transfer (smFRET) measurements <cit.>. 
  
  In an smFRET measurement, a molecule is labeled with a donor and an acceptor dye and the energy transfer from the donor to the acceptor is measured <cit.> (Figure 1A). The efficiency of energy transfer, termed the FRET efficiency E, can be determined by measuring the donor and acceptor emission intensities or by measuring the lifetime of the excited state of the donor dye <cit.>. Importantly, E depends on the distance between the donor and acceptor dye r via E =  R_0^6 / (R_0^6 + r^6 ), where R_0 is a constant called the Förster distance <cit.>. E is most sensitive to distance changes around R_0. Commonly used donor and acceptor pairs have R_0 values between 4 and 6 nm. Thus, the sensitivity range of FRET matches with the typical size of biomolecules, making smFRET ideal to use as a 'spectroscopic ruler' <cit.> to study distances between or within proteins <cit.>, nucleic acids <cit.> or other molecules. The distance information can be used to uncover the sampled conformational states of the molecule and the rates of interconversion between these states <cit.>.
  
  
  
  
  In a solution-based smFRET measurement (also termed diffusion-based smFRET or freely-diffusion smFRET) a very low concentration (1- 50 pM) of donor and acceptor labeled molecules is used <cit.> (Figure 1). The low concentration ensures that at most only one molecule is present in the excitation volume, which is created by the focused laser beam of a confocal microscope. During the diffusional transit through the excitation volume, the donor is excited, and the labeled molecule generates a short burst of photons. The detected photons can the used for the determination of E. In an alternating laser excitation (ALEX) <cit.> or pulsed interleaved excitation (PIE) <cit.> scheme, the presence of the acceptor dye is checked by using an additional laser.  Because of the short observation time on each individual molecule (typically around 1 ms), only a limited number of photons can be detected (typically tens or hundreds of photons per burst) leading to large uncertainties in the measured E values. For correct statistical inference, the E values of a large number of molecules, typically hundreds or thousands, are determined in a single measurement. When the measurement is performed under equilibrium conditions (Figure 2A), a histogram for the E values can be made and fitted with a mixture of multiple Gaussian distributions (e.g. <cit.>), for which each Gaussian distribution reflects a conformational state of the molecule. That is, in the mixture model, the mean of the Gaussian distribution represents the donor-acceptor distance of the conformational state and the mixture weight corresponds to its relative abundance in the population. When the measurement is performed under non-equilibrium conditions (Figure 2B), the E distribution could vary over time, and fitting with a time-independent mixture model would, in general, be inappropriate. 
  
  Some examples of non-equilibrium processes that have been studied with solution-based smFRET are: the association and dissociation of the nucleosome core particle (NCP) <cit.>, the dissociation of the ribosome recycling factor ATP-binding cassette E1 (ABCE1) from the small ribosomal subunit <cit.>, promoter escape of the bacterial RNA polymerase (RNAP) <cit.> and the stepping action of molecular motors <cit.>. In these measurements, the system was brought out of equilibrium by diluting the sample into a new buffer environment, in which also the measurement is performed. Subsequently, fluorescence bursts were detected while the molecules diffuse through the excitation volume and they relaxes back to equilibrium.   
  
  In this paper, we provide an analysis framework for the analysis of non-equilibrium solution-based smFRET data. Given a non-equilibrium model for the conformational dynamics and a model for the conformation-dependent FRET efficiency distribution, the likelihood function could be constructed. The likelihood function can be used for various statistical inference problems, including finding the best model parameters that describe the data. The model parameters, such as the rate constants of the non-equilibrium model and the average FRET efficiencies of the different conformational states, were found by maximizing the likelihood function. In this paper, we will provide the formulation of this likelihood function and we show that the parameter inference problem can be solved efficiently with the Expectation-Maximization (EM) algorithm <cit.>.
  
  
  
  
  
  § MODEL DEFINITION
  
  In this section, we define the quantities of our model for a solution-based smFRET measurement. Suppose a measurement starts at time 0 and stops at time T (T>0). During this period a random number of molecules are detected while they diffuse through the excitation volume of a focused laser beam. We denote the number of detected molecules by the random variable N(T). N(T) is a counting process, meaning that N(T) is non-negative, integer valued and non-decreasing in T. For notational convenience we write N(T) simply as N and understand its explicit dependence on T. We label each detected molecule by i (i ∈{1,…,N}) based on its detection time as follows. The times the N molecules are detected is given by the sequence of random variables T_1,…,T_N, where T_i belongs to molecule i and are ordered according to 0 ≤ T_1< ⋯ <T_N ≤ T. We assume throughout the work that the time spend in the excitation volume (typically 1 ms) is short compared to the conformational dynamics of the molecule. Thus, if T_i^e and T_i^l are the time the molecule enters and leaves the excitation volume, respectively, then T_i can be defined anywhere in between T_i^e and T_i^l (e.g. (T_i^l-T_i^e) /2). We also assume that the concentration of labeled molecules is low, so that the probability of having two molecules simultaneously in the excitation volume is negligible.  
  
  Most biomolecules, such as proteins, exist  in a countable number of conformational states that are separated by kinetic barriers <cit.>. Suppose we want to describe a population of identical molecules that can exist in K different conformational states (K ≥ 1), which are numbered from 1 to K. The conformational state of molecule i is a random variable Z_i with a discrete state space 1, … ,K. Due to measurement noise, samples from Z_i are hidden to the observer and are thus the latent variables of our model. The N detected molecules give rise to the sequence Z_1 …,Z_N. We expect that the vast majority of bursts originate from different molecules, so therefore we model the Z_i to be independent. In a non-equilibrium measurement, the probability that a molecule is in a particular conformational state varies over time. Therefore, the distribution of Z_i is conditionally dependent on T_i. This property makes Z_1, … ,Z_N a non-stationary random process. The conditional distribution of Z_i given T_i is denoted by p(z | t; θ)  P(Z_i=z | T_i=t), where z ∈{ 1,…,K } and θ is a vector of model parameters that we want to determine from the data. For example, θ could represent a (set of) conversion rate(s) from one conformational state to another. Note that if the population of molecules is in equilibrium, then the distribution Z_i is independent of time. Then Z_1,…,Z_N is a sequence of independent and identically distributed random variables and is thus a stationary random process. 
  
  As stated above, the conformational states of the molecules cannot be directly observed. However, different conformational states can give rise to different donor-acceptor distances, which can be determined by measuring their FRET efficiency. To avoid problems with parameter identifiability, we assume that each conformational state gives a unique donor-acceptor distance. However, the differences between the distances of the conformational states can be arbitrarily small. The FRET efficiency of molecule i is denoted by the random variable E_i. In principle, the value of E_i lies between 0 and 1. For practical reasons, this assumption is often relaxed and the range of E_i is simply extended to the real line. The N detected molecules give a sequence E_1,…,E_N. The E_i are taken to be independent, because they are expected to originate from different molecules. In a non-equilibrium measurement, the E_i are not identically distributed (see below), so E_1,…,E_N is a non-stationary random process.
  
  We define the distribution of E_i conditional on Z_i by p(e | z; η_z )  P(E_i=e | Z_i=z), where η_z is a vector of parameters that we want to infer from the data. In general, a good approximation for this conformation-dependent FRET efficiency distribution is a Gaussian distribution, which is parameterized by its mean μ_z and standard deviation σ_z. The parameter vector η_z is in this case η_z={μ_z,σ_z }. In the Gaussian approximation, μ_z represents the average FRET efficiency of the zth conformational state and is related to the donor-acceptor distance via μ_z =  R_0^6 / (R_0^6+r^6 ). The σ_z determines the width of the Gaussian distribution and is influenced by the measurement noise (typically σ_z ≈ 0.05-0.10). In principle other distributions could also be considered, such as a beta distribution <cit.>, but then the computation of the maximum likelihood estimate for η_z is more complex. 
  
  We define the distribution of E_i conditional on T_i by p(e|t; Ω)  P(E_i=e|T_i=t), where Ω is a vector of parameters we define below. By using the sum rule of probability, we obtain 
  
      p(e|t; Ω) = ∑_z=1^K P(E_i=e,Z_i=z|T_i=t)
  
  By using the fact that the FRET efficiency depends on the donor-acceptor distance (and R_0), and thus on the conformational state of the molecule, and not directly on time, we get
  
      p(e|t; Ω) =    ∑_z=1^K P(E_i=e|Z_i=z,T_i=t) P(Z_i=z|T_i=t) 
      
      =    ∑_z=1^K P(E_i=e|Z_i=z) P(Z_i=z|T_i=t)
  
  so
  
      p(e|t; Ω) = ∑_z=1^K p(e | z ; η_z) p(z | t ; θ)
  
  From Eq. <ref> it follows that the parameter vector Ω is {θ,η_1,…,η_K}. p(e|t; Ω) gives the probability that a molecule detected at time t has a FRET efficiency of e, which, according to Eq. <ref>, is equal to the sum of probabilities that the FRET efficiency is generated by the different conformational states (at time t). An example of p(e|t; Ω), p(e | z ; η_z) and p(z | t ; θ) is given in Figure 3. We will further elaborate on how to deduce the model parameters from the measurement data and after that, a few examples will be given for illustration purposes.
  
  
  
  
  
  § LIKELIHOOD ESTIMATION
  
  We use the method of maximum likelihood estimation to infer the model parameters from the data. Maximum likelihood estimation is a statistical inference method that has often been applied in the analysis of single-molecule data <cit.>. From the model, the likelihood function can be constructed and maximized with respect to the unknown model parameters. This yields the maximum likelihood estimates for the parameters, which have the desired property that they converge (almost surely) to the true parameters if the amount of data becomes large and if the statistical model is correct. 
  
  The dataset of a solution-based smFRET measurement can be represented by D={(e_1,t_1 ),…,(e_n,t_n)}, where e_i, t_i and n are a sample of E_i, T_i and N, respectively (Figure 1B). The parameters we want to estimate from the dataset are represented by Ω.
  
  By using the fact that the observations on each molecule are independent, the likelihood function can be expressed as the following product
  
      L(Ω ; D)' =     P(E_1=e_1, …,E_n=e_n,T_1=t_1, …,T_n=t_n;Ω) 
      
      =     P(T_1=t_1, …,T_n=t_n) ∏_i=1^n p(e_i | t_i ; Ω )
  
  where P(T_1=t_1, …,T_n=t_n) is the probability of detecting the n molecules at time t_1,…,t_n. Because this probability does not depend on Ω, we can work with the conditional likelihood
  
      L(Ω ; D) = ∏_i=1^n p(e_i | t_i ; Ω ) = ∏_i=1^n∑_z=1^K p(e_i | z ; η_z) p(z | t_i ; θ)
  
  where for the last step we used Eq. <ref>. For computational convenience, we take the natural logarithm of L(Ω;D) 
  
      l(Ω ; D) = log L(Ω ; D)  = ∑_i=1^nlog∑_z=1^K p(e_i | z ; η_z) p(z | t_i ; θ)
  
  The maximum likelihood estimates for θ and η_1,…,η_K are found from maximizing l(Ω ; D)
  
      Ω̂  =  _Ω l(Ω ; D)
  
  where Ω̂ = {θ̂, η̂_z,…,η̂_K }. Unfortunately, the maximization of l(Ω;D) is, in general, quite difficult because of the summation inside the logarithm. We will show that the maximization can often be simplified when the EM algorithm is used <cit.>. The EM algorithm is an iterative method that can simplify the optimization problem when the model involves latent variables (e.g. as in hidden Markov models <cit.>).
  
  The EM algorithm starts with an initial estimate for the parameters Ω^(0) that are iteratively updated via two subsequent steps, termed the Expectation (E) and Maximization (M) step. After k iterations, we use our current estimate of Ω, i.e. Ω^(k)={θ^(k),η_1^(k),…,η_K^(k)}, in the E step to calculate
  
      Q(Ω | Ω^(k)) =    𝔼_Z | E, T, Ω^(k) [l(Ω;Z,E|T)] 
      
      =    ∑_i=1^n∑_z=1^K p(z | e_i , t_i; Ω^(k)) logp(e_i,z | t_i; Ω)
      
      =    ∑_i=1^n∑_z=1^K p(z | e_i , t_i; Ω^(k)) log( p(e_i | z; η_z) p(z|t_i ; θ) )
  
  where l(Ω;Z,E|T) is the log-likelihood of Z={Z_1, …, Z_n} and E={E_1, …, E_n} conditional on T={T_1, …, T_n}, and 𝔼_Z | E,T, Ω^(k)[·] denotes the expectation with respect to the conditional distribution of Z given E, T and the current estimates of the parameters Ω^(k). In Eq. <ref>, p(e_i | z; η_z ) and p(z| t_i ; θ) are the models for the conformation-dependent FRET efficiency and the non-equilibrium conformational dynamics, respectively (Section 2), and p(z | e_i,t_i; Ω^(k)) is calculated with Bayes rule 
  
      p(z | e_i,t_i; Ω^(k)) = p(e_i | z; η_z^(k)) p(z|t_i ; θ^(k))/∑_z=1^K p(e_i | z; η_z^(k)) p(z|t_i ; θ^(k))
  
  In the M step, the parameters are updated by maximizing Q(Ω | Ω^(k)) with respect to Ω
  
      Ω^(k+1)  =  _Ω Q(Ω | Ω^(k))
  
  The E and M step are repeated until convergence is reached. One reasonable convergence test would be to stop iterating when Q(Ω^(k+1) | Ω^(k)) - Q(Ω^(k) | Ω^(k)) < ϵ, where ϵ is a threshold value. 
  
  The optimization problem of Eq. <ref> can be expressed as separate optimization problems, because Q(Ω|Ω^(k)) can be expressed as the following sum
  
      Q(Ω | Ω^(k)) = h(θ | Ω^(k)) + ∑_z=1^K g(η_z | Ω^(k))
  
  where
  
      h(θ | Ω^(k)) = ∑_i=1^n∑_z=1^K p(z | e_i , t_i; Ω^(k)) logp(z|t_i ; θ)
  
  and
  
      g(η_z | Ω^(k)) = ∑_i=1^n p(z | e_i , t_i; Ω^(k)) logp(e_i | z; η_z)
  
  Thus, the parameters that maximize Q(Ω|Ω^(k)) are found from maximizing each of the individual terms in Eq. <ref>
  
      θ^(k+1)  =  _θ h(θ| Ω^(k))
  
  and
  
      η_z^(k+1)  =  _η_z g(η_z| Ω^(k))
  
  for z ∈{1,…, K }.
  
  Up to this point, we did not make any explicit assumption about the distributions p(e_i |z; η_z ) and p(z|t_i; θ), and is thus applicable to any model for the conformation-dependent FRET efficiency and the non-equilibrium conformational dynamics. It shows how the maximum likelihood estimates can be found with the EM algorithm, and that the optimization problem can be expressed as K+1 smaller optimization tasks involving the separate optimization for θ and the K η_z parameter vectors. This result is important from a computational perspective. Without the EM algorithm, the maximization of the likelihood function of Eq. <ref> involves K+1  parameter vectors. However, by using the EM algorithm, the maximum is found by solving K+1 separate optimization problems, which, in general, drastically simplifies the problem.
  
  
  
  § GAUSSIAN FRET EFFICIENCY DISTRIBUTION
  
  Before we look at some illustrative examples of the above analysis framework, we first consider how the M step can be implemented in case the conformation-dependent FRET efficiency distribution is Gaussian. In this case, p(e_i|z;μ_z, σ_z) = 𝒩(e_i ; μ_z, σ_z), where
  
      𝒩(e_i ; μ_z, σ_z) = 1/√(2 πσ_z^2) e^-(e_i - μ_z)^2/2 σ_z^2
  
  The Gaussian distribution is parameterized by its mean μ_z and standard deviation σ_z, so η_z={μ_z,σ_z}. As long as the FRET efficiency does not lie close to 0 or 1, the Gaussian approximation is often a very good approximation of the FRET efficiency distribution (of a single distance). 
  
  In the M step, we need to maximize g(η_z| Ω^(k)) = g(μ_z, σ_z| Ω^(k)) with respect to μ_z and σ_z. From Eq. <ref> we have
  
      g(μ_z, σ_z | Ω^(k))  =    ∑_i=1^n p(z | e_i , t_i; Ω^(k)) log𝒩(e_i ; μ_z, σ_z)
      
      =     - ∑_i=1^n p(z | e_i , t_i; Ω^(k)) ( (e_i - μ_z)^2/2 σ_z^2 + log√(2 πσ_z^2))
  
  In Eq. <ref>, p(z|e_i,t_i;Ω^(k))  is determined during the E step via Eq. <ref> and will be considered in more detail in the next section. Maximizing Eq. <ref> with respect to μ_z and σ_z gives
  
      μ_z^(k+1) = ∑_i=1^n p(z|e_i,t_i;Ω^(k)) e_i/∑_i=1^n p(z|e_i,t_i;Ω^(k))
  
  and 
  
      σ_z^(k+1) = √(∑_i=1^n p(z|e_i,t_i;Ω^(k)) ( e_i - μ_z^(k+1))^2/∑_i=1^n p(z|e_i,t_i;Ω^(k)))
  
  In conclusion, the parameters μ_z and σ_z are updated during the M step by determining the weighted sample average and the weighted sample standard deviation of the measured FRET efficiencies. 
  
  
  
  
  § ILLUSTRATIVE EXAMPLES
  
  In this section we consider simple examples of models for non-equilibrium conformational dynamics, with a focus on the calculation of the E and M step. Throughout this section we assume that the conformation-dependent FRET efficiencies are Gaussian, and its implementation for the M step was covered in the previous section. 
  
  
  
  
   §.§ Example 1: Exponential two-state model
  
  Suppose a molecule can be in two conformations, denoted by state 1 and state 2. At the beginning of the measurement, all the molecules start  in state 1 and go randomly to state 2 with a rate constant λ. If the rate constant is independent of time, then the lifetime of state 1 is exponentially distributed. Exponential lifetimes are commonly observed in smFRET data <cit.>. We assume that molecules in state 2 cannot convert back to state 1. This process could, for example, correspond to the folding of a protein molecule. In this process, all molecules are initially in the unfolded state (state 1) and fold with a rate constant λ to the folded state (state 2). In this particular example, the process is irreversible, so the folded molecule cannot convert back to the unfolded state. Under these conditions, the model for the conformational dynamics is 
  
      p(z | t; λ) = {[     e^-λ t      z = 1; 1 - e^-λ t      z = 2 ].
  
  with λ>0. In this example the K+1 vectors (where K = 2) are θ={λ}, η_1={μ_1,σ_1 } and η_2={μ_2,σ_2 }. And thus the set of model parameters is Ω={λ,μ_1,σ_1,μ_2,σ_2 }. In Figure 4A, p(z | t; λ) is plotted for λ =0.01 s^-1, and in Figure 4B, p(e | t; Ω) is plotted by using Eq. <ref>, where p(z | t; η_z) is a Gaussian distribution with μ_1=0.5, σ_1=0.07, μ_2=0.7 and σ_2=0.06.
  
  We can use the EM algorithm to find the maximum likelihood estimates for Ω. For the E step calculation in Eq. <ref>, we use Eq. <ref> constructed with our current parameter estimates λ^(k), μ_1^(k), σ_1^(k), μ_2^(k) and σ_2^(k) to calculate 
  
      p(z | e_i, t_i; Ω^(k)) = {[        e^-λ^(k) t_i𝒩(e_i ; μ_1^(k), σ_1^(k)) / e^-λ^(k) t_i 𝒩 (e_i ; μ_1^(k), σ_1^(k)) + ( 1 - e^-λ^(k) t_i) 𝒩(e_i ; μ_2^(k), σ_2^(k))                                                                                                                                  z = 1; ( 1 - e^-λ^(k) t_i) 𝒩(e_i ; μ_2^(k), σ_2^(k)) / e^-λ^(k) t_i 𝒩(e_i ; μ_1^(k), σ_1^(k)) + ( 1 - e^-λ^(k) t_i) 𝒩(e_i ; μ_2^(k), σ_2^(k))                                                                                                                                  z = 2 ].
  
  where 𝒩(e_i ; μ_1^(k), σ_1^(k)) and 𝒩(e_i ; μ_2^(k), σ_2^(k)) are defined by Eq. <ref>. With p(z|e_i,t_i; Ω^(k) ) calculated, we proceed to the M step and update the parameters. The M step for μ_1, σ_1, μ_2 and σ_2 are given by Eq. <ref> and <ref> and using p(z|e_i,t_i;Ω^(k) ) as given by Eq. <ref>. To update λ, we use Eq. <ref>, <ref> and <ref> to obtain
  
      h(λ | Ω^(k))  =    ∑_i=1^n∑_z=1^K p(z | e_i , t_i; Ω^(k)) logp(z|t_i ; λ)
      
      =    ∑_i=1^n  - p(z = 1| e_i , t_i ; Ω^(k)) λ t_i + p(z = 2| e_i , t_i ; Ω^(k)) log( 1  - e^-λ t_i)
  
  To give the next iteration value λ^(k+1), Eq. <ref> needs to be maximized with respect to λ. No closed form expression for λ^(k+1) exists. Fortunately, h(λ | Ω^(k)) is a concave function of λ, because
  
      ∂^2 h(λ | Ω^(k))/∂^2 λ = - ∑_i=1^np(z = 2| e_i , t_i ; Ω^(k)) t_i^2 e^λ t_i/( e^λ t_i -1 )^2 < 0
  
  for every λ>0, so the maximum of h(λ | Ω^(k)) can be found relatively straightforwardly with any numerical optimization algorithm. 
  
  Note that in this particular example, the multi-parameter optimization problem of Eq. <ref> involves 5 parameters and is reduced by the EM algorithm to 3 separate and simpler optimization problems. 
  
  
  
  
  
   §.§ Example 2: Exponential two-state model with offset
  
  We modify example 1, by assuming that at the beginning of the measurement some molecules within the population are already in state 2. Let π denote the fraction of molecules being in state 1 at the beginning of the measurement and 1-π the fraction of molecules in state 2. Then the model for the conformational dynamics is 
  
      p(z | t; λ,π) = {[     π e^-λ t        z = 1; 1 - π e^-λ t        z = 2 ].
  
  with λ > 0 and 0<π<1.  In this example, θ={λ,π}, η _z={μ_z,σ} and Ω={λ,π,μ_1,σ_1,μ_2,σ_2 }.  In Figure 4C, p(z | t; λ) is plotted for λ =0.01 s^-1 and π=0.8.  In Figure 4D, p(e | t; Ω) is plotted by using Eq. <ref>, where p(z | t; η_z) is a Gaussian distribution with μ_1=0.5, σ_1=0.07, μ_2=0.7 and σ_2=0.06.
  
  We use the EM algorithm to find the maximum likelihood estimate for Ω. During the E step, we use Eq. <ref> and <ref> together with our current parameter estimates λ^(k), π^(k), μ_1^(k), σ_1^(k), μ_2^(k) and σ_2^(k), to calculate
  
      p(z | e_i, t_i; Ω^(k)) = {[        π^(k) e^-λ^(k) t_i𝒩(e_i ; μ_1^(k), σ_1^(k)) /π^(k) e^-λ^(k) t_i 𝒩(e_i ; μ_1^(k), σ_1^(k)) + ( 1 - π^(k) e^-λ^(k) t_i) 𝒩(e_i ; μ_2^(k), σ_2^(k))                                                                                                                                                  z = 1; ( 1 - π^(k) e^-λ^(k) t_i) 𝒩(e_i ; μ_2^(k), σ_2^(k)) /π^(k) e^-λ^(k) t_i 𝒩(e_i ; μ_1^(k), σ_1^(k)) + ( 1 - π^(k)e^-λ^(k) t_i) 𝒩(e_i ; μ_2^(k), σ_2^(k))                                                                                                                                                  z = 2 ].
  
  The M step for μ_1 ,σ_1, μ_2 and σ_2 are given by Eq. <ref> and <ref> and using p(z|e_i,t_i;Ω^(k) ) as given by Eq. <ref>. To update λ and π, we use Eq. <ref>, <ref> and <ref> to calculate
  
      h(λ, π | Ω^(k))  =    ∑_i=1^n∑_z=1^K p(z | e_i , t_i; Ω^(k)) logp(z|t_i ; λ,π)
       
      =    ∑_i=1^n p(z = 1| e_i , t_i ; Ω^(k)) ( logπ - λ t_i  )  + p(z = 2| e_i , t_i ; Ω^(k)) log( 1  - π e^-λ t_i)
  
  which can be maximized numerically to give λ^(k+1) and π^(k+1). Note that the complex multi-parameter optimization problem of Eq. <ref> is reduced by the EM algorithm to 3 separate and simpler optimization problems.
  
  
  
   §.§ Example 3: Non-exponential two-state model
  
  Similar as in the previous two examples, we consider here that a molecule can be in two conformational states, denoted by state 1 and state 2. At the beginning of the measurement all molecules are in state 1, and transit randomly to state 2. The process is irreversible, so once a molecule is in state 2 it cannot convert back to state 1. However, in contrast to the previously discussed examples, instead of the exponential lifetime of state 1, we consider here the following simple model for non-equilibrium conformational dynamics
  
      p(z | t; τ) = {[ τ/τ +t  z = 1; t/τ +t  z = 2 ].
  
  with τ>0. In this model, θ={τ}, η_1={μ_1,σ_1 }, η_2={μ_2,σ_2 } so Ω={τ,μ_1,σ_1,μ_2,σ_2 }. The model can be interpreted as the irreversible transition from state 1 to state 2, occurring with a (hazard) rate that depends on time as -∂/∂ tlog( p(z = 1 | t; τ) ) = 1/(τ + t). Thus, the transition becomes slower as time increases, whereas in example 1 and 2 it was independent of time. In Figure 4E, p(z | t; τ) is plotted for τ =100 s, and in Figure 4F, p(e | t; Ω) is plotted by using Eq. <ref>, with p(z | t; η_z) being Gaussian with μ_1=0.5, σ_1=0.07, μ_2=0.7 and σ_2=0.06.
  
  Irrespective of the precise molecular interpretation, we can use the EM algorithm to find the maximum likelihood estimate for Ω. In the E step, we use our current parameter estimates τ^(k), μ_1^(k), σ_1^(k), μ_2^(k) and σ_2^(k) together with Eq. <ref> and <ref> to calculate 
  
      p(z | e_i, t_i; Ω^(k)) = {[ τ^(k)𝒩(e_i ; μ_1^(k), σ_1^(k)) /τ^(k)𝒩(e_i ; μ_1^(k), σ_1^(k)) + t_i 𝒩(e_i ; μ_2^(k), σ_2^(k))                                                                                          z = 1;  t_i 𝒩(e_i ; μ_2^(k), σ_2^(k)) /τ^(k)𝒩(e_i ; μ_1^(k), σ_1^(k)) + t_i 𝒩(e_i ; μ_2^(k), σ_2^(k))                                                                                          z = 2 ].
  
  The M step for μ_1, σ_1, μ_2 and σ_2 are given by Eq. <ref> and <ref> and using p(z | e_i, t_i; Ω^(k)) as given by Eq. <ref>. To update τ, we combine Eq. <ref>, <ref> and <ref> to give
  
      h(τ | Ω^(k))  =    ∑_i=1^n∑_z=1^K p(z | e_i , t_i; Ω^(k)) logp(z|t_i ; τ)
       
      =    ∑_i=1^n p(z = 1| e_i , t_i ; Ω^(k)) log( τ/τ + t_i)  + p(z = 2| e_i , t_i ; Ω^(k)) log( t_i/τ + t_i)
      
      =    ∑_i=1^n p(z = 1| e_i , t_i ; Ω^(k)) logτ - log(τ + t_i) + p(z = 2| e_i , t_i ; Ω^(k)) log t_i
  
  Maximizing h(τ | Ω^(k)) with respect to τ gives 
  
      τ^(k+1)  =  _τ( ∑_i=1^n p(z = 1| e_i , t_i ; Ω^(k)) logτ - log(τ + t_i) )
  
  which can be solved relatively straightforwardly with any numerical method. 
  
  
  
  
   §.§ Example 4: K-state model at equilibrium
  
  If a population of molecules are in equilibrium, then the stochastic process is stationary. Then, the probabilities to be in one of the conformational states do not depend on time. Let in this case p(z | t_i ; θ )=p(z;θ )=w_z. 
  Suppose a population of molecules are in equilibrium and can acquire K conformational states. In this example θ={w_1, … ,w_K }, η_z={μ_z, σ_z} and Ω={w_1, …,w_K,μ_1,…,μ_K,σ_1,…,σ_K }. We use the EM algorithm to find the maximum likelihood estimate for Ω. During the E step, we use Eq. <ref> together with our current parameter estimates Ω^(k), to calculate p(z | e_i, t_i ; Ω^(k)) = p(z | e_i ; Ω^(k)), where
  
      p(z | e_i ; Ω^(k)) = w_z^(k)𝒩(e_i; μ_z^(k), σ_z^(k))/∑_z=1^K w_z^(k)𝒩(e_i; μ_z^(k), σ_z^(k))
  
  The M step for μ_z and σ_z are given by Eq. <ref> and <ref> and using p(z | e_i ; Ω^(k)) as given by Eq. <ref>. To update w_z, we combine Eq. <ref> and <ref>
  
      h(w_1, …,w_K| Ω^(k))  =    ∑_i=1^n∑_z=1^K p(z | e_i ; Ω^(k)) logw_z
      
      =    ∑_z=1^K[ ∑_i=1^n p(z | e_i ; Ω^(k)) ] log w_z
  
  This needs to be maximized with respect to w_1, … ,w_K subject to the constraint ∑_z=1^K w_z = 1. This optimization problem is equivalent to finding the maximum likelihood estimates of a multinomial distribution, so  
  
      w_z^(k+1) = 1/n∑_i=1^n p(z | e_i ; Ω^(k))
  
  It may not come as a surprise that, when the process is in equilibrium and the conformation-dependent FRET efficiency distributions are Gaussian, we retrieve the E and M step of a Gaussian Mixture Model (GMM). 
  
  
  
  § MODEL SELECTION
  
  An important part of statistical inference through maximum likelihood is the selection of the appropriate model for the data. Here, we propose a simple scheme to visualize how well the model with the estimated parameters fit the data. The idea is to compare a histogram constructed over a certain time interval with the expected distribution of the model over that interval. Other measures could also be used for model selection, such as the Bayesian information criterion (BIC) or the Akaike information criterion (AIC), but those will not be considered here.
  
  To make our idea concrete, we first partition the time period [0,T] into m non-overlapping time intervals, which we denote by 𝒯_j, j ∈{1,…, m}.  For j =1, the time interval 𝒯_1 is defined as 𝒯_1=[0,τ_1^1] with τ_1^1>0. For j ∈{2,…, m-1}, the intervals are defined as 𝒯_j=(τ_j^0,τ_j^1] where τ_j-1^1=τ_j^0<τ_j^1=τ_j+1^0. For j=m, 𝒯_m=(τ_m^0,T] with τ_m^0<T. Next, we construct a histogram for the FRET efficiency of all molecules that have been detected within the time interval T_j. In Figure 2 an example is shown for m = 2 and the corresponding histograms over 𝒯_1=[0 s,400 s] and 𝒯_2=(400 s,800 s]. We can compare this histogram with the FRET efficiency distribution over the time interval 𝒯_j, which we denote by p_𝒯_j(E)  and is equal to
  
      p_𝒯_j(E) = 1/N_j∑_t_i ∈𝒯_j p(E | t_i ; Ω̂ )
  
  where E is the FRET efficiency and p(E | t_i ; Ω̂ ) is given by Eq. <ref>. The normalization constant N_j is equal to the number of detected molecules in 𝒯_j and have the property that ∑_j=1^m N_j = n. The time intervals should be chosen such that N_j ≥ 1 for every interval. By using Eq. <ref>, Eq. <ref> can be expressed as 
  
      p_𝒯_j(E) = 1/N_j∑_t_i ∈𝒯_j∑_z=1^K p(E | z ; η̂_z) p(z | t_i ; θ̂)
  
  were η̂_z and θ̂ are the maximum likelihood estimates of η_z and θ, respectively. Eq. <ref> depends on the selected model for the conformation-dependent FRET efficiency distribution and the selected model for the conformational dynamics, which are both parameterized by the maximum likelihood estimates as obtained from the complete dataset. The quality of the fit can now be visualized by examining the discrepancy between the histograms and the model prediction of Eq. <ref> (see next section for an example). In principle, the discrepancy could be made precise by any measure of goodness of fit, such as the chi-squared test. 
  
  
  
  § NUMERICAL CALCULATION
  
  In this section, we test the maximum likelihood method on simulated data. To do this, we simulated data with typical experimental parameters and compared these parameters to the maximum likelihood estimates that are obtained from the simulated data. 
  
  In all our simulations, we simulated the data of 3000 bursts. The timepoints the bursts are detected is modeled to be Poisson process with a rate 10 s^-1. In this case, the average time between two subsequent bursts is 100 ms and the average measurement time is 5 min (3000 · 100 ms). The typical residence time of a molecule in the excitation volume is 1 ms, so the probability to have two molecules simultaneously in the excitation volume is 10^-2· 10^-2 = 10^-4 and can thus be ignored.
  
  First, we simulated the FRET efficiency data according to example 1 of Section 5.1. In this irreversible two-state model the lifetime of state 1 is exponentially distributed with an average lifetime of λ^-1 and the conformation-dependent FRET efficiency distributions are a Gaussian distribution for both states. In the simulation, we chose a lifetime of 100 s, so λ = 0.01 s^-1. The average FRET efficiency of state 1 and 2 were set to 0.5 and 0.7, respectively, and the standard deviation to 0.07 and 0.06 for state 1 and 2, respectively. This model is plotted in Figure 4A-B. In Figure 5A, the simulated dataset is shown. To extract the parameters from the simulated data, we maximized the likelihood function by using the EM algorithm as described in Section 5.1. The initial parameters were set to λ^(0) = 0.025 s^-1, μ_1^(0)=0.35, μ_2^(0)=0.9 and σ_1^(0)=σ_2^(0)=0.12 and were updated until Q(Ω^(k+1)|Ω^(k))-Q(Ω^(k)|Ω^(k)) was smaller than 10^-5. In Figure 5B-D, the parameters values during each iteration of the EM algorithm are depicted and shows convergence to values close to the true parameters. To compare the model fit with the data, we constructed 3 FRET efficiency histograms over a time interval of 100 s (Figure 5E). The histograms agree well with the model prediction of Eq. <ref>.
  
  
  To evaluate in more detail the performance of the maximum likelihood method, we simulated 100 datasets, each containing data of 3000 bursts. The average and standard deviation of the estimated parameters are shown in Table 1, which shows that the maximum likelihood estimates indeed converge to the true parameters with a typical error (standard error of the mean) of <1%. 
  
  In Table 2 and 3, the simulation results are shown for alternative models of the conformational dynamics. The data of Table 2, is based on Eq. <ref> with λ = 0.01 s^-1 and π = 0.8 and Table 3 is based on Eq. <ref> with τ = 50 s. The performance of the maximum likelihood method was evaluated by simulating 100 datasets, each consisting of 3000 bursts. The EM algorithm, as described in Section 5.2 and 5.3, was used to estimate the model parameters from the simulated data. Table 2 and 3 show that the maximum likelihood estimates are in perfect agreement to the true parameters. Thus, the simulations indicate that the maximum likelihood method can be used to efficiently obtain the model parameters from the data.
  
  
  
  
  
  
  
  
  
  § CONCLUDING REMARKS 
  
  In this paper, we present a framework for the analysis of non-equilibrium solution-based smFRET data. The framework enables obtaining parameter estimates for the conformation-dependent FRET efficiency distribution as well as the distribution that describes the non-equilibrium conformational dynamics. Parameter estimates were found by maximizing the appropriate likelihood function. By making use of the EM algorithm, the maximization task is significantly simplified. We provided illustrative examples of models for the conformational dynamics, which represent various simple non-equilibrium processes. More complex models could be implemented in a similar way.
  
  Information about the conformational dynamics of single-molecules has predominately been obtained through the analysis of immobilized molecules. In these surface-based smFRET measurement, individual molecules are attached to the surface of a glass slide and followed over time by using confocal scanning microscopy <cit.> or total internal reflection microscopy <cit.>. This data could give direct insight into the order of events and the kinetics associated with each step. However, one disadvantage of surface-based measurements is that the molecules need to be specifically immobilized without effecting their function, which can sometimes be challenging, especially for large macromolecular complexes. Immobilization might also be inappropriate if processes like protein folding and chaperone interactions are studied. Besides the technological developments in microfluidic mixing to obtain dynamic information from freely-diffusing molecules with smFRET <cit.>, information about the dynamics can also be obtained with conventional solution-based smFRET measurements as considered in this work <cit.>. Analysis methods for solution-based smFRET measurements have been developed to extract information about fast conformational dynamics, which would occur on timescales similar to the residence time of the molecule inside the excitation volume (around 1 ms) <cit.>. The analysis framework provided in this paper would be appropriate to study slower processes, which would occur on or above the sub-minute timescale. This limit might be significantly lowered due to the recent development of high-throughput smFRET systems <cit.>. 
  
  We note that alternative analysis methods can also be used. For instance, a histogram for the FRET efficiencies could be constructed for different non-overlapping time intervals (see also Section 6). The histograms could then be fitted with a time-independent GMM (Section 5.4), similarly to what would be done in an equilibrium situation. One could equate the estimates for the mixture weights to the probability of acquiring the corresponding conformational states at a timepoint that is equal to the center of the time interval. However, this approximation would only be valid if the distribution to acquire the different conformational states is approximately constant over the time interval. If this approximation can be made, a course-grained picture of the conformational dynamics can be obtained, which can be further analyzed. However, because a sufficient number of observations are needed to properly estimate the mixture weights, a minimum duration for each time interval would be around 3 min for simple systems. Therefore, this time interval method would only be appropriate to analyse very slow processes that occur on time scales much slower than 3 min (say 20 min or above).
  
  This work focuses on obtaining parameter estimates from the data by maximizing the likelihood function. However, further work is needed to establish the errors on those estimates. In principle, multiple sources of error or uncertainty should be considered. First, if the model for the data is incorrect, the wrong likelihood function is maximized. This is the subject of model selection. It would be good practice to test different (nested) models for the data and select the most appropriate model by, for example, the likelihood ratio test. Alternatively, the BIC or the closely related AIC could be used for model selection. Moreover, the discrepancy between the model prediction and the data could be examined as described in Section 6. Secondly, errors can arise because of the limited amount of data. A way to deal with this type of uncertainty is by bootstrapping or by evaluating the Fisher Information. Thirdly, the EM algorithm might converge to a local maximum of the likelihood function instead of the global maximum. Varying the starting conditions of the algorithm might help in finding the global maximum. 
  
  Maximum likelihood estimators have the desired invariance property, i.e., if θ̂ is the maximum likelihood estimator for θ, then g(θ̂) is the maximum likelihood estimator for g(θ). This property can be used to estimate other physical properties from the data. For instance, once the maximum likelihood estimates for the average FRET efficiencies are obtained, i.e. μ_z in case of a Gaussian distribution, then the maximum likelihood estimates for the donor-acceptor distance r could obtained by solving μ_z=  R_0^6 / (R_0^6+r^6 ). This could give important insight in to the molecular structure of the conformational states, especially if multiple measurements are combined with different donor and acceptor labeling positions. 
  
  
  
  § ACKNOWLEDGEMENTS
  
  This work was done in the Poolman lab at the University of Groningen. I would like to thank Monique Wiertsema and Bert Poolman for critically reading of the manuscript. I also would like to thank Monique Wiertsema for the help with the simulations.
  
  
  ieeetr
  
  
  