- Email: [email protected]
Error estimation in the phase plane method of multi-exponential decay analysis
1 June 1998
Optics Communications 151 Ž1998. 313–320
Full length article
Error estimation in the phase plane method of multi-exponential decay analysis )
Eugene G. Novikov
Systems Analysis Department, Byelorussian State UniÕersity, F. Skarina AÕenue 4, Minsk 220050, Belarus Received 20 October 1997; revised 20 January 1998; accepted 22 January 1998
Abstract The problems of error estimation for the parameters of a multi-exponential fluorescence decay law are considered. Two algorithms for the determination of the covariance matrix in the phase plane method, based on the analytical transformations and bootstrap approach, are suggested. Simulation experiments have been performed in order to show the applicability of the proposed algorithms for the determination of the fluorescence decay parameters and their standard deviations. q 1998 Published by Elsevier Science B.V. All rights reserved. Keywords: Fluorescence decay fitting; Multi-exponential decay law; Phase plane method; Prony method; Marquardt nonlinear least squares; Covariance matrix; Variance; Bootstrap
1. Introduction The development history of software algorithms for fluorescence decay analysis is fairly rich w1–4x, but many different aspects of this problem are still of particular interest w5,6x, in particular to develop more accurate, more reliable and faster fitting methods. In many applications fluorescence decay can be adequately approximated by a sum of exponentials w7x: n
IŽt. s
Ý ak exp Žytrt k . ,
probe, excited by a train of short light pulses, is detected as discrete numbers of photons, collected in a finite number of the time channels of a multichannel analyzer. The finite width of an excitation pulse and the delays in a detector make the true decay curve I Ž t . Žimpulse response of a fluorescent probe. convoluted with the instrumental function g Ž t . Žwhich itself is a convolution of an excitation pulse and detector response.. The detected fluorescence decay f Ž t . is then determined by a convolution integral:
Ž1.
ks 1
where n is the number of exponentials, a s a k 4ks1, . . . , n and t s t k 4ks1, . . . , n are the amplitudes and decay constants which have to be defined. The methods of time-resolved fluorescence spectroscopy w4,7x are widely used to obtain information about the time behaviour of fluorescence. One of the more commonly used techniques is time-correlated single photon counting w1x, where the fluorescence response of a
)
E-mail: [email protected]
f Žt. s
t
H0 I Ž t y z . g Ž z . d z ,
g Ž z . s 0,
z F 0.
Ž2.
A variety of methods for fluorescence data fitting has been exercised with different level of success w8,9x. Among them two groups are clearly prominent: iterative nonlinear least squares ŽNLLS. methods, directly optimizing some final criteria Žusually x 2 . Žsimplex method w2,10x, Marquardt method w11x, etc.., and non-iterative methods, based on the transformation of the detected decay to a linear algebraic set of equations with respect to unknown parameters ŽLaplace transforms method w12x, method of moments
0030-4018r98r$19.00 q 1998 Published by Elsevier Science B.V. All rights reserved. PII S 0 0 3 0 - 4 0 1 8 Ž 9 8 . 0 0 0 5 0 - 9
314
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
w13x, method of modulating functions w14x, phase plane method w2,15,16x, Prony’s method w5,17,18x, etc... The methods in the first group, granting high quality of fitting, are relatively slow and need initial guesses for the parameters to be known. When the initial guesses are far from the true values, the convergence may be very slow and time consuming. The non-iterative methods are not so accurate in fitting but extremely fast and normally do not need initial guesses. However, analysis of the existing techniques w8,9x reveals strict influence of the extra factors on the accuracy of fitting Žparameter of transformation for Laplace transform method, cut-off ratio for method of moments, kind of modulation function for method of modulating functions.. Our experience shows that the phase plane method is more stable in that respect, providing acceptable accuracy and high speed of processing. The common way of fluorescence decay analysis may be composed of non-iterative routine, generating initial guesses, and iterative procedure, starting the search with those initial values. The phase plane method is a good choice for the first step and Marquardt NLLS for the second one. One note should be made, concerning the principal sameness of the Prony’s method and phase plane method. Both of them yield the equation of a linear regression with respect to the unknown coefficients, which can be fit by the methods of linear least-squares minimization. The difference is in the way of transformation of the convolution integral Ž2. with the multi-exponential fluorescence response Ž1. into the equation of a linear regression. Prony’s method uses the technique of Z transformations w17,18x, and the phase plane method is based on the integration of the convolution Ž2. with the running upper limit w15,16x. A significant problem, that always accompanies the problem of point estimation, involves determination of the variance and confidence interval for the estimated parameters. This problem is prevalent for most methods mentioned in this paper. Comparably extensive statistical theory has been developed only for the NLLS methods w10,11,19x. A way of estimating the amplitudes and decay time variances has been proposed also for the method of moments w13x. But cut-off problems, which make this method very unstable, especially for multi-exponential decay fitting, become even more pronounced for the calculation of the variances. In the present paper two algorithms of error estimation for the phase plane method are proposed and compared to the estimators, provided by the Marquardt NLLS method. The further considerations are also appropriate for the Prony’s method.
equation Ž2. with the multi-exponential fluorescence response Ž1. to a linear equation with constant coefficients: n
n
Ý c k Fk Ž x . s Ý rk Gk Ž x . , ks 0
Ž3.
ks1
where F0 Ž x . s f Ž x ., x
ky1
Fk Ž x . s
H0 f Ž t .Ž x y t .
s
H0 g Ž t .Ž x y t .
x
r Ž k y 1 . ! d t ,Gk Ž x .
ky1
r Ž k y 1 . ! d t ,k s 1, . . . ,n.
Ž4. This transformation can also be applied to the analysis of multi-exponential decay without convolution Ž g Ž t . is the delta function.. In this case functions G k Ž x . take the form: G k Ž x . s x ky1r Ž k y 1 . !,
k s 1, . . . ,n.
Ž5.
Eq. Ž3. is a linear functional relationship with respect to the coefficients c k , k s 0, . . . ,n and rk , k s 1, . . . ,n, which must be fitted. The decay times t k can be obtained as roots of the following polynomial: n
Ý Žy1. k c kt kyn s 0,
Ž6.
ks 0
and the amplitudes a k are evaluated from the set of linear equations: n
rk s
k
Ý
Ý Žy1. kyj c jy1tmjyk ,
am
ms 1
Ž7.
js1
where t k , k s 1, . . . ,n are known from Eq. Ž6.. The analytical solution of the set Ž7. gives: ny 1
am s
Ý Žy1. k rkq1tmk ks0
5½
n
Ł Ž 1 y tmrt j .
js1, j/m
m s 1, . . . ,n.
5
,
Ž8.
As already mentioned, the curves f Ž x . and g Ž x . are detected in the finite number N of time channels: f ) Ž x i . and g ) Ž x i ., i s 1, . . . , N. Introducing the following matrices: F ) s Fk) Ž x i . 4 ks0 , . . . , n ; is1 , . . . , N ; G ) s G k) Ž x i . 4 ks1 , . . . , n ; is1 , . . . , N
Ž9.
and vectors c s c k 4 ks0 , . . . , n ;r s r k 4 ks1 , . . . , n ,
Ž 10.
Eq. Ž3. can be rewritten as: 2. Phase plane method In this section the phase plane method is shortly outlined. It is based on the transformation of the convolution
c F ) s rG ) .
Ž 11 .
Taking into account statistical distortions, naturally presented in the measured data, one has: e s c F ) I rG ) ,
Ž 12.
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
where e is the vector of residuals, which is, according to the central limit theorem w20x, Gaussian with mean 0 and covariance function: cov Ž e . s ² eeX : s V.
Ž 13.
The coefficients of linear relationship Ž3. c s c k 4ks0, . . . , n and r s rk 4ks1, . . . , n can be obtained by the maximum likelihood method w5,21x, maximizing the following probability density function: P Ž e . A exp Ž yeVy1 eXr2 . ,
Ž 14.
or, similarly, minimizing the quadratic form w22x:
eigenvalue and corresponding eigenvector u of the symn metric matrix R, when Ý2ks0 u 2k s 1 w21x. In our case it is certainly not correct, as the elements of F and G are highly correlated, and their variances are considerably different. On the contrary, since integrating Eq. Ž4. normally smoothes statistical distortions in the data, the relative variance of the elements Fk) and G k) , k s 1, . . . ,n of the vectors F and G is much smaller than the relative variance of the zero element F0) of the vector F. Therefore, we may put c 0 s 1 in Eq. Ž3., thus, getting the equation of a linear regression: n
X
eV I1 eX s Ž c F ) I rG ) . Vy1 Ž c F ) I rG ) .
f )Ž x. s
s u R u s min,
Ž 15.
X
where u s c, r 4 is the vector-string, consisting of two sub-strings, and Rs
F ) V I 1 F )X F ) V I 1 G )X
F ) V I 1 G )X G ) V I 1 G )X
Ž 16.
is the matrix, formed by four sub-matrices of lower dimensions; the prime denotes the transpose operation. To avoid trivial solution Ž c k s 0 and r k s 0., one may set, for n example, Ý2ks0 u 2k s Ý nks0 c k2 q Ý nks1 r k2 s 1, or u 0 s c 0 s 1. Differentiating Eq. Ž15. with respect to the elements of u results in a set of equations, that can be solved numerically. The solution gives the maximum likelihood estimators of the coefficients c and r in Eq. Ž3.. According to the invariance principle w23x, the parameters a and t , calculated from the vectors c and r via Eqs. Ž6. and Ž8. respectively, are estimators, providing the maximum to the likelihood function Ž14., as well. The problem of the maximum likelihood estimation is mainly related to the complex structure of the matrix V, and, particularly, to the dependence of the elements of V on the unknown coefficients c and r. Iterative algorithms w24x, or introducing a reasonable approximation for V w5x are used in order to solve this problem. Both approaches require to calculate the inverse matrix Vy1, the dimension of which corresponds to the number of channels N and may be very large Žhundreds or thousands of channels.. This gives rise to problems of numerical implementation, resulting in loss of accuracy. Moreover, the first approach may cause obstacles, naturally accompanying the iterative algorithms, like arriving into local minima and increased time of calculation. In the second one, proper approximation for the matrix V is not always obvious. Two special cases, when minimization of Eq. Ž15. is straightforward, should be mentioned at this point. If all elements of the vectors F and G are independent of each other and identically distributed with the same variance, then the problem of minimization of the quadratic form Ž15. is equivalent to the determination of the minimal
n
Ý c k Fk) Ž x . y Ý rkGk) Ž x . , ks 1
X
315
Ž 17.
ks1
where functions Fk) Ž x . and G k) Ž x ., k s 1, . . . ,n are assumed to be noiseless compared to F0) , and, consequently, statistically independent of each other. In this case, maximum likelihood estimators of the coefficients c k and rk , k s 1, . . . ,n are approximately equivalent to the estimators, found by the linear least squares method w21,22x, directly applied to the linear regression equation Ž17.. This approximation allows for easy analytical calculation of the inverse matrix of V, which is diagonal in that case: y1 y1 ) Vy1 s Viy1 j 4 ;Vi j s 0,i / j;Vi j s 1rf Ž x i . ,i s j.
Ž 18. Testing, undertaken in the Section 6, shows that this approximation ensures extremely fast and reasonably accurate fittings. Hence, further we shall assume that c 0 s 1, and vector c has the form: c s c k 4ks1, . . . , n .
3. Analytical estimation of the covariance matrix The covariance matrix of the maximum likelihood estimators takes the form w22,23x: Ry1 s
K c ,c 4 K r ,c 4
K c ,r 4 , K r ,r 4
Ž 19.
where K c,c4 and K r,r 4 are the auto-covariance matrices of the coefficients c and r respectively, and K c,r 4 and K r,c4 are the cross-covariance matrices of c and r. In this section we consider a way of analytical transformation of the obtained covariance matrix Ž19. of r and c into the covariance matrix of a and t . First, the auto-covariance matrices of the decay times t and amplitudes a are calculated. Then, the cross-covariance matrix of t and a can be found. Ži. As the decay times are only dependent on the vector c, one can obtain the truncated Taylor expansion of t with respect to c: Dt s Tc D c ,
Ž 20 .
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
316
where D c s D c k 4ks1, . . . , n , Dt s Dt k 4ks1, . . . , n and matrix Tc is derived from Eq. Ž6. as:
½
Tc s Et Ž c . rEcs Ž ytm .
K a, a4 s A r K r ,r 4 AX r q A r K r ,t 4 AXt
kq1
q A t K t ,r 4 AX r q A t K t ,t 4 AXt ,
ny 1
5
Ý Ž i q 1. c i Žytm . i
r
is1
.
Ž 21 .
ks1 , . . . , n ; ms1 , . . . , n
Ž 22.
Žii. The amplitudes a are dependent on both decay times t and coefficients r via Eq. Ž8.. Thus, truncated Taylor expansion of a with respect to t and r contains two components: D a s A r D r q A t Dt ,
Ž 23.
where D a s D a k 4ks1, . . . , n , D r s D rk 4ks1, . . . , n , and the matrices of derivatives A r and A t are obtained from Eq. Ž8.: p
½
A r s Ea Ž r ,t . rErs Ž y1 . Ž ytq . n
Ł Ž 1 y tqrt j .
js1, j/q
A t s EaE w Ž r ,t . rEts
Ž 28 .
Žiii. Derivation of the cross-covariance matrix of the decay times and amplitudes is straightforward: K a,t 4 s A r K r ,t 4 q A t K t ,t 4 .
K t ,t 4 s Tc K c ,c 4 TcX .
½
5
Ž 27.
where K r ,t 4 s Tc K r ,c 4 .
Combining Eqs. Ž20. and Ž21. and applying the standard error propagation principle w11x, the covariance matrix of the decay times takes the form:
r
Combining Eqs. Ž23. – Ž26. and applying the standard error propagation principle, the covariance matrix of the amplitude becomes:
Ž 29.
Eqs. Ž22., Ž27. and Ž29. immediately lead to the covariance matrix of the parameters a and t from the covariance matrix of the coefficients c and r, obtained by the maximum likelihood method. The variances of the amplitudes sA2 a k 4 and decay times sA2 t k 4 are represented by the diagonal elements of the matrices K a, a4 and Kt ,t 4 respectively:
sA2 a k 4 s K a k ,a k 4 , sA2 t k 4 s K t k ,t k 4 ,k s 1, . . . ,n. Ž 30. Software implementation of these transformations is relatively simple and very efficient. Examples of the estimation of the variance for two- and three-exponential decays are presented in Section 6.
p
,
Ž 24.
qs1 , . . . , n ; ps1 , . . . , n
4. Marquardt NLLS method
ny 1
Ý Žytq . kq1 rkq1
In this section, we shortly describe the Marquardt NLLS method w11x. If we denote the following vectors:
ks 0
F ) s f ) Ž t i . 4 is1 , . . . , N ,
r t p Ž t p y tq .
F s f Ž t i . 4 is1 , . . . , N ,
Ž 31 .
)Ž
n
=
5
Ł Ž 1 y tqrt j .
js1, j/q
, qs1 , . . . , n ; ps1 , . . . , n ; q/p
Ž 25.
where f t i ., i s 1, . . . , N is the number of photons, collected in the ith time channel of the multichannel analyzer Žexperimental fluorescence decay., the estimators of the amplitudes a, and decay times t in Eq. Ž1. are obtained by minimization of the following quadratic form: X
x 2 Ž a,t . s w F ) y F x W w F ) y F x s min, A t s Ea Ž r ,t . rEts
½
y
Ý rkq1Žytq .
k
where W is the diagonal matrix of weights: W s w Ž t i .4is1, . . . , N . The weight factor w Ž t i . depends on the experimental statistics, and for the Poisson statistics of the number of photons in a channel it becomes: w Ž t i . s 1rf ) Ž t i .. Function f Ž t i . s f Ž t i , a,t . is nonlinear with respect to the parameters a, t , and can be linearized by expanding in Taylor series. Keeping in the series only components with first order derivatives, one obtains:
ks 0 n
= y1rtq q
1r Ž tm y tq .
Ý ms 0 , m/q
n
r
Ł Ž 1 y tqrt j .
js1, j/q
Ž 32 .
ny 1
5
. qs1 , . . . , n ; ps1 , . . . , n ; qsp
Ž 26.
F ) s F0 q D aFa q Dt Ft ,
Ž 33.
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
where F0 s f Ž t i , a 0 ,t 0 . 4 is1 , . . . , N , Fa s EF Ž a 0 ,t 0 . rE as =g Ž z . d z
5
½H
ti
0
Ž 34.
exp y Ž t i y z . rt 0 k 4 ,
Ž 35.
ks1 , . . . , n ; is1 , . . . , N
Ft s EF Ž a 0 ,t 0 . rEts
½H
ti
0
a0 k Ž t i y z . rt 02k 4
=exp y Ž t i y z . rt 0 k 4 g Ž z . d z
5
, ks1 , . . . , n ; is1 , . . . , N
Ž 36. D a s D a k 4ks1, . . . , n and Dt s Dt k 4ks1, . . . , n are the vectors of the parameter increments, a 0 s a0 k 4ks1, . . . , n and t 0 s t 0 k 4ks1, . . . , n are the initial values of the parameters. Applying the linear least squares method to Eq. Ž33. yields the following set of linear equations with respect to the parameters increments D a, Dt : FaW Ž F ) y F0 .
X
FtW Ž F ) y F0 . s
X
Ž 1 q l . FaWFaX FtWFaX
FaWFtX
Ž1 q l.
FtWFtX
D aX , Dt X
Ž 37.
where the prime denotes the transpose operation. New vectors a and t are obtained: a s a 0 q D a and t s t 0 q Dt , and a new value of x 2 Ž a,t . is calculated. If x 2 Ž a,t . - x 2 Ž a 0 ,t 0 ., then new guesses a and t are accepted Ž a 0 s a, t 0 s t . and l decreases; if x 2 Ž a,t . G x 2 Ž a 0 ,t 0 ., new guesses are rejected and l increases. This iterative procedure is repeated until some final criteria will be satisfied. The covariance matrix K of the estimated parameters a and t is represented by the inverse matrix of the set Ž37. at the last iteration Žwhere l is supposed to be close to zero., provided that the global minimum with the true values of parameters has been achieved w11x. The variances of the amplitudes sA2 a k 4 and decay times sA2 t kyn 4 are represented by the diagonal elements of the matrix K:
sA2 a k 4 s K k k ,
k s 1, . . . ,n,
sA2 t kyn 4 s K k k ,
k s n q 1, . . . ,2 n.
Ž 38.
5. Bootstrap estimation of covariance matrix In recent years simulation techniques were proved to be very powerful for solving statistical problems, such as estimation of a bias, variance, confidence interval, etc.
317
w25x. This approach has been recently developed for fluorescence decay analysis w26x. In general, bootstrap approach in the application to the fluorescence decay analysis consists in computer generation of intensity histograms, provided that the number of counts in a histogram channel is a Gaussian or Poisson random value with mean and variance equal to the experimentally obtained number of counts in this channel. Every generated histogram Žwith the unique realization of statistical noise. is fitted and the results are stored over the predefined number of runs. Afterwards, the required statistical characteristics of the parameters can be calculated Žmean value, covariance matrix, etc... The statistical accuracy of the estimation is solely determined by the number of runs. In this paper let us call this algorithm ‘global bootstrap’. The main advantage of this approach is that it can be used for any method of analysis and any model Žnot only multi-exponential. of fluorescence decay. However, if the fitting method is relatively slow Žlike NLLS., it may take much time for obtaining reasonable statistical accuracy of the estimators. On the contrary, for non-iterative methods the fitting time is very low, and multiple repetition of simulation runs is not extremely time consuming. One more modification of the bootstrap algorithm for the phase plane method can be suggested. It is based on the fact that maximum likelihood estimators are asymptotically Gaussian w23x. First, applying the phase plane method ŽSection 2. to the experimentally obtained data, one gets the vector of estimators u Žthe vectors a and t can be readily found from Eqs. Ž6. and Ž8. and the covariance matrix Ry1 ŽEqs. Ž16. and Ž19... Taking into account Gaussian statistics of the elements of u, information, obtained at the first step, is adequate for further generating of the random Gaussian realizations u ) of the vector u directly, assuming that the vector of mean values is u and the covariance matrix is Ry1. That generating procedure can be performed in the following way w27x: 1. Find matrix B, such that BBX s Ry1 Žtriangular decomposition of Ry1 ., where the prime denotes the transpose operation. 2. Generate vector z s z k 4ks1, . . . ,2 n , of independent Gaussian variables with zero mean and variance, set to 1. 3. Calculate new Gaussian realization u ) of the vector u as: u ) s zBX q u. The vector u ) can then be used to calculate the decay times and amplitudes according to Eqs. Ž6. and Ž8.. As earlier, the results are stored over the predefined number of runs. Actually, this algorithm Žnamed ‘local bootstrap’. is another way of transformating the covariance matrix Ry1 of r and c into the covariance matrix of a and t . It avoids truncation of the Taylor series and calculations of the derivatives ŽSection 3. and can be recommended for more complicated, than sum of exponentials, decay laws. On the
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
318
other hand, this procedure is much faster than the global bootstrap, where every run starts from the simulation of the whole fluorescence decay.
6. Simulation and testing Simulation experiments have been performed in order to show the applicability of the proposed algorithms for the determination of the fluorescence decay parameters and their standard deviations. The configuration of simulation experiments is arranged, as described in Refs. w5,28x. The number of time channels is 256 or 512 on the detection range equal to 256 relative units, thus, ensuring the width of a channel as 1 or 0.5 relative units respectively. For the instrumental function g Ž t . the following expression is used: 2
g Ž t . s exp Ž y Ž t y ² t : . rS 2 ,
Ž 39.
where ² t : s 30 and S s 4. The counts of the function f Ž t . are calculated as a digital convolution of the discrete representation of g Ž t . with the true decay I Ž t .. The number of counts in the ith channel is the Poisson random variable with mean value equal to the exact value of the intensity in the ith channel. The number of counts in the
peak channel is set to either 10 4 or 5 = 10 3 for both instrumental function and fluorescence response. In every experiment the estimation has been held for 100 runs of the generated fluorescence and instrumental functions, each one with different realization of statistical noise. The results of every estimation were stored for the calculation of the mean value of the parameters: 100
MG t k 4 s
100
Ý it kr100,
MG a k 4 s
is1
Ý ia kr100,
and their variances: 100
sG2 t k 4 s
Ý it k2r100 y MG t k 4 2 , is1 100
sG2 a k 4 s
Ý ia2k r100 y MG ak 4 2 ,
Ž 41.
is1
where it k and i a k are the estimators, obtained at the ith run, k s 1, . . . ,n. This procedure is similar to the global bootstrap ŽSection 5., only instead of the experimentally obtained number of counts in a channel, the true Žtheoretically calculated. value is used for simulation. Moreover, at each run we calculate the standard deviation, obtained by the analytical sA ŽEq. Ž30. is used for the phase plane method, and Eq. Ž38. is used for the Marquardt NLLS
Table 1 Double-exponential fitting by the phase plane ŽPP. and Marquardt NLLS methods Method
M, s
Ž 40.
is1
Set I
Set II
a1 s 1
t 1 s 10
a2 s 1
t 2 s 20
a1 s 1
t 1 s 50
a2 s 1
t 2 s 150
Convoluted data PP MG sG M sA 4 s sA 4 M s B 4 s s B 4 NLLS MG sG M sA 4 s sA 4
1.01 0.066 0.032 0.0063 0.032 0.0064 0.93 0.044 0.040 0.0035
9.84 0.55 0.31 0.020 0.31 0.020 9.61 0.37 0.35 0.014
1.01 0.075 0.038 0.0065 0.038 0.066 1.08 0.048 0.044 0.0035
20.0 0.32 0.164 0.037 0.165 0.038 19.60 0.22 0.18 0.015
0.98 0.079 0.074 0.011 0.075 0.010 0.99 0.067 0.067 0.0069
49.1 2.88 2.11 0.21 2.13 0.21 49.50 2.42 2.54 0.14
1.03 0.083 0.072 0.011 0.070 0.010 1.01 0.071 0.071 0.0069
149 7.29 7.01 2.07 7.37 2.39 149.00 6.30 6.11 1.37
Data without convolution PP MG sG M sA 4 s sA 4 M s B 4 s s B 4 NLLS MG sG M sA 4 s sA 4
1.01 0.074 0.042 0.0080 0.043 0.0080 0.92 0.057 0.049 0.0035
10.10 0.59 0.40 0.027 0.41 0.028 9.60 0.39 0.39 0.016
0.99 0.084 0.050 0.0083 0.051 0.084 1.08 0.060 0.052 0.0036
20.00 0.38 0.21 0.049 0.22 0.051 19.50 0.28 0.21 0.019
1.01 0.072 0.066 0.0086 0.066 0.0084 1.00 0.055 0.051 0.0044
50.20 2.91 2.64 0.18 2.67 0.18 50.00 2.08 2.06 0.088
0.99 0.078 0.072 0.0087 0.072 0.0085 0.99 0.058 0.055 0.0044
151.00 6.62 6.15 1.56 6.50 1.72 150.00 4.89 4.57 0.81
Number of time channels is 256, peak counts 10 4 . Designations: M is the mean value of the estimator, s is the standard square deviation of the estimator.
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
319
Table 2 Triple-exponential fitting by the phase plane ŽPP. and Marquardt NLLS methods for shorter decay times N, P 256, 10
4
Method
M, s
a1 s 4
t1 s5
a2 s 2
t 2 s 15
a3 s 1
t 3 s 40
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4 MG sG M sA 4 s sA 4
4.21 0.27 0.088 0.019 0.092 0.019 3.91 0.23 0.20 0.034
5.07 0.40 0.19 0.023 0.19 0.023 4.91 0.33 0.31 0.037
1.92 0.23 0.10 0.025 0.11 0.026 2.06 0.22 0.20 0.040
15.70 1.81 0.53 0.19 0.56 0.20 14.50 1.13 1.01 0.18
0.98 0.090 0.021 0.012 0.021 0.013 1.05 0.056 0.051 0.0086
40.30 0.99 0.22 0.17 0.23 0.19 39.30 0.59 0.51 0.090
NLLS
512, 10 4
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.06 0.21 0.066 0.0084 0.069 0.0086
5.16 0.34 0.14 0.010 0.14 0.010
1.93 0.21 0.076 0.011 0.078 0.011
15.50 1.22 0.37 0.10 0.38 0.11
0.98 0.055 0.014 0.0048 0.014 0.0050
40.20 0.61 0.14 0.054 0.14 0.056
512, 5 = 10 3
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.12 0.26 0.10 0.021 0.11 0.021
5.52 0.43 0.21 0.023 0.21 0.023
1.77 0.24 0.11 0.028 0.11 0.029
16.50 1.95 0.69 0.24 0.72 0.26
0.95 0.093 0.026 0.015 0.027 0.016
40.60 1.07 0.28 0.20 0.30 0.23
Designations: N is the number of channels; P is the number of counts in the peak channel; M is the mean value of the estimator; s is the standard square deviation of the estimator.
Table 3 Triple-exponential fitting by the phase plane and Marquardt NLLS methods for longer decay times N, P 256, 10
4
Method
M, s
a1 s 4
t 1 s 10
a2 s 2
t 2 s 40
a3 s 1
t 3 s 120
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4 MG sG M sA 4 s sA 4
4.00 0.20 0.11 0.019 0.11 0.019 4.01 0.15 0.14 0.022
10.10 0.66 0.38 0.035 0.38 0.035 10.10 0.50 0.41 0.040
2.04 0.13 0.077 0.026 0.077 0.022 2.02 0.096 0.11 0.030
41.50 6.66 3.77 1.15 3.89 1.11 41.00 5.69 4.92 0.95
0.94 0.26 0.15 0.051 0.15 0.047 0.97 0.22 0.20 0.045
133 57.30 26.3 10.8 28.2 9.4 129 53.7 30.5 8.7
NLLS
512, 10 4
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.04 0.14 0.077 0.0098 0.080 0.0099
10.10 0.43 0.27 0.017 0.29 0.017
2.05 0.094 0.049 0.014 0.049 0.013
40.80 4.79 2.58 0.55 2.66 0.54
0.97 0.20 0.10 0.026 0.105 0.025
125 16.0 7.82 2.39 8.20 2.83
512, 5 = 10 3
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.06 0.17 0.11 0.018 0.11 0.017
10.40 0.60 0.37 0.032 0.38 0.032
2.05 0.13 0.080 0.026 0.079 0.022
43.60 6.15 4.11 0.97 4.22 0.93
0.87 0.24 0.16 0.043 0.16 0.038
137.00 41.8 21.21 6.54 23.90 6.02
Designations: N is the number of channels; P is the number of counts in the peak channel; M is the mean value of the estimator; s is the standard square deviation of the estimator.
320
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
method. and local bootstrap s B ŽSection 5. techniques. The results are also stored and the mean value Ž M sA 4 and M s B4. and standard deviation of the standard deviation estimator Ž s sA 4 and s s B4. are calculated. The results of the simulation experiments are given in Tables 1–3. We test fit of double and triple exponential decays with shorter decay times: t 1 s 10, t 2 s 20 Ždouble exponential decay: Table 1, set I., and t 1 s 5, t 2 s 15 t 3 s 40 Žtriple exponential decay: Table 2.; and longer decay times: t 1 s 50, t 2 s 150 Ždouble exponential decay: Table 1, set II., and t 1 s 10, t 2 s 40 t 3 s 120 Žtriple exponential decay: Table 3.. Fitting has been performed also for the undistorted by the convolution with the instrumental function g Ž t . Ž g Ž t . is the delta function. double exponential decay ŽTable 1.. For the triple exponential decay the influence of the number of time channels and noise level on the accuracy of estimation is explored ŽTables 2 and 3.. Simulation results show that the phase plane method mostly gives a reasonable accuracy of the estimation. Parameters, obtained by the phase plane algorithm, have been used as initial guesses for the iterative Marquardt NLLS procedure, ensuring better x 2 value in all cases. Better x 2 may not always correspond to the less error of the parameters, but, as a rule, Marquardt NLLS provides lower variance for the parameters, compared to the phase plane method. From the tables it is clearly seen that analytical standard deviation M sA 4 and standard deviation s , obtained by the global bootstrap, for the Marquardt NLLS procedure, perfectly correspond almost in all cases. Only triple exponential deconvolution slightly decreases the analytical estimators sA ŽTables 2 and 3., because the statistical noise, introduced by the instrumental function and not accounted by the method, becomes more prominent in that case. The same reason for the ambiguity holds for the phase plane method. As long as the covariance matrix V in Eq. Ž15. incorporates only the noise, produced by the pure fluorescence decay f ) Ž x . ŽEq. Ž18.. and not by the noise and covariance of the integrals: Fk) Ž x ., G k) Ž x ., k s 1, . . . ,n, the estimators of the variance, obtained by the analytical or local bootstrap approaches are smaller than the true values, which are, of course, dependent on all sources of noise. However, the order of magnitude for both estimators is the same and analytical estimators can be used for rough approximation of the variance for the phase plane method Žkeeping in mind that real values are higher that the analytically predicted ones.. Besides the initially incorrect matrix V, both methods introduce additional distortions into the estimation. Analytical estimators suffer from errors of truncation in the Taylor series. Local bootstrap estimators are statistically distorted, when the number of runs is insufficient. By increasing the number of runs, this drawback can always
be avoided. However, this may increase the processing time. Comparing the results of the analytical and local bootstrap estimations, good correspondence is obtained. As software implementation of the analytical approach is much faster, it may be recommended as a method of choice for multi-exponential decay analysis. References w1x D.V. O’Connor, D. Philips, Time-correlated Single Photon Counting, Academic Press, London, 1984. w2x J.N. Demas, Excited State Lifetime Measurements, Academic Press, New York, 1983. w3x D.F. Eaton, Pure Appl. Chem. 62 Ž1990. 1631. w4x A.R. Holzwarth, Meth. Enzymol. 246 Ž1995. 334. w5x J. Enderlein, R. Erdmann, Optics Comm. 134 Ž1997. 371. w6x Z. Zhang, K.T.V. Grattan, Y. Hu, A.W. Palmer, B.T. Meggitt, Rev. Sci. Instrum. 67 Ž1996. 2590. w7x J.R. Lakowicz, Principles of Fluorescence Spectroscopy, Plenum, New York, 1983. w8x D.V. O’Connor, W.R. Ware, J.C. Andre, J. Phys. Chem. 83 Ž1979. 1333. w9x V.V. Apanasovich, E.G. Novikov, J. Appl. Spectrosc. 56 Ž1992. 538. w10x J.A. Nelder, R. Mead, Comput. J. 7 Ž1965. 308. w11x P.R. Bevington, Data Reduction and Error Analysis for the Physical Sciences, McGraw-Hill, New York, 1969. w12x M. Ameloot, Meth. Enzymol. 210 Ž1992. 279. w13x E.W. Small, L.J. Libertini, D.W. Brown, J.R. Small, SPIE Proc. 1054 Ž1989. 36. w14x B. Valeur, Chem. Phys. 30 Ž1978. 85. w15x J.Y. Jezequel, M. Bouchy, J.C. Andre, Anal. Chem. 54 Ž1982. 2199. w16x V.V. Apanasovich, E.G. Novikov, Optics Comm. 78 Ž1990. 279. w17x J. Szamosi, Z.A. Schelly, J. Phys. Chem. 88 Ž1984. 3197. w18x J. Tang, J.R. Norris, Nucl. Instrum. Meth. Phys. Res. A 273 Ž1988. 338. w19x G.R. Phillips, E.M. Eyring, Anal. Chem. 60 Ž1988. 738. w20x D.T. Gillespie, Am. J. Phys. 51 Ž1983. 520. w21x M.G. Kendall, A. Stuart, The Advanced Theory of Statistics, vol. 2, McMillan, New York, 1979. w22x S.A. Aivazyan, I.S. Yenyukov, L.D. Meshalkin, Applied Statistics, Study of Relationships, Finansy Statistika, Moscow, 1985. w23x S. Zacks, The Theory Of Statistical Inference, Wiley, New York, 1971. w24x M.H. Pesaran, L.J. Slater, Dynamic Regression: Theory and Algortithms, Wiley, New York, 1984. w25x B.E. Efron, The Jackknife, The Bootstrap and Other Resampling Plans, SIAM, Philadelphia, PA, 1982. w26x M. Straume, S.G. Fraiser-Cadoret, M. Johnson, in: J.R. Lakowicz ŽEd.., Topics in Fluorescence Spectroscopy, Plenum, New York, 1991, p. 177. w27x J.H. Stapleton, Linear Statistical Models, Wiley, New York, 1995. w28x V.V. Apanasovich, E.G. Novikov, SPIE Proc. 2340 Ž1994. 59.
Optics Communications 151 Ž1998. 313–320
Full length article
Error estimation in the phase plane method of multi-exponential decay analysis )
Eugene G. Novikov
Systems Analysis Department, Byelorussian State UniÕersity, F. Skarina AÕenue 4, Minsk 220050, Belarus Received 20 October 1997; revised 20 January 1998; accepted 22 January 1998
Abstract The problems of error estimation for the parameters of a multi-exponential fluorescence decay law are considered. Two algorithms for the determination of the covariance matrix in the phase plane method, based on the analytical transformations and bootstrap approach, are suggested. Simulation experiments have been performed in order to show the applicability of the proposed algorithms for the determination of the fluorescence decay parameters and their standard deviations. q 1998 Published by Elsevier Science B.V. All rights reserved. Keywords: Fluorescence decay fitting; Multi-exponential decay law; Phase plane method; Prony method; Marquardt nonlinear least squares; Covariance matrix; Variance; Bootstrap
1. Introduction The development history of software algorithms for fluorescence decay analysis is fairly rich w1–4x, but many different aspects of this problem are still of particular interest w5,6x, in particular to develop more accurate, more reliable and faster fitting methods. In many applications fluorescence decay can be adequately approximated by a sum of exponentials w7x: n
IŽt. s
Ý ak exp Žytrt k . ,
probe, excited by a train of short light pulses, is detected as discrete numbers of photons, collected in a finite number of the time channels of a multichannel analyzer. The finite width of an excitation pulse and the delays in a detector make the true decay curve I Ž t . Žimpulse response of a fluorescent probe. convoluted with the instrumental function g Ž t . Žwhich itself is a convolution of an excitation pulse and detector response.. The detected fluorescence decay f Ž t . is then determined by a convolution integral:
Ž1.
ks 1
where n is the number of exponentials, a s a k 4ks1, . . . , n and t s t k 4ks1, . . . , n are the amplitudes and decay constants which have to be defined. The methods of time-resolved fluorescence spectroscopy w4,7x are widely used to obtain information about the time behaviour of fluorescence. One of the more commonly used techniques is time-correlated single photon counting w1x, where the fluorescence response of a
)
E-mail: [email protected]
f Žt. s
t
H0 I Ž t y z . g Ž z . d z ,
g Ž z . s 0,
z F 0.
Ž2.
A variety of methods for fluorescence data fitting has been exercised with different level of success w8,9x. Among them two groups are clearly prominent: iterative nonlinear least squares ŽNLLS. methods, directly optimizing some final criteria Žusually x 2 . Žsimplex method w2,10x, Marquardt method w11x, etc.., and non-iterative methods, based on the transformation of the detected decay to a linear algebraic set of equations with respect to unknown parameters ŽLaplace transforms method w12x, method of moments
0030-4018r98r$19.00 q 1998 Published by Elsevier Science B.V. All rights reserved. PII S 0 0 3 0 - 4 0 1 8 Ž 9 8 . 0 0 0 5 0 - 9
314
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
w13x, method of modulating functions w14x, phase plane method w2,15,16x, Prony’s method w5,17,18x, etc... The methods in the first group, granting high quality of fitting, are relatively slow and need initial guesses for the parameters to be known. When the initial guesses are far from the true values, the convergence may be very slow and time consuming. The non-iterative methods are not so accurate in fitting but extremely fast and normally do not need initial guesses. However, analysis of the existing techniques w8,9x reveals strict influence of the extra factors on the accuracy of fitting Žparameter of transformation for Laplace transform method, cut-off ratio for method of moments, kind of modulation function for method of modulating functions.. Our experience shows that the phase plane method is more stable in that respect, providing acceptable accuracy and high speed of processing. The common way of fluorescence decay analysis may be composed of non-iterative routine, generating initial guesses, and iterative procedure, starting the search with those initial values. The phase plane method is a good choice for the first step and Marquardt NLLS for the second one. One note should be made, concerning the principal sameness of the Prony’s method and phase plane method. Both of them yield the equation of a linear regression with respect to the unknown coefficients, which can be fit by the methods of linear least-squares minimization. The difference is in the way of transformation of the convolution integral Ž2. with the multi-exponential fluorescence response Ž1. into the equation of a linear regression. Prony’s method uses the technique of Z transformations w17,18x, and the phase plane method is based on the integration of the convolution Ž2. with the running upper limit w15,16x. A significant problem, that always accompanies the problem of point estimation, involves determination of the variance and confidence interval for the estimated parameters. This problem is prevalent for most methods mentioned in this paper. Comparably extensive statistical theory has been developed only for the NLLS methods w10,11,19x. A way of estimating the amplitudes and decay time variances has been proposed also for the method of moments w13x. But cut-off problems, which make this method very unstable, especially for multi-exponential decay fitting, become even more pronounced for the calculation of the variances. In the present paper two algorithms of error estimation for the phase plane method are proposed and compared to the estimators, provided by the Marquardt NLLS method. The further considerations are also appropriate for the Prony’s method.
equation Ž2. with the multi-exponential fluorescence response Ž1. to a linear equation with constant coefficients: n
n
Ý c k Fk Ž x . s Ý rk Gk Ž x . , ks 0
Ž3.
ks1
where F0 Ž x . s f Ž x ., x
ky1
Fk Ž x . s
H0 f Ž t .Ž x y t .
s
H0 g Ž t .Ž x y t .
x
r Ž k y 1 . ! d t ,Gk Ž x .
ky1
r Ž k y 1 . ! d t ,k s 1, . . . ,n.
Ž4. This transformation can also be applied to the analysis of multi-exponential decay without convolution Ž g Ž t . is the delta function.. In this case functions G k Ž x . take the form: G k Ž x . s x ky1r Ž k y 1 . !,
k s 1, . . . ,n.
Ž5.
Eq. Ž3. is a linear functional relationship with respect to the coefficients c k , k s 0, . . . ,n and rk , k s 1, . . . ,n, which must be fitted. The decay times t k can be obtained as roots of the following polynomial: n
Ý Žy1. k c kt kyn s 0,
Ž6.
ks 0
and the amplitudes a k are evaluated from the set of linear equations: n
rk s
k
Ý
Ý Žy1. kyj c jy1tmjyk ,
am
ms 1
Ž7.
js1
where t k , k s 1, . . . ,n are known from Eq. Ž6.. The analytical solution of the set Ž7. gives: ny 1
am s
Ý Žy1. k rkq1tmk ks0
5½
n
Ł Ž 1 y tmrt j .
js1, j/m
m s 1, . . . ,n.
5
,
Ž8.
As already mentioned, the curves f Ž x . and g Ž x . are detected in the finite number N of time channels: f ) Ž x i . and g ) Ž x i ., i s 1, . . . , N. Introducing the following matrices: F ) s Fk) Ž x i . 4 ks0 , . . . , n ; is1 , . . . , N ; G ) s G k) Ž x i . 4 ks1 , . . . , n ; is1 , . . . , N
Ž9.
and vectors c s c k 4 ks0 , . . . , n ;r s r k 4 ks1 , . . . , n ,
Ž 10.
Eq. Ž3. can be rewritten as: 2. Phase plane method In this section the phase plane method is shortly outlined. It is based on the transformation of the convolution
c F ) s rG ) .
Ž 11 .
Taking into account statistical distortions, naturally presented in the measured data, one has: e s c F ) I rG ) ,
Ž 12.
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
where e is the vector of residuals, which is, according to the central limit theorem w20x, Gaussian with mean 0 and covariance function: cov Ž e . s ² eeX : s V.
Ž 13.
The coefficients of linear relationship Ž3. c s c k 4ks0, . . . , n and r s rk 4ks1, . . . , n can be obtained by the maximum likelihood method w5,21x, maximizing the following probability density function: P Ž e . A exp Ž yeVy1 eXr2 . ,
Ž 14.
or, similarly, minimizing the quadratic form w22x:
eigenvalue and corresponding eigenvector u of the symn metric matrix R, when Ý2ks0 u 2k s 1 w21x. In our case it is certainly not correct, as the elements of F and G are highly correlated, and their variances are considerably different. On the contrary, since integrating Eq. Ž4. normally smoothes statistical distortions in the data, the relative variance of the elements Fk) and G k) , k s 1, . . . ,n of the vectors F and G is much smaller than the relative variance of the zero element F0) of the vector F. Therefore, we may put c 0 s 1 in Eq. Ž3., thus, getting the equation of a linear regression: n
X
eV I1 eX s Ž c F ) I rG ) . Vy1 Ž c F ) I rG ) .
f )Ž x. s
s u R u s min,
Ž 15.
X
where u s c, r 4 is the vector-string, consisting of two sub-strings, and Rs
F ) V I 1 F )X F ) V I 1 G )X
F ) V I 1 G )X G ) V I 1 G )X
Ž 16.
is the matrix, formed by four sub-matrices of lower dimensions; the prime denotes the transpose operation. To avoid trivial solution Ž c k s 0 and r k s 0., one may set, for n example, Ý2ks0 u 2k s Ý nks0 c k2 q Ý nks1 r k2 s 1, or u 0 s c 0 s 1. Differentiating Eq. Ž15. with respect to the elements of u results in a set of equations, that can be solved numerically. The solution gives the maximum likelihood estimators of the coefficients c and r in Eq. Ž3.. According to the invariance principle w23x, the parameters a and t , calculated from the vectors c and r via Eqs. Ž6. and Ž8. respectively, are estimators, providing the maximum to the likelihood function Ž14., as well. The problem of the maximum likelihood estimation is mainly related to the complex structure of the matrix V, and, particularly, to the dependence of the elements of V on the unknown coefficients c and r. Iterative algorithms w24x, or introducing a reasonable approximation for V w5x are used in order to solve this problem. Both approaches require to calculate the inverse matrix Vy1, the dimension of which corresponds to the number of channels N and may be very large Žhundreds or thousands of channels.. This gives rise to problems of numerical implementation, resulting in loss of accuracy. Moreover, the first approach may cause obstacles, naturally accompanying the iterative algorithms, like arriving into local minima and increased time of calculation. In the second one, proper approximation for the matrix V is not always obvious. Two special cases, when minimization of Eq. Ž15. is straightforward, should be mentioned at this point. If all elements of the vectors F and G are independent of each other and identically distributed with the same variance, then the problem of minimization of the quadratic form Ž15. is equivalent to the determination of the minimal
n
Ý c k Fk) Ž x . y Ý rkGk) Ž x . , ks 1
X
315
Ž 17.
ks1
where functions Fk) Ž x . and G k) Ž x ., k s 1, . . . ,n are assumed to be noiseless compared to F0) , and, consequently, statistically independent of each other. In this case, maximum likelihood estimators of the coefficients c k and rk , k s 1, . . . ,n are approximately equivalent to the estimators, found by the linear least squares method w21,22x, directly applied to the linear regression equation Ž17.. This approximation allows for easy analytical calculation of the inverse matrix of V, which is diagonal in that case: y1 y1 ) Vy1 s Viy1 j 4 ;Vi j s 0,i / j;Vi j s 1rf Ž x i . ,i s j.
Ž 18. Testing, undertaken in the Section 6, shows that this approximation ensures extremely fast and reasonably accurate fittings. Hence, further we shall assume that c 0 s 1, and vector c has the form: c s c k 4ks1, . . . , n .
3. Analytical estimation of the covariance matrix The covariance matrix of the maximum likelihood estimators takes the form w22,23x: Ry1 s
K c ,c 4 K r ,c 4
K c ,r 4 , K r ,r 4
Ž 19.
where K c,c4 and K r,r 4 are the auto-covariance matrices of the coefficients c and r respectively, and K c,r 4 and K r,c4 are the cross-covariance matrices of c and r. In this section we consider a way of analytical transformation of the obtained covariance matrix Ž19. of r and c into the covariance matrix of a and t . First, the auto-covariance matrices of the decay times t and amplitudes a are calculated. Then, the cross-covariance matrix of t and a can be found. Ži. As the decay times are only dependent on the vector c, one can obtain the truncated Taylor expansion of t with respect to c: Dt s Tc D c ,
Ž 20 .
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
316
where D c s D c k 4ks1, . . . , n , Dt s Dt k 4ks1, . . . , n and matrix Tc is derived from Eq. Ž6. as:
½
Tc s Et Ž c . rEcs Ž ytm .
K a, a4 s A r K r ,r 4 AX r q A r K r ,t 4 AXt
kq1
q A t K t ,r 4 AX r q A t K t ,t 4 AXt ,
ny 1
5
Ý Ž i q 1. c i Žytm . i
r
is1
.
Ž 21 .
ks1 , . . . , n ; ms1 , . . . , n
Ž 22.
Žii. The amplitudes a are dependent on both decay times t and coefficients r via Eq. Ž8.. Thus, truncated Taylor expansion of a with respect to t and r contains two components: D a s A r D r q A t Dt ,
Ž 23.
where D a s D a k 4ks1, . . . , n , D r s D rk 4ks1, . . . , n , and the matrices of derivatives A r and A t are obtained from Eq. Ž8.: p
½
A r s Ea Ž r ,t . rErs Ž y1 . Ž ytq . n
Ł Ž 1 y tqrt j .
js1, j/q
A t s EaE w Ž r ,t . rEts
Ž 28 .
Žiii. Derivation of the cross-covariance matrix of the decay times and amplitudes is straightforward: K a,t 4 s A r K r ,t 4 q A t K t ,t 4 .
K t ,t 4 s Tc K c ,c 4 TcX .
½
5
Ž 27.
where K r ,t 4 s Tc K r ,c 4 .
Combining Eqs. Ž20. and Ž21. and applying the standard error propagation principle w11x, the covariance matrix of the decay times takes the form:
r
Combining Eqs. Ž23. – Ž26. and applying the standard error propagation principle, the covariance matrix of the amplitude becomes:
Ž 29.
Eqs. Ž22., Ž27. and Ž29. immediately lead to the covariance matrix of the parameters a and t from the covariance matrix of the coefficients c and r, obtained by the maximum likelihood method. The variances of the amplitudes sA2 a k 4 and decay times sA2 t k 4 are represented by the diagonal elements of the matrices K a, a4 and Kt ,t 4 respectively:
sA2 a k 4 s K a k ,a k 4 , sA2 t k 4 s K t k ,t k 4 ,k s 1, . . . ,n. Ž 30. Software implementation of these transformations is relatively simple and very efficient. Examples of the estimation of the variance for two- and three-exponential decays are presented in Section 6.
p
,
Ž 24.
qs1 , . . . , n ; ps1 , . . . , n
4. Marquardt NLLS method
ny 1
Ý Žytq . kq1 rkq1
In this section, we shortly describe the Marquardt NLLS method w11x. If we denote the following vectors:
ks 0
F ) s f ) Ž t i . 4 is1 , . . . , N ,
r t p Ž t p y tq .
F s f Ž t i . 4 is1 , . . . , N ,
Ž 31 .
)Ž
n
=
5
Ł Ž 1 y tqrt j .
js1, j/q
, qs1 , . . . , n ; ps1 , . . . , n ; q/p
Ž 25.
where f t i ., i s 1, . . . , N is the number of photons, collected in the ith time channel of the multichannel analyzer Žexperimental fluorescence decay., the estimators of the amplitudes a, and decay times t in Eq. Ž1. are obtained by minimization of the following quadratic form: X
x 2 Ž a,t . s w F ) y F x W w F ) y F x s min, A t s Ea Ž r ,t . rEts
½
y
Ý rkq1Žytq .
k
where W is the diagonal matrix of weights: W s w Ž t i .4is1, . . . , N . The weight factor w Ž t i . depends on the experimental statistics, and for the Poisson statistics of the number of photons in a channel it becomes: w Ž t i . s 1rf ) Ž t i .. Function f Ž t i . s f Ž t i , a,t . is nonlinear with respect to the parameters a, t , and can be linearized by expanding in Taylor series. Keeping in the series only components with first order derivatives, one obtains:
ks 0 n
= y1rtq q
1r Ž tm y tq .
Ý ms 0 , m/q
n
r
Ł Ž 1 y tqrt j .
js1, j/q
Ž 32 .
ny 1
5
. qs1 , . . . , n ; ps1 , . . . , n ; qsp
Ž 26.
F ) s F0 q D aFa q Dt Ft ,
Ž 33.
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
where F0 s f Ž t i , a 0 ,t 0 . 4 is1 , . . . , N , Fa s EF Ž a 0 ,t 0 . rE as =g Ž z . d z
5
½H
ti
0
Ž 34.
exp y Ž t i y z . rt 0 k 4 ,
Ž 35.
ks1 , . . . , n ; is1 , . . . , N
Ft s EF Ž a 0 ,t 0 . rEts
½H
ti
0
a0 k Ž t i y z . rt 02k 4
=exp y Ž t i y z . rt 0 k 4 g Ž z . d z
5
, ks1 , . . . , n ; is1 , . . . , N
Ž 36. D a s D a k 4ks1, . . . , n and Dt s Dt k 4ks1, . . . , n are the vectors of the parameter increments, a 0 s a0 k 4ks1, . . . , n and t 0 s t 0 k 4ks1, . . . , n are the initial values of the parameters. Applying the linear least squares method to Eq. Ž33. yields the following set of linear equations with respect to the parameters increments D a, Dt : FaW Ž F ) y F0 .
X
FtW Ž F ) y F0 . s
X
Ž 1 q l . FaWFaX FtWFaX
FaWFtX
Ž1 q l.
FtWFtX
D aX , Dt X
Ž 37.
where the prime denotes the transpose operation. New vectors a and t are obtained: a s a 0 q D a and t s t 0 q Dt , and a new value of x 2 Ž a,t . is calculated. If x 2 Ž a,t . - x 2 Ž a 0 ,t 0 ., then new guesses a and t are accepted Ž a 0 s a, t 0 s t . and l decreases; if x 2 Ž a,t . G x 2 Ž a 0 ,t 0 ., new guesses are rejected and l increases. This iterative procedure is repeated until some final criteria will be satisfied. The covariance matrix K of the estimated parameters a and t is represented by the inverse matrix of the set Ž37. at the last iteration Žwhere l is supposed to be close to zero., provided that the global minimum with the true values of parameters has been achieved w11x. The variances of the amplitudes sA2 a k 4 and decay times sA2 t kyn 4 are represented by the diagonal elements of the matrix K:
sA2 a k 4 s K k k ,
k s 1, . . . ,n,
sA2 t kyn 4 s K k k ,
k s n q 1, . . . ,2 n.
Ž 38.
5. Bootstrap estimation of covariance matrix In recent years simulation techniques were proved to be very powerful for solving statistical problems, such as estimation of a bias, variance, confidence interval, etc.
317
w25x. This approach has been recently developed for fluorescence decay analysis w26x. In general, bootstrap approach in the application to the fluorescence decay analysis consists in computer generation of intensity histograms, provided that the number of counts in a histogram channel is a Gaussian or Poisson random value with mean and variance equal to the experimentally obtained number of counts in this channel. Every generated histogram Žwith the unique realization of statistical noise. is fitted and the results are stored over the predefined number of runs. Afterwards, the required statistical characteristics of the parameters can be calculated Žmean value, covariance matrix, etc... The statistical accuracy of the estimation is solely determined by the number of runs. In this paper let us call this algorithm ‘global bootstrap’. The main advantage of this approach is that it can be used for any method of analysis and any model Žnot only multi-exponential. of fluorescence decay. However, if the fitting method is relatively slow Žlike NLLS., it may take much time for obtaining reasonable statistical accuracy of the estimators. On the contrary, for non-iterative methods the fitting time is very low, and multiple repetition of simulation runs is not extremely time consuming. One more modification of the bootstrap algorithm for the phase plane method can be suggested. It is based on the fact that maximum likelihood estimators are asymptotically Gaussian w23x. First, applying the phase plane method ŽSection 2. to the experimentally obtained data, one gets the vector of estimators u Žthe vectors a and t can be readily found from Eqs. Ž6. and Ž8. and the covariance matrix Ry1 ŽEqs. Ž16. and Ž19... Taking into account Gaussian statistics of the elements of u, information, obtained at the first step, is adequate for further generating of the random Gaussian realizations u ) of the vector u directly, assuming that the vector of mean values is u and the covariance matrix is Ry1. That generating procedure can be performed in the following way w27x: 1. Find matrix B, such that BBX s Ry1 Žtriangular decomposition of Ry1 ., where the prime denotes the transpose operation. 2. Generate vector z s z k 4ks1, . . . ,2 n , of independent Gaussian variables with zero mean and variance, set to 1. 3. Calculate new Gaussian realization u ) of the vector u as: u ) s zBX q u. The vector u ) can then be used to calculate the decay times and amplitudes according to Eqs. Ž6. and Ž8.. As earlier, the results are stored over the predefined number of runs. Actually, this algorithm Žnamed ‘local bootstrap’. is another way of transformating the covariance matrix Ry1 of r and c into the covariance matrix of a and t . It avoids truncation of the Taylor series and calculations of the derivatives ŽSection 3. and can be recommended for more complicated, than sum of exponentials, decay laws. On the
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
318
other hand, this procedure is much faster than the global bootstrap, where every run starts from the simulation of the whole fluorescence decay.
6. Simulation and testing Simulation experiments have been performed in order to show the applicability of the proposed algorithms for the determination of the fluorescence decay parameters and their standard deviations. The configuration of simulation experiments is arranged, as described in Refs. w5,28x. The number of time channels is 256 or 512 on the detection range equal to 256 relative units, thus, ensuring the width of a channel as 1 or 0.5 relative units respectively. For the instrumental function g Ž t . the following expression is used: 2
g Ž t . s exp Ž y Ž t y ² t : . rS 2 ,
Ž 39.
where ² t : s 30 and S s 4. The counts of the function f Ž t . are calculated as a digital convolution of the discrete representation of g Ž t . with the true decay I Ž t .. The number of counts in the ith channel is the Poisson random variable with mean value equal to the exact value of the intensity in the ith channel. The number of counts in the
peak channel is set to either 10 4 or 5 = 10 3 for both instrumental function and fluorescence response. In every experiment the estimation has been held for 100 runs of the generated fluorescence and instrumental functions, each one with different realization of statistical noise. The results of every estimation were stored for the calculation of the mean value of the parameters: 100
MG t k 4 s
100
Ý it kr100,
MG a k 4 s
is1
Ý ia kr100,
and their variances: 100
sG2 t k 4 s
Ý it k2r100 y MG t k 4 2 , is1 100
sG2 a k 4 s
Ý ia2k r100 y MG ak 4 2 ,
Ž 41.
is1
where it k and i a k are the estimators, obtained at the ith run, k s 1, . . . ,n. This procedure is similar to the global bootstrap ŽSection 5., only instead of the experimentally obtained number of counts in a channel, the true Žtheoretically calculated. value is used for simulation. Moreover, at each run we calculate the standard deviation, obtained by the analytical sA ŽEq. Ž30. is used for the phase plane method, and Eq. Ž38. is used for the Marquardt NLLS
Table 1 Double-exponential fitting by the phase plane ŽPP. and Marquardt NLLS methods Method
M, s
Ž 40.
is1
Set I
Set II
a1 s 1
t 1 s 10
a2 s 1
t 2 s 20
a1 s 1
t 1 s 50
a2 s 1
t 2 s 150
Convoluted data PP MG sG M sA 4 s sA 4 M s B 4 s s B 4 NLLS MG sG M sA 4 s sA 4
1.01 0.066 0.032 0.0063 0.032 0.0064 0.93 0.044 0.040 0.0035
9.84 0.55 0.31 0.020 0.31 0.020 9.61 0.37 0.35 0.014
1.01 0.075 0.038 0.0065 0.038 0.066 1.08 0.048 0.044 0.0035
20.0 0.32 0.164 0.037 0.165 0.038 19.60 0.22 0.18 0.015
0.98 0.079 0.074 0.011 0.075 0.010 0.99 0.067 0.067 0.0069
49.1 2.88 2.11 0.21 2.13 0.21 49.50 2.42 2.54 0.14
1.03 0.083 0.072 0.011 0.070 0.010 1.01 0.071 0.071 0.0069
149 7.29 7.01 2.07 7.37 2.39 149.00 6.30 6.11 1.37
Data without convolution PP MG sG M sA 4 s sA 4 M s B 4 s s B 4 NLLS MG sG M sA 4 s sA 4
1.01 0.074 0.042 0.0080 0.043 0.0080 0.92 0.057 0.049 0.0035
10.10 0.59 0.40 0.027 0.41 0.028 9.60 0.39 0.39 0.016
0.99 0.084 0.050 0.0083 0.051 0.084 1.08 0.060 0.052 0.0036
20.00 0.38 0.21 0.049 0.22 0.051 19.50 0.28 0.21 0.019
1.01 0.072 0.066 0.0086 0.066 0.0084 1.00 0.055 0.051 0.0044
50.20 2.91 2.64 0.18 2.67 0.18 50.00 2.08 2.06 0.088
0.99 0.078 0.072 0.0087 0.072 0.0085 0.99 0.058 0.055 0.0044
151.00 6.62 6.15 1.56 6.50 1.72 150.00 4.89 4.57 0.81
Number of time channels is 256, peak counts 10 4 . Designations: M is the mean value of the estimator, s is the standard square deviation of the estimator.
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
319
Table 2 Triple-exponential fitting by the phase plane ŽPP. and Marquardt NLLS methods for shorter decay times N, P 256, 10
4
Method
M, s
a1 s 4
t1 s5
a2 s 2
t 2 s 15
a3 s 1
t 3 s 40
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4 MG sG M sA 4 s sA 4
4.21 0.27 0.088 0.019 0.092 0.019 3.91 0.23 0.20 0.034
5.07 0.40 0.19 0.023 0.19 0.023 4.91 0.33 0.31 0.037
1.92 0.23 0.10 0.025 0.11 0.026 2.06 0.22 0.20 0.040
15.70 1.81 0.53 0.19 0.56 0.20 14.50 1.13 1.01 0.18
0.98 0.090 0.021 0.012 0.021 0.013 1.05 0.056 0.051 0.0086
40.30 0.99 0.22 0.17 0.23 0.19 39.30 0.59 0.51 0.090
NLLS
512, 10 4
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.06 0.21 0.066 0.0084 0.069 0.0086
5.16 0.34 0.14 0.010 0.14 0.010
1.93 0.21 0.076 0.011 0.078 0.011
15.50 1.22 0.37 0.10 0.38 0.11
0.98 0.055 0.014 0.0048 0.014 0.0050
40.20 0.61 0.14 0.054 0.14 0.056
512, 5 = 10 3
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.12 0.26 0.10 0.021 0.11 0.021
5.52 0.43 0.21 0.023 0.21 0.023
1.77 0.24 0.11 0.028 0.11 0.029
16.50 1.95 0.69 0.24 0.72 0.26
0.95 0.093 0.026 0.015 0.027 0.016
40.60 1.07 0.28 0.20 0.30 0.23
Designations: N is the number of channels; P is the number of counts in the peak channel; M is the mean value of the estimator; s is the standard square deviation of the estimator.
Table 3 Triple-exponential fitting by the phase plane and Marquardt NLLS methods for longer decay times N, P 256, 10
4
Method
M, s
a1 s 4
t 1 s 10
a2 s 2
t 2 s 40
a3 s 1
t 3 s 120
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4 MG sG M sA 4 s sA 4
4.00 0.20 0.11 0.019 0.11 0.019 4.01 0.15 0.14 0.022
10.10 0.66 0.38 0.035 0.38 0.035 10.10 0.50 0.41 0.040
2.04 0.13 0.077 0.026 0.077 0.022 2.02 0.096 0.11 0.030
41.50 6.66 3.77 1.15 3.89 1.11 41.00 5.69 4.92 0.95
0.94 0.26 0.15 0.051 0.15 0.047 0.97 0.22 0.20 0.045
133 57.30 26.3 10.8 28.2 9.4 129 53.7 30.5 8.7
NLLS
512, 10 4
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.04 0.14 0.077 0.0098 0.080 0.0099
10.10 0.43 0.27 0.017 0.29 0.017
2.05 0.094 0.049 0.014 0.049 0.013
40.80 4.79 2.58 0.55 2.66 0.54
0.97 0.20 0.10 0.026 0.105 0.025
125 16.0 7.82 2.39 8.20 2.83
512, 5 = 10 3
PP
MG sG M sA 4 s sA 4 M s B 4 s s B 4
4.06 0.17 0.11 0.018 0.11 0.017
10.40 0.60 0.37 0.032 0.38 0.032
2.05 0.13 0.080 0.026 0.079 0.022
43.60 6.15 4.11 0.97 4.22 0.93
0.87 0.24 0.16 0.043 0.16 0.038
137.00 41.8 21.21 6.54 23.90 6.02
Designations: N is the number of channels; P is the number of counts in the peak channel; M is the mean value of the estimator; s is the standard square deviation of the estimator.
320
E.G. NoÕikoÕr Optics Communications 151 (1998) 313–320
method. and local bootstrap s B ŽSection 5. techniques. The results are also stored and the mean value Ž M sA 4 and M s B4. and standard deviation of the standard deviation estimator Ž s sA 4 and s s B4. are calculated. The results of the simulation experiments are given in Tables 1–3. We test fit of double and triple exponential decays with shorter decay times: t 1 s 10, t 2 s 20 Ždouble exponential decay: Table 1, set I., and t 1 s 5, t 2 s 15 t 3 s 40 Žtriple exponential decay: Table 2.; and longer decay times: t 1 s 50, t 2 s 150 Ždouble exponential decay: Table 1, set II., and t 1 s 10, t 2 s 40 t 3 s 120 Žtriple exponential decay: Table 3.. Fitting has been performed also for the undistorted by the convolution with the instrumental function g Ž t . Ž g Ž t . is the delta function. double exponential decay ŽTable 1.. For the triple exponential decay the influence of the number of time channels and noise level on the accuracy of estimation is explored ŽTables 2 and 3.. Simulation results show that the phase plane method mostly gives a reasonable accuracy of the estimation. Parameters, obtained by the phase plane algorithm, have been used as initial guesses for the iterative Marquardt NLLS procedure, ensuring better x 2 value in all cases. Better x 2 may not always correspond to the less error of the parameters, but, as a rule, Marquardt NLLS provides lower variance for the parameters, compared to the phase plane method. From the tables it is clearly seen that analytical standard deviation M sA 4 and standard deviation s , obtained by the global bootstrap, for the Marquardt NLLS procedure, perfectly correspond almost in all cases. Only triple exponential deconvolution slightly decreases the analytical estimators sA ŽTables 2 and 3., because the statistical noise, introduced by the instrumental function and not accounted by the method, becomes more prominent in that case. The same reason for the ambiguity holds for the phase plane method. As long as the covariance matrix V in Eq. Ž15. incorporates only the noise, produced by the pure fluorescence decay f ) Ž x . ŽEq. Ž18.. and not by the noise and covariance of the integrals: Fk) Ž x ., G k) Ž x ., k s 1, . . . ,n, the estimators of the variance, obtained by the analytical or local bootstrap approaches are smaller than the true values, which are, of course, dependent on all sources of noise. However, the order of magnitude for both estimators is the same and analytical estimators can be used for rough approximation of the variance for the phase plane method Žkeeping in mind that real values are higher that the analytically predicted ones.. Besides the initially incorrect matrix V, both methods introduce additional distortions into the estimation. Analytical estimators suffer from errors of truncation in the Taylor series. Local bootstrap estimators are statistically distorted, when the number of runs is insufficient. By increasing the number of runs, this drawback can always
be avoided. However, this may increase the processing time. Comparing the results of the analytical and local bootstrap estimations, good correspondence is obtained. As software implementation of the analytical approach is much faster, it may be recommended as a method of choice for multi-exponential decay analysis. References w1x D.V. O’Connor, D. Philips, Time-correlated Single Photon Counting, Academic Press, London, 1984. w2x J.N. Demas, Excited State Lifetime Measurements, Academic Press, New York, 1983. w3x D.F. Eaton, Pure Appl. Chem. 62 Ž1990. 1631. w4x A.R. Holzwarth, Meth. Enzymol. 246 Ž1995. 334. w5x J. Enderlein, R. Erdmann, Optics Comm. 134 Ž1997. 371. w6x Z. Zhang, K.T.V. Grattan, Y. Hu, A.W. Palmer, B.T. Meggitt, Rev. Sci. Instrum. 67 Ž1996. 2590. w7x J.R. Lakowicz, Principles of Fluorescence Spectroscopy, Plenum, New York, 1983. w8x D.V. O’Connor, W.R. Ware, J.C. Andre, J. Phys. Chem. 83 Ž1979. 1333. w9x V.V. Apanasovich, E.G. Novikov, J. Appl. Spectrosc. 56 Ž1992. 538. w10x J.A. Nelder, R. Mead, Comput. J. 7 Ž1965. 308. w11x P.R. Bevington, Data Reduction and Error Analysis for the Physical Sciences, McGraw-Hill, New York, 1969. w12x M. Ameloot, Meth. Enzymol. 210 Ž1992. 279. w13x E.W. Small, L.J. Libertini, D.W. Brown, J.R. Small, SPIE Proc. 1054 Ž1989. 36. w14x B. Valeur, Chem. Phys. 30 Ž1978. 85. w15x J.Y. Jezequel, M. Bouchy, J.C. Andre, Anal. Chem. 54 Ž1982. 2199. w16x V.V. Apanasovich, E.G. Novikov, Optics Comm. 78 Ž1990. 279. w17x J. Szamosi, Z.A. Schelly, J. Phys. Chem. 88 Ž1984. 3197. w18x J. Tang, J.R. Norris, Nucl. Instrum. Meth. Phys. Res. A 273 Ž1988. 338. w19x G.R. Phillips, E.M. Eyring, Anal. Chem. 60 Ž1988. 738. w20x D.T. Gillespie, Am. J. Phys. 51 Ž1983. 520. w21x M.G. Kendall, A. Stuart, The Advanced Theory of Statistics, vol. 2, McMillan, New York, 1979. w22x S.A. Aivazyan, I.S. Yenyukov, L.D. Meshalkin, Applied Statistics, Study of Relationships, Finansy Statistika, Moscow, 1985. w23x S. Zacks, The Theory Of Statistical Inference, Wiley, New York, 1971. w24x M.H. Pesaran, L.J. Slater, Dynamic Regression: Theory and Algortithms, Wiley, New York, 1984. w25x B.E. Efron, The Jackknife, The Bootstrap and Other Resampling Plans, SIAM, Philadelphia, PA, 1982. w26x M. Straume, S.G. Fraiser-Cadoret, M. Johnson, in: J.R. Lakowicz ŽEd.., Topics in Fluorescence Spectroscopy, Plenum, New York, 1991, p. 177. w27x J.H. Stapleton, Linear Statistical Models, Wiley, New York, 1995. w28x V.V. Apanasovich, E.G. Novikov, SPIE Proc. 2340 Ž1994. 59.