INTRODUCTION

TOP ABSTRACT INTRODUCTION THE MODEL DATA SIMULATIONS AND ANALYSIS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX REFERENCES

Molecular motions in proteins occur over a broad range in time (picoseconds to nanoseconds) and considerable effort has been and is being expended to detect and quantify them in the hope of understanding the role of intra-protein dynamics in protein function. The advent of molecular dynamics simulations (Van Gunsteren et al., 1993) has introduced a new imperative for accurate determination of the rates and amplitudes of backbone and side chain mobility. Current computational capability restricts simulations to at most a few nanoseconds in the majority of cases. Fortuitously, the magnetic resonance (NMR, EPR) and fluorescence techniques most commonly used to measure protein dynamics are optimally suited for quantitation of picosecond to nanosecond events. Fluorescence anisotropy decay measurements have been widely applied, in part because of the relative ease and low cost of the measurements and the increased sophistication and precision possible with modern laser-based time-resolved fluorescence techniques. However, analysis of fluorescence anisotropy decay data to extract meaningful parameters of motion for the construction of physical models of protein dynamics is always unavoidably compromised by the convolution of the instrument response function with the fluorescence decay data per se. The problem is particularly difficult for the recovery of events occurring in the low picosecond regime (<100 ps). Accordingly, reliable methods of data analysis are essential if the recovered parameters are to yield credible physical insight. In this regard, the discretization scheme used in the analysis of time-correlated single photon counting data for deconvoluting the instrument response function (IRF) is particularly important. This problem was recognized in the seminal paper of Grinvald and Steinberg (1974), and subsequently by McKinnon et al. (1977) who introduced a discretization scheme based on a local linear approximation of the IRF. A similar approach was described by Wahl (1979) in his pioneering work on the numerical problems related specifically to anisotropy decay analysis.

Recently, progress has been made in designing discretization schemes for convolution integrals which provide more accurate estimation of decay parameters (Holtom, 1990; Vee et al., 1993; Bajzer et al., 1995). These new discretization schemes have been tested using simulations in the simpler case of lifetime estimation from fluorescence intensity decay data. It was found that discretizations based on quadratic (QA) and cubic (CA) local approximations result in more accurate determination of short lifetimes than the discretizations based on a linear approximation (LA) (Vee et al., 1993; Bajzer et al., 1995). This result is not entirely expected because, depending on the level of noise by which the function is corrupted, its integral is sometimes more accurately estimated by lower-order approximations. Another comparative study (Periasamy, 1988) has demonstrated that the widely used applied Grinvald-Steinberg (GS) discretization scheme (Grinvald and Steinberg, 1974) yields less accurate results than discretization based on a local linear approximation.

The purpose of the present study is to investigate what can be gained in terms of the accuracy of parameter estimation in cases of complex anisotropy decays when more elaborate discretization schemes (QA or CA) are used. In particular, the ability of various discretization schemes to recover very short rotational relaxation times and lifetimes will be addressed. Additionally, a procedure that may allow discrimination between homogeneous and heterogeneous anisotropy decays will be discussed. This latter problem is known to be very difficult (Ludescher et al., 1987; Bialik et al., 1998). However, we show that the judicious use of various information criteria in combination with more elaborate discretization schemes facilitates the discrimination of these two decay models. In discriminating heterogeneous anisotropy decay from the analytically simpler homogeneous decay, Bialik et al. (1998) have suggested an approach based on carefully designed physical constraints and statistical criteria of goodness-of-fit. We consider the idea of applying Akaike (AIC), Bayesian (BIC), and Hannan-Quinn information criteria (HQIC), (cf. Walter and Pronzato, 1997; Davidian and Gallant, 1993). Our results suggest that the AIC is the best criterion for this purpose and could complement the approach of Bialik et al. (1998) and the analytical approach of Ludescher et al. (1987).

To our knowledge, previous simulations for anisotropy decay analysis (cf. Wahl, 1979; Ludescher et al., 1987; Bialik et al., 1998) did not address the issue of improving the accuracy in the estimation of the rotational relaxation times by applying more accurate discretization schemes. We demonstrate that QA and CA discretizations (in an improved version of those described in Bajzer et al., 1995) are better at recovering short lifetimes and short rotational relaxation times than GS and LA discretizations. We have also addressed the issue of whether it is more favorable to first determine lifetimes and their amplitudes from total intensity decay data and subsequently by fixing these parameters, estimate rotational relaxation times with their preexponentials from parallely and perpendicularly polarized intensity components, or perform simultaneous estimation of all parameters from the polarized intensity components. The results of our simulation suggest that the former option is preferable in some instances.

	THE MODEL

TOP ABSTRACT INTRODUCTION THE MODEL DATA SIMULATIONS AND ANALYSIS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX REFERENCES

Fluorescence anisotropy decay r(t) in optically heterogeneous molecules is defined as follows (Rigler and Ehrenberg, 1973).

(1)

where g is an optical correction factor and

(2)

are the parallel and perpendicular polarized components of the intensity decay consisting of contributions of L distinct fluorescent sites  that give rise to L anisotropies r, each modeled by a polyexponential function (see, e.g., Szabo, 1984; Döring et al., 1997):

(3)

Here _j are rotational relaxation times, _j are the corresponding preexponential factors, and r is the residual anisotropy. The total fluorescence intensity for a given fluorescent site  is modeled by a polyexponential function:

(4)

where _k are lifetimes and A_k the corresponding amplitudes. This heterogeneity may reflect the effect of N conformers that, however, are not heterogeneous with respect to anisotropy decay, or it may reflect decay of an excited donor by a donor-specific relaxation process and by energy transfer to surrounding acceptors (Bajzer and Prendergast, 1993).

The intensity components I and I can be expressed in terms of the anisotropy r and the total intensity I (see Eqs. 2-4):

(5)

In time-correlated photon counting measurement the intensity decay components I and I are related to the actual number of counts in channel i by convolution with the instrument response function denoted by R(x) (see O'Connor and Phillips, 1984):

(6)

where the subscript s denotes the  or  component,  is the zero-time shift, _s is the light-scattering correction parameter, b is the constant background, and

(7)

with h being the channel width.

Equations 3-5 imply that I and I are given by polyexponential functions with LN(M + 1) exponential components. Consequently, the procedure for discretization of the convolution integrals in Eq. 6 is identical to that used for the analysis of the total intensity decay. The previously developed discretization method (Bajzer et al., 1995) is modified in the present paper to better account for discretization errors in the first few channels, where they are most significant. In essence, we have now rigorously taken into account the fact that R(0) = 0 and that R(t) is a continuous function defined to be zero for t  0 and positive for t > 0. The two conditions cannot be satisfied for the LA approximation simultaneously, and if only one is satisfied there is no significant decrease of discretization errors. For QA and CA both conditions can be satisfied and the corresponding reduction in the discretization error significantly improves the accuracy of parameter estimation. Technical details of the modifications are presented in the Appendix.

	DATA SIMULATIONS AND ANALYSIS

TOP ABSTRACT INTRODUCTION THE MODEL DATA SIMULATIONS AND ANALYSIS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX REFERENCES

The relevance of the various discretizations with respect to the accuracy of parameter recovery for fluorescence anisotropy was tested using simulated data for the parallel and perpendicular counts F and F. A synthetic analytical IRF function of the form R(t) = t²e^t was used (see O'Connor and Phillips, 1984). The parameter  was chosen to yield a given full width at half-maximum (FWHM), and  was evaluated from the condition max_iR_i = C, where C is the number of counts in the peak channel for F. We have simulated experimental conditions commonly found in our laboratory: FWHM = 0.050 ns, C = 2 × 10⁴ counts, channel width h = 0.010 ns, and the number of channels n = 2048. The average background per channel was estimated from the count average of 10 channels in which the number of counts was simulated by a constant b corresponding to two counts per channel, subsequently corrupted with Poisson noise. We assume that scatter parameter is small and neglect the difference between and , i.e., we assume  =  = .

Counts F and F were first analytically calculated from Eqs. 6 and 7, by using Eqs. 3-5 and the synthetic analytical IRF. The scaling was chosen to yield C counts in the peak channels. To obtain synthetic noisy data, Poisson noise was added to F and F yielding the "measured" counts C and C, respectively. These data were analyzed by the maximum likelihood method (Bajzer et al., 1991; Bajzer and Prendergast, 1992), which minimizes the Poisson deviance that includes both parallel and perpendicular components:

(8)

Here the expressions for F and F are evaluated from the discretized forms of Eq. 6 as given in the Appendix (see also Bajzer et al., 1995) and as implemented in the GLOBALS software package (Beechem and Gratton, 1988; Beechem et al., 1990) in which Grinvald-Steinberg discretization is used. The IRF counts R_i involved in the discretized forms of Eq. 6 are determined from the synthetic R(t) by using Eq. 7 followed by corruption with Poisson noise _i: _i =  R(t)dt + _i.

Minimization of the Poisson deviance is performed with respect to the model parameters: _k, A_k, _j, _j, r, g, , . A novel strategy for minimization was used; it is based on a combination of n-dimensional and specially designed two-dimensional minimizations performed by using the Nelder-Mead simplex algorithm (Walters et al., 1991; Press et al., 1992), which does not require derivatives (see SIMPS by Bajzer and Penzar 1998 at www.mathworks.com/ftp/optimsv5.shtml). When the "minimum" is obtained, it is verified using a modified Levenberg-Marquardt algorithm implemented in a finely tuned subroutine that calculates derivatives numerically (Morris, 1981). The program is coded in FORTRAN 77 and all calculations were performed in double precision.

The goal of anisotropy decay data analysis is the recovery of the decay parameters: lifetimes (_k), amplitudes (A_k), rotational relaxation times (_j), and the preexponential factors _j. Consequently, the success of a given discretization scheme used in data analysis should be measured by the accuracy achieved in determining those parameters. The accuracy was judged by the coefficient of variation: CV(p) = SD(p)/ and the relative bias, RB(p) = |p  |/p for each decay parameter p. Mean values, , of the decay parameters and the corresponding standard deviations SD(p) were estimated from the values of the parameters obtained by analysis of 51 synthetic data sets that were generated with the same model parameters and differed only by the Poisson noise realization. The relative bias should be specifically sensitive to the accuracy of discretization because, in general, bias reflects deviations from the ideal model, which in our case exist in principle because of discretization. When parameters are recovered from measured data, the coefficient of variation can be estimated (either from the corresponding Hessian matrix or by Monte Carlo simulations). However, the relative bias cannot be estimated because we do not know the true values of the parameters, but if we assume that the recovered values are the true values, then those can be used to generate a number of simulated data sets as suggested above. These sets should be then analyzed, and the relative bias determined. If for a given parameter the relative bias is high, it is likely that this parameter has not been accurately determined from real data. Such an analysis is certainly more informative than just judging the value of deviance (theoretically distributed as ² distribution) and estimating the coefficient of variation, which is usually done.

In our effort to facilitate the discrimination of the heterogeneous from the homogeneous model, we propose to use information criteria for model selection defined in general as:

(9)

(cf. Davidian and Gallant, 1993; Walter and Pronzato, 1997). Here the first term represents the maximum of the log-likelihood function L(n, p, ) for n data points and the model with  free parameters, achieved for the best-fit parameter vector p = p_m. This term describes how well the given model describes the data. The second term is the penalty term for the number of involved free parameters. It takes into account that the model with more free parameters is not necessarily the more adequate model. Various weighting factors w(n) have been derived in the literature. For the Akaike information criterion (AIC) w(n) = 1, for the Bayesian information criterion (BIC) w(n) = (log n)/2, and for the Hannan-Quinn (HQIC) information criterion w(n) = log(log n) (cf. Davidian and Gallant, 1993). It can be shown that in the case of Poisson distribution, a log-likelihood function is related to Poisson deviance: L(n, p) = D/2 + C_n, where C_n for a given set of data points does not depend on free model parameters or their number. Therefore, C_n can be neglected in application of selection criteria that require finding the maximum of the log-likelihood function (or a minimum of the deviance). The model for which an information criterion achieves a lower value should be considered the preferred model by that criterion.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION THE MODEL DATA SIMULATIONS AND ANALYSIS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX REFERENCES

We have performed a great number of simulations and here we present results that illustrate the main findings. First, we pay attention to the homogeneous model, which is overwhelmingly used in the literature to interpret fluorescence anisotropy decay data. The question always remains whether in some cases it would have been more appropriate to interpret data in terms of a more general heterogeneous model. Unfortunately, such an investigation has been rarely conducted. Most authors seem to be satisfied with how data appear to be fitted by a proposed homogeneous decay. Second, we address the question of how various discretizations influence the accuracy of the recovered parameters in the heterogeneous case and discuss to what extent the heterogeneity is detectable by use of the maximum likelihood method in conjunction with information criteria for model selection.

Homogeneous anisotropy decay

In a number of simulations with GS and LA discretizations, we have found that LA discretization most often yields less accurate results for short rotational relaxation times than the GS discretization when the zero-time shift is not zero. This is illustrated in Table 1, which shows the relative biases (in percent) for the two shorter relaxation times ₁ and ₂ and for various zero-time shifts. (In this section on homogeneous anisotropy decay the notation is simplified in the following way: _{k  k1}, A_k  A_k1, _{j  j1}, _{j  j1}, r  r₁.) The results for QA and CA discretizations are also shown in Table 1 to emphasize their advantage over GS and LA discretizations.

The results of Table 1 suggest that for negligible , LA discretization is more accurate than the GS discretization method. This was also concluded by Periasamy (1988), who conducted simulations in the much simpler case of fluorescence lifetime estimation from total intensity data, and with  = 0. However, when   0, LA discretization is clearly less accurate than GS discretization. The discretization scheme of Wahl (1979) (WA) leads to a formula essentially equivalent to the one obtained for LA discretization (Bajzer et al., 1995); the only difference is that Wahl's approximation for the first channel is apparently more elaborate (see the Appendix). The result of Table 1 for WA and  = 0 suggests that the numerical consequence of this variation is remarkable, as it yields a significantly more accurate estimation. However, for very small zero-time shifts Wahl's discretization is less accurate than either GS or LA. For zero-time shifts closer to one channel width, both LA and WA discretizations are substantially inaccurate and unreliable. It should be noted that WA discretization is developed for  = 0 only, and that LA for   0 formally does not depend on : it varies with  only for negative zero-time shifts (Bajzer et al., 1995). In principle, one could formulate a discretization scheme based on Wahl's approach, which would include   0. Such an extension can be proposed in a number of ways; to find out which is the most advantageous would require a separate study, which may not be warranted as the quadratic and cubic approximations are likely to yield a more accurate parameter estimation.

In Table 1 only the relative bias for rotational relaxation times ₁ and ₂ are shown. Accuracy in the estimation of other decay parameters is also important for interpretation of data. For all discretizations except LA and WA, the relative bias for the other decay parameters was in an acceptable range of 15% (an exception was the case of GS for  = 0.008 ns with RB(₁) = 30%). As noted before, the accuracy of parameter estimation is also judged by the coefficient of variance. For all parameters and discretizations except LA and WA for  = 0.008, 0.010, CV was smaller than 15%. The results obtained clearly suggest that LA and WA are inferior to GS, QA, and CA discretizations. Therefore, in the following studies we will compare only the latter three discretizations.

An enormous number of simulations would be required to systematically cover all the relevant values of parameters and all combinations of the number of lifetimes and rotational relaxation times. A more modest, but achievable goal is to perform simulations with the sample of parameters found in actual analyses of fluorescence anisotropy decay data. In Table 2 we list decay parameters from a sample of literature data. As we could not always find the experimental conditions for the chosen literature data, we performed simulations under the experimental conditions commonly found in our laboratory (as specified in the previous section) and with r = 0,  = 0.05,  = 0.008 ns. Of course, with these arbitrary chosen conditions, the accuracy in parameter estimation may not reflect the accuracy achieved in the actual analyses of the original experimental data analyzed in the references in Table 2. The results of our simulations are presented in Fig. 1. RB and CV for only those examples and parameters are shown for which at least one of these statistics is >15% for any of the three discretizations compared. Clearly, except for a single outstanding parameter (A₃ in example (f)) QA and CA outperform GS in all decay parameters, for some very dramatically. This is also mostly true for parameters not shown, and for examples (c) and (g); however, in these cases accuracy in estimation is reasonably satisfactory (often far below 15% in CV and RB) for all discretizations. It is noteworthy that in the examples for which all parameters were estimated sufficiently accurately there was either just one long rotational relaxation time component (example (c)) or two components (example (g)) with a shorter component much longer than short components of any other example. This confirms our experience that the discretization method used will matter most when the anisotropy decay is characterized by at least two components, with at least one being short (i.e., of the order of 100 ps).

View larger version (16K): [in this window] [in a new window]   FIGURE 1   Accuracy of estimation of decay parameters as measured by the relative bias and the coefficient of variation for discretizations GS, QA, CA. The values of the parameters are chosen according to Table 2. Note that when the value of RB or CV is higher than 100%, the actual percentage is shown in parentheses. For other details see the text.

In the next set of simulations we sought the shortest possible rotational relaxation time that can be recovered under the requirement of minimal accuracy for lifetimes, amplitudes, rotational relaxation times, and the corresponding preexponential factors. The fluorescence lifetimes were chosen to be in the nanosecond range. The longer rotational relaxation time (₂) was also chosen to be in the nanosecond range, while the value of the shorter rotational relaxation time ₁ was lowered in our simulations until the relative bias or the coefficient of variation for any of the decay parameters reached a value between 14 and 16%. Knowing the decay parameters within 15% uncertainty is usually sufficient for interpretation.

In Fig. 2 A a simpler case is shown, namely when both the residual anisotropy and the scattering correction are negligible; consequently these parameters were fixed to zero and were not included as variables in the minimization. In this case, the standard GS discretization is quite stable with respect to the effects of zero-time shift change, yielding the shortest recoverable '₁ value of ~20 ps. Quadratic discretization yields shorter '₁ values for moderate zero-time shifts, but much longer '₁ values for positive zero-time shifts of ~2.5 channel widths and somewhat longer for the same negative zero-time shift. Cubic discretization clearly provides the best results with a stable '₁ value of ~14 ps.

View larger version (33K): [in this window] [in a new window]   FIGURE 2   Comparison of discretizations for anisotropy decay with a short rotational relaxation time and lifetimes in nanosecond range. Displayed are the minimal values '1 of the shorter rotational relaxation time for which the coefficient of variation and the relative bias for all decay parameters did not exceed 16%, and at least for one parameter was not <14%. Simulations were based on the following decay parameters: 1 = 1 ns; 2 = 3.7 ns; A1:A2 = 1, 1 in nanoseconds, displayed; 2 = 1.5 ns; 1 = 0.1; 2 = 0.25. The zero-time shift /h is expressed in fractions of channel width h = 0.01 ns. A and B differ in values of residual anisotropy and scattering correction parameter. The missing bar in B for /h = 2.5 indicates that the 1 value, which would meet the required criteria of accuracy in estimation, was not found.

When the residual anisotropy is not negligible and the scattering correction parameter is  = 0.05 (Fig. 2 B), all discretizations yield results for '₁ that vary significantly with the zero-time shift. Cubic discretization is again the most favorable. Application of GS discretization results in considerably longer '₁ values for smaller zero-time shifts, which probably reflects the effect of correlation among the parameters , r, and . A similar behavior is also observed for cubic discretization, while quadratic discretization yields poor results only for positive , especially for /h = 2.5, where no '₁ values that met the required criteria of accuracy were found. In general, the results of Fig. 2 suggest that to recover very short rotational relaxation times, a discretization scheme based on the cubic approximation is best.

In our third set of simulations, we chose a more complex decay with three lifetimes and three rotational relaxation times. The shortest lifetime was varied from 20 to 200 ps, while the other two lifetimes were in the nanosecond range. The two longest rotational relaxation times were also in the nanosecond range. The rotational relaxation time ₁ was chosen to be successively shorter until the required accuracy was achieved, as in the simulations above. The rotational relaxation lifetime '₁ obtained in this search is presented in Fig. 3 for various ₁, zero-time shifts, and discretization schemes. Also, for each set of parameters two cases were considered, namely when the lifetimes with their corresponding amplitudes are fixed and when they are free. The first case occurs when the fluorescence lifetimes and amplitudes have been determined from the total fluorescence intensity decay in a separate analysis.

View larger version (21K): [in this window] [in a new window]   FIGURE 3   Comparison of discretizations for anisotropy decay with a short rotational relaxation time and short lifetimes. Displayed are the minimal values '1 of the shortest rotational relaxation time found, as in Table 2. Simulations were based on the following decay parameters: 1 in nanoseconds, displayed; 2 = 1 ns; 3 = 5 ns; A1:A2:A3 = 10:8:20, 1 in nanoseconds, displayed; 2 = 0.4 ns; 3 = 3 ns; 1 = 3 = 1; 2 = 0.15. Zero-time shift as in Table 2, r =  = 0. Missing bars designate the situation when that value of 1 which would meet the required criteria for accuracy was not found. For the meaning of "fixed" and "free" lifetimes, see text.

An inspection of Fig. 3 reveals several interesting findings. 1) Using predetermined lifetimes and amplitudes is for GS discretization in some cases preferable (short ₁) and in some not (longer ₁ and /h  0). For the cubic and quadratic approximation the results are in most cases less sensitive to whether lifetimes and amplitudes are predetermined or not. Yet, for negative zero-time shifts, fixing the lifetimes and amplitudes results in shorter '₁ values when cubic discretization is applied. 2) For the shortest lifetime ₁ = 0.020 ns and all parameters free, GS discretization yielded estimates of parameters that were in all cases less accurately determined than allowed. The cubic approximation was distinctly better than the GS discretization, as it provided reasonably short '₁ for   0. The quadratic discretization was also better then GS discretization for /h = 0, 0.5. 3) Irrespective of the magnitude of the lifetime ₁, the required accuracy was not achieved for the GS discretization when all parameters were free and /h = 1.5. For /h = 0.5, the required accuracy was achieved only for the longest lifetime ₁ considered.

The data from Fig. 3 suggest, in general, that short rotational relaxation times and short lifetimes (in the range of tens of picoseconds) are most reliably recovered by using the cubic discretization. Quadratic approximation is preferable to the standard GS scheme. From the results obtained there is evidence that one can recover shorter rotational relaxation times when the lifetimes with the corresponding amplitudes are predetermined in a separate analysis. This is especially true for GS discretization, while for QA and CA this seems important only when one of the lifetimes is very short (e.g., 20 ps).

Heterogeneous anisotropy decay

When the anisotropy decay is heterogeneous the data analysis, in general, becomes more delicate, as there are more parameters to be determined and their correlation could be stronger. We have considered the fairly complex case where each of two fluorescent sites (L = 2) is characterized by a single lifetime component (N = 1) and two rotational relaxation components (M = 2) and two different residual anisotropies. This is essentially a generalization of the homogeneous decay exemplified in Table 1, examples (a) and (g), and in Fig. 2 B.

In an analogy with simulations yielding results of Fig. 2, we searched for the shortest rotational relaxation time '₁₁ that can be recovered under the requirement of minimal accuracy for decay parameters. Our simulations clearly indicated that RB and CV for most of the decay parameters were higher than for the homogeneous case. Therefore, we raised the uncertainty tolerance, i.e., we lowered the value of the shortest rotational relaxation time until RB or CV for any decay parameter reached a value between 21 and 24%. The results are presented in Fig. 4. GS discretization for /h = 0.5, 1.5 failed to yield sufficiently accurate parameter estimates for any ₁₁, and all discretization failed for /h = 0.5. Quite unexpectedly, QA and CA failed for /h = 2.5, while with GS discretization a 40 ps rotational relaxation time can be recovered within the specified accuracy. However, in general the results of Fig. 4 indicate the advantage of cubic and quadratic discretizations.

View larger version (16K): [in this window] [in a new window]   FIGURE 4   Comparison of discretizations for heterogeneous anisotropy decay with a short rotational relaxation time and lifetimes in the nanosecond range. Displayed are the minimal values '11 of the shortest rotational relaxation time for which the coefficient of variation and the relative bias for all decay parameters did not exceed 24%, and at least for one parameter was not <21%. Simulations were based on the following decay parameters: 11 = 1 ns; 12 = 3.7 ns; A11:A12 = 1, 11 in nanoseconds; 21 = 1.5 ns; 12 = 0.5 ns; 22 = 2.5 ns; 11 = 0.1; 21 = 0.25; 12 = 0.2; 22 = 0.1; r1 = 0.02; r2 = 0.03. The zero-time shift /h is expressed in fractions of channel width h = 0.01 ns,  = 0.05. Missing bars correspond to the situation where the value of 11 for which the required criteria of accuracy in estimation would be met, was not found.

The issue of discretization appears to be important in discriminating between the heterogeneous and homogeneous cases. We have synthesized data for the heterogeneous model and analyzed them, once implying the same heterogeneous model and then the corresponding homogeneous model (L = 1, N = 2, M = 2). For each of 51 simulated analyses for a given set of decay parameters, we calculated model selection criteria AIC, BIC, and HQIC and determined the probability that a given criterion falsely indicates the homogeneous model as preferable. In 15 different parameter sets, each analyzed by using GS, QA, and CA, we found that AIC was always a better model selection criterion than BIC or HQIC, i.e., the probability P_AIC that AIC would indicate a homogeneous model was always smaller than or equal to the corresponding probabilities for BIC and HQIC. In fact, even in cases where heterogeneity was rather marginal (e.g., ₁₁ = 0.2, ₂₁ = 0.25, ₁₂ = 2, ₂₂ = 2.5) P_AIC was not larger than 0.34, while P_BIC was typically above 0.5. In all but one case P_AIC was smaller for QA and CA discretization than for GS discretization, clearly indicating the advantage of quadratic and cubic discretizations in discriminating the heterogeneous from the homogeneous model (see examples in Table 3, last column).

Even when the Akaike model selection criterion faithfully justifies the heterogeneous model, the question remains whether in fact the two rotational relaxation times can be considered different when uncertainties in parameter estimation are taken into account. For that reason we used a t-test to find out when it is justifiable to consider two close rotational relaxation times statistically different. Illustrative examples presented in Table 3 suggest that CA discretization is the best in statistically discriminating different corresponding rotational relaxation times. In evaluating how close these times are we used the ratio measure defined as R = ₂₁/₁₁ + ₂₂/₁₂. If there is no heterogeneity R = 2, and for the heterogeneous case R > 2. The ratio measure for the means of recovered parameters, R', is much closer to R for QA and especially CA than for GS. This again provides evidence that cubic discretization is advantageous. The same conclusion can be obtained when considering the median of the average relative error for anisotropy parameters defined in Table 3.

	CONCLUDING REMARKS

TOP ABSTRACT INTRODUCTION THE MODEL DATA SIMULATIONS AND ANALYSIS RESULTS AND DISCUSSION CONCLUDING REMARKS APPENDIX REFERENCES

This paper has addressed the hitherto unexplored issue regarding the influence of different discretization schemes on the accuracy of estimation of anisotropy decay parameters. We have found that the improved discretization schemes presented in this paper allow reliable and accurate estimation of rotational relaxation times, even those of the order of only a few channel widths, providing that the FWHM of IRF is five channel widths. Although the cubic approximation of the instrument response count function is found most reliable, the quadratic approximation is often satisfactory. The commonly used Grinvald-Steinberg discretization is less adequate for short relaxation times and can yield estimates that are twice shorter or longer than the actual value. However, for longer rotational relaxation times (equivalent to 30 channel widths or more), GS discretization becomes more accurate and provides estimates that will not impede interpretation. In retrospect, those data analyses of complex anisotropy decays in the literature performed by GS discretization, which did not yield rotational relaxation times shorter than 30 channel widths, can be considered sufficiently accurate.

In an effort to address the difficult problem of discriminating between homogeneous and heterogeneous anisotropy decay we have found that the Akaike information criterion is more useful and reliable in model selection than Bayesian and Hannan-Quinn information criteria. It is shown that even in the case when slightly heterogeneous decay appears homogeneous due to statistical uncertainty in parameters, the Akaike criterion is capable of detecting heterogeneity. We found that the discretization scheme also matters in discriminating between homogeneous and heterogeneous anisotropy decay. Again, when cubic discretization was applied it was possible to detect heterogeneous anisotropy decay more often than for quadratic and even more so for GS discretization.

The results presented here were obtained for a fixed channel width. However, these results are still valid if all time variables and parameters are scaled to be in the same proportion with some other channel width.

Based on the results obtained and our experience from numerous simulations, we recommend the following general steps for reliable estimation of fluorescence anisotropy decay parameters.

1.	Use the software that has discretizations based on quadratic and cubic approximations of the IRF count function;
2.	First estimate lifetimes and amplitudes from the total fluorescence intensity decay data. Use cubic discretization. Estimate coefficients of variation by using Monte Carlo simulations;
3.	Assuming that obtained lifetimes and amplitudes are fixed parameters, estimate the rotational relaxation times with the corresponding preexponential factors for the homogeneous anisotropy decay. Use cubic discretization;
4.	Perform an equivalent analysis to 3) for various heterogeneous models. By applying the Akaike information criterion, determine whether any heterogeneous model is more adequate than the homogeneous model. Methods proposed by Bialik et al. (1998) can be further used to discern the most adequate decay model;
5.	Estimate the coefficients of variation in anisotropy decay parameters for the most preferred model, using the Monte Carlo simulation;
6.	Estimate anisotropy decay parameters for the preferred model by using the quadratic discretization. If obtained parameters differ from those obtained by the cubic discretization so much that the interpretation becomes questionable, perform simulations along the lines suggested in this paper and, based on the estimated accuracy in parameters (indicated by the coefficient of variation and the relative bias), infer which of the two discretizations is more adequate;
7.	As a final check one may perform yet another fit to data in which all decay parameters are free, using the most preferred model and the most adequate discretization.