Implementation and Evaluation of Data Analysis Strategies for Time

---

Article pubs.acs.org/ac

Implementation and Evaluation of Data Analysis Strategies for TimeResolved Optical Spectroscopy Chavdar Slavov,* Helvi Hartmann,† and Josef Wachtveitl* Institute of Physical and Theoretical Chemistry, Goethe University, Frankfurt am Main, Germany S Supporting Information *

ABSTRACT: Time-resolved optical spectroscopy plays a key role in illuminating the mechanisms of many fundamental processes in physics, chemistry, and biology. However, to extract the essential information from the highly complex timeresolved data, advanced data analysis techniques are required. Here we present the implementation strategies and the evaluation of the familiar global lifetime and target analysis as well as the not so widely adopted lifetime distribution analysis (LDA). Furthermore, we demonstrate the implementation of analysis strategies dealing with a number of artifacts inherently present in data from ultrafast optical experiments. The focus of the work is placed on LDA as it allows invaluable exploration depth of the kinetic information contained in the experimental data. We establish a clear regularization procedure for the use of LDA in ultrafast optical spectroscopy and evaluate the performance of a number of factors that play a role in the reliable reconstruction of lifetime distributions. Our results show that the optimal regularization factor can be determined well with the L-curve and the generalized cross-validation techniques. Moreover, the performance evaluations indicate that the most efficient regularization norm is the identity matrix. The analytical procedures described in this work can be readily implemented and used for the analysis of any time-resolved data.

T

Typically time-resolved data is subjected to global analysis, which is a simultaneous analysis of multiple kinetic traces recorded in dependence of wavelength or other experimental variables.11,16 The most common global analysis procedure is global lifetime analysis (GLA), where the data, S(t), is approximated by a discrete sum-of-exponentials function (eq 1). The global parameters spanning the dataset are the lifetimes, τ, while the pre-exponential amplitudes, A, are determined for each kinetic trace. GLA results in the so-called decay-associated spectra (DAS), where the pre-exponential amplitudes for each lifetime component are plotted against an experimental variable, e.g., detection wavelength, λi. The DAS is a compact representation of the kinetic information in the data. GLA can be viewed as a special case of the more powerful global target analysis (GTA), where kinetic models are directly tested for compliance with the data. GTA is essential because kinetic modeling is the only quantitative way to reveal reaction mechanisms. Even the model-independent analysis approaches, which are powerful in determining the number of components involved in the reaction, resort to kinetic modeling to determine the reaction mechanism.14 Regrettably, in a significant number of ultrafast studies, the data analysis stops with the GLA and at most a qualitative reaction scheme based on the DAS interpretation is

ime-resolved optical spectroscopy plays a key role in investigating the mechanisms of a large variety of photophysical and photochemical reactions and of numerous biological processes.1−9 Together with the improvement of the experimental techniques, essential part of the progress in the field is the development of procedures that allow in-depth analysis of the acquired data. Excellent examples of such advanced tools are “TIMP/Glotaran”,10,11 “Ultrafast toolbox”,12 and “MCRALS”.13,14 Despite some specific differences that should be accounted for, e.g., noise statistics, presence of experimental artifacts, etc., the basic analysis strategies for time-resolved data are generally the same for all data acquisition techniques. Those strategies can be classified as model-independent and modeldependent.14 In model-independent strategies, the decomposition of the experimental data into the time-dependent profiles of the different reacting species is performed using factor analysis methods.15 The decomposition is done without a priori assuming an explicit mathematical function. On the contrary, in modeldependent strategies always a specific reaction mechanism expressed in terms of a mathematical function is used.11 The data analysis program presented here uses predominantly modeldependent and related approaches that will be further summarized (see ref 14 for a detailed description of the modelindependent strategies in ultrafast spectroscopy). n

S(t , λexc , λi) =

Received: November 6, 2014 Accepted: January 15, 2015 Published: January 15, 2015

∑ Aj(τ , λexc , λi) exp(−t /τj) (1)

j=1 © 2015 American Chemical Society

2328

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry

the integral in eq 3 needs to be discretized into a quasicontinuous sum of n exponential functions similar to eq 1 but with n typically >50. This discretization does not simply reduce the method to GLA, where a priori a discrete sum of exponentials is assumed, but still allows recovery of lifetime distributions. In effect, the LDA approach becomes model-independent and naturally takes care of complicated analysis problems including nonexponential or stretched exponential kinetics. The advantage of LDA becomes apparent when it is compared with GLA for the description of complex kinetics.28,30 In contrast to GLA, which typically suggests a number of similar solutions with DAS relatively complex for interpretation, LDA compresses the kinetic information into a distribution map, which gives a comprehensive overview of the kinetics. Here we present the development of OPTIMUS, a new, versatile data analysis tool for time-resolved optical spectroscopy. Implemented here are all the methods discussed above, i.e., GLA, GTA, and LDA, as well as procedures for analysis of specific ultrafast spectroscopy features, like wavelength-dependence of time zero11 and coherent artifacts.31 The focus of the work is placed on the LDA method applied to ultrafast spectroscopy. To demonstrate the capabilities of the method, we have conducted a detailed analysis of its performance in dependence of a variety of factors. The motivation for developing OPTIMUS was to conceive a tool that is capable of promptly executing a variety of analysis strategies and thus allowing the researcher to focus on exploring different hypotheses rather than to struggle with analysis difficulties of technical nature. The module OPTIMUSGLA will be available as a standalone application with an intuitive graphical user interface (GUI) for free download at ref 32. At present, GTA and LDA are command line based; however, dedicated GUIs are planned. The capabilities of the LDA module were demonstrated in a recent study on the complex excited state dynamics of a pyrene-adenosine derivative.30

proposed. However, because of the often very high complexity of the studied reactions, such an approach is relatively speculative and rather the proper procedure is to test actual kinetic models directly for data conformity. Although, the latter step is complicated to accomplish, the methodology is generally established.10−12,16−18 In essence, GTA consists of devising a kinetic reaction scheme based on previous knowledge and certain assumptions. The kinetic scheme is then represented by a system of differential equations. Here the discussion will be limited to first-order kinetics, where the time-dependent populations of the involved species are described by homogeneous first order differential equations (X(t) is the vector of the time-dependent population of the involved species, T is the kinetic transfer matrix):

d X(t ) = TX(t ) dt

(2)

The solution of this system of differential equations is a sumof-exponentials function like in eq 1, with the difference that now the lifetimes, τ, depend in a complex manner, defined by the kinetic model, on the kinetic rates. It is now apparent why the GLA is a special case of GTA, as it represents the solution of a system of linear independent differential equations, which describe a kinetic scheme with parallel, noninteracting states with equal initial concentrations. The major challenge with kinetic modeling is the identification of the correct model.11 GTA’s nature is such that many models can potentially fit the data, moreover, even for the same model a great number of feasible solutions are possible. Thus, additional constraints are needed to identify the most reasonable model. Such constrains are based on previous knowledge, e.g., spectral properties of the components, rates for particular reactions, etc. In addition, the kinetic model of choice should not violate any of the known physical laws. A kinetic model that describes well the data and fulfills the aforementioned conditions can be considered right only until proven otherwise. Despite the identifiability challenge, kinetic modeling remains the only method that can reveal reaction mechanisms. Recently, another powerful approach, lifetime density/ distribution analysis (LDA), gained popularity in fields where complex kinetics on a microsecond or longer timescale is studied, e.g., protein folding,19−22 rhodopsin cycle,23 etc. Although, the LDA method uses exponential functions for describing the data, it lies at the interface of model-dependent and modelindependent methods. Despite the general availability of LDA algorithms24−26 established more than 20 years ago, currently the approach is not massively adopted in ultrafast spectroscopy and to our knowledge only the Holzwarth group has used it extensively.27,28 This is surprising since the complexity of the studied samples and processes is generally beyond the discrete sum-of-exponentials frame (GLA). The great number of spatial, energetic, and temporal degrees of freedom of the studied molecule ensembles and their matrix produces a continuous distribution of individual exponential decays (eq 3). In this respect, the approximation of the data profile by a discrete set of exponentials (GLA) must be interpreted as a curve parametrization, concealing potential contribution of a larger number of underlying decays.26,29 S(t , λexc , λi) =

∫0

■

MATERIALS AND METHODS OPTIMUS was coded using the MATLAB programming language (Mathworks, Inc.) and takes advantage of the available toolboxes, Optimization and Global optimization, Statistics, Parallel computing, GUIDE (GUI building), and MATLAB compiler. The program was tested on a standard 32-bit, Windows 7 PC (Intel Core i3-530 CPU, 4 GB RAM). The structure of OPTIMUS is modular with a common computation routine for all analysis methods (Figure 1) to allow future extensions without significant alterations of the core. The fitting of the data is realized by minimizing the square of the residual norm (eq 4, Sexp(t), experimental data; Sfit(t), fit), using the nonlinear least-squares MATLAB function lsqnonlin and employing the trust-region-reflective algorithm. f (t ) = min

2 2

(4)

Often for such a multidimensional problem, the solution surface contains several local minima, in which a gradient based algorithm could be trapped. To avoid such a scenario and to find the global minimum, one needs to explore the solution space by starting the optimization from different points. Such an option is implemented here based on the algorithms provided by the Global optimization toolbox. Typically the multistart function is used, which generates a specified number of starting sets within the given boundary conditions and automatically performs the nonlinear least-squares optimization for each starting set. Furthermore, OPTIMUS employs the Parallel computing

∞

Φ(τ , λexc , λi) exp( −t /τ ) dτ

Sexp(t ) − Sfit(t )

(3)

The function S(t) represents the Laplace transform of the spectral distribution function, Φ(τ). For computational reasons, 2329

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry

S3−S5 in the Supporting Information)31,33 where the CA is approximated by a function composed of a Gaussian and/or its first and second derivative. Global Lifetime Analysis. Since GLA is a special case of GTA, it generally does not need a separate implementation. However, for simplicity reasons and computational speed, OPTIMUS has a separate GLA module. In GLA, a model function for each detection wavelength is constructed by combining a sum of convoluted with the IRF exponential functions (eq S1 in the Supporting Information) with, if desired, the functions accounting for the CA (eqs S3−S5 in the Supporting Information). The nonlinear parameters that are optimized in the gradient search are the lifetimes, the IRF widths (varied independently for each wavelength, approximated by a polynomial function, or fixed if known) and the time zero position (eq S2 in the Supporting Information). The optimization is initiated with some starting values, used to calculate each of the basic functions (eqs S1−S5 in the Supporting Information). The results are then used to form the columns of a guess matrix, Sguess. The contribution (amplitudes) of each function to the experimental data, sexp, is found by linear fitting using the Moore-Penrose pseudoinverse34,35 (MATLAB pinv function) of Sguess,S+guess.

Figure 1. Block-diagram of the structure and operation of OPTIMUS.

+ a = Sguess sexp

toolbox to speed up the solution space exploration by distributing the computations to different CPUs. General Considerations. The analysis of ultrafast timeresolved data is often hindered by experimental artifacts (see the Supporting Information for details on their treatment) which need to be mentioned before describing the implementation of GLA, GTA, and LDA. Despite the significant advances in the experimental techniques, the time-resolution of the setups is finite and determined by the instrument response. Hence, the recorded data, Sexp, do not represent the pure molecular response to a Dirac δ-function input but rather contain also the instrument response function (IRF). Hence the detected signal, Sexp, is a convolution, ⊗, of the pure molecular response, eventually approximated by a sum-of-exponents function in GLA, GTA, and LDA and the IRF (eq 5). The reconstruction of the pure molecular response involves solving the ill-posed deconvolution operation, which is a tedious task and thus the problem is addressed by performing iterative reconvolution of the IRF with the simulated molecular response (Supporting Information, eq S1).

and sfit = Sguessa

(7)

The fit to the experimental data, sfit, is given by eq 7, and the fit parameters are optimized iteratively to minimize the square of the residual norm (eq 4). The elements of a are plotted against the detection wavelength to yield the DAS. Global Target Analysis. The implementation of GTA follows the standard linear system theory approach.18 Typically the envisioned kinetic models are expressed as a system of firstorder differential equations like eq 2 but extended to take into account the IRF, IRF(t), and the initial concentrations vector, ε (ε(λexc), if the initial concentrations depend on the excitation wavelength):36 d X(t ) = TX(t ) + ε(λexc)IRF(t ) dt

(8)

The kinetic transfer matrix T is formed by n

Tjk = kjk − δjk(kjj +

n

S(t , λexc , λi) =

(6)

∑ kmk), m=1

∑ Aj(τj , λexc , λi) exp(−t /τj) ⊗ IRF(t )

j , k = 1, .... , n (9)

where δjk is the Kronecker delta. The off-diagonal elements of the transfer matrix, Tjk, represent the transition rates (kjk) between the compartments, while the diagonal elements, Tjj, describe the sum of all decay rates of a compartment. Such a system of differential equations can be solved analytically (eq 10). The eigendecomposition of the transfer matrix yields T = UΛU−1, and hence exp(Tt) = U exp(Λt)U−1 (Λ, diagonal matrix with the eigenvalues; U, eigenvector matrix of the transfer matrix, T). The eigenvalues are related to the measured lifetimes via Λjj = −1/τj, while the eigenvector matrix indicates the contribution of a particular compartment to the lifetime of a process. A negative value in this matrix corresponds to the rise of a compartment’s population, while a positive value represents the depopulation of a compartment.

j=1

(5)

Wavelength dependence (chirp) of the time-zero position of the recorded signals, caused by group velocity dispersion and self-phase modulation of light, is another experimental artifact present in picosecond experiments. The chirp can be modeled by a polynomial function11 (Supporting Information, eq S2). The employment of ultrashort excitation pulses (∼100 fs or shorter) in transient absorption spectroscopy often leads to distortions in the recorded signal at time-zero position called coherent artifact (CA). The CA is caused by stimulated Raman scattering and/or by pump-induced temporal changes in the medium refractive index, which lead to spectral changes in the weak probe via cross-phase modulation.31 An elegant approach to treat the CA contributions is introduced in OPTIMUS (eqs

X(t ) = exp(Tt ) ⊗ ε(λexc)IRF(t ) 2330

(10) DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry

balance between the regularization and the residual norms. The TR problem can be formulated similarly to eq 6:

Finally, an expression for the time-dependence of the experimental signal can be derived (SASk, kth species spectrum):

+ a = Sguess , αsexp

S(t , λexc , λi) n

= (∑ SASk (λi)(U −1ε(λexc))j Ujk) ∑ exp(Λjjt ) ⊗ IRF(t ) k=1

where S+guess,α is the Tikhonov regularized inverse:

j=1

+ T 2 T −1 T Sguess , α = (SguessSguess + α L L) Sguess

(11)

It is apparent that the expression is similar to eq 5. However, in GTA the lifetimes and the amplitudes of the model function are no longer pure mathematical parameters but are now related to physical parameters, like rate constants of particular processes and spectra of particular compartments, respectively. An important parameter that emerges from the kinetic description is the species-associated spectrum (SAS). Physically SAS represent the stationary spectra of the model compartments as if they were measured separately. SAS are related to DAS through the first summation in eq 11, which is equal to the preexponential amplitude in eq 5. In effect, at a given detection wavelength, the DAS of a particular lifetime component is a linear combination of the SAS of the compartments.

∑ bjkSASk (λi) k=1

(12)

In OPTIMUS, the user needs to set up a kinetic rate matrix and an excitation vector (initial concentrations vector) describing the envisioned model along with the parameters for fitting the different experimental artifacts. The fitting of the experimental data is carried out similarly to GLA with the difference that now the exponentials are calculated using the result from the eigendecomposition (eq 11). Note that in GTA the amplitudes calculated using eq 6 represent the SAS and NOT the DAS. The DAS can be calculated via eq 12. Lifetime Density/Distribution Analysis. The general idea of LDA was described in the introduction, and thus this part is concentrated on the technical details of the method implementation. Obtaining the pre-exponential function Φ(τ,λexc,λi) (eq 3) or rather its discrete (eq 1) but quasicontinuous form, Aj(τ,λexc,λi), j = 1,...,n > 50, is equivalent to performing the inverse Laplace transformation of the data, S(t,λexc,λi). However, this operation is an ill-posed problem and hence cannot be performed in a trivial manner; even minimal noise levels will hinder simple numerical procedures such as the one described for GLA from recovering the true lifetime distribution and will induce significant artificial amplitude oscillations. Instead, ill-posed problems need to be reformulated to include additional information (e.g., smoothness), which improves their conditioning so that a stable numerical solution can be obtained.37 The procedure is called regularization, with Tikhonov regularization (TR)38 and its variants being most widely used. TR is implemented in OPTIMUS-LDA for iterative inverse Laplace transformation of the experimental data to obtain the lifetime distributions. In general form37 and formulated for the problems discussed here, TR can be expressed as min{ Sguessaα − sexp

2 2

+ α 2 La α 22 }

(15)

The regularization (eq 13) simply extends the nonlinear minimization problem (eq 4) and thus all methods (GLA, LDA, and GTA) share one computational routine (Figure 1), and only the preparation of the Sguess matrix is special for each method. In LDA Sguess is calculated using a constant (not altered during the minimization), quasi-continuous set of lifetimes (typically ∼100) distributed evenly on a decimal log scale. This is necessary to account for kinetics occurring on several timescales (e.g., femtosecond to nanosecond). Therefore, the implementation of LDA and GLA differs only by the setting up of the lifetimes and the introduction of a regularization term, α∥Laα∥2. The aforesaid options (ID, FD, SD) for the operator, L, are implemented in OPTIMUS. The fitting of the experimental artifacts is as described above. Selection of the Regularization Parameter. The most challenging step in LDA is the selection of an optimal regularization parameter, α, defining the regularization strength. If α is too large, the regularized solution will not fit the data well and the residual norm will be large. On the contrary, if α is too small the residual norm will be small, but significant overfitting will occur, which will render the solution meaningless. Three methods for optimal selection of α are implemented here: Lcurve,39,40 minimal product method (MPM),41 and generalized cross-validation (GCV).42 L-Curve Criterion. This graphical method is the most common for optimal α estimation. The L-curve is a decimal log−log plot of the smoothing norm, ∥Laα∥2, versus the corresponding residual norm, ∥Sguessaα − Sexp∥2.39,40 Typically, the L-curve represents the compromise between the two quantities, and the corner of the curve appears at the optimal α. Here, the method is realized by finding the regularized solution of the LDA problem for a large number of α values spread on a log 10 scale. The calculations at each α are performed using parallel computing and, if desired, global optimization. The Lcurve is plotted at the end of the calculations. Minimal Product Method. The MPM originates from the brain research field,41 and it is relatively new and not widely used. The method is a modification of the L-curve, where the product, 7 , between the smoothing norm, ∥Laα∥2, and the residual norm, ∥Sguessaα − sexp∥2, is plotted for each α value.41 The minimum of the 7(α) curve corresponds to the optimal α value. The method could be advantageous when the corner of the L-curve is not welldefined (e.g., at low S/N ratios). In such cases, the increased curvature of 7(α), as compared to the L-curve, could be helpful in determining the optimal α. The computational strategy for this method is the same as for the L-curve and only the final graphical output is distinct. Generalized Cross-Validation. GCV42 was introduced here following its successful implementation for similar applications.22,23 GCV is a predictive method where the square of the

n

DASj (λi) =

(14)

n

(13)

where ∥Laα∥2 is the discrete smoothing norm with L being the identity matrix (ID), a discrete approximation of a derivative operator (typically first (FD) or second (SD)) or a diagonal weighting matrix.37 The regularization parameter, α, controls the

regularized residual norm, ;(α) =

Sguessaα − sexp

2

, is

2 37,43

minimized based on how well α predicts missing data values. The GCV function is defined as 2331

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry

.(α) =

Sguessaα − sexp

(n −

2

2 2 + trace(SguessSguess, α)

)

(16)

In OPTIMUS, .(α) can be searched using global optimization and parallel computing and the result is presented as a .(α) or ;(α) plot, where the minimum gives the optimal α.

■

RESULTS AND DISCUSSION The performance of the GLA, GTA, and LDA modules was evaluated using simulated datasets, resembling real experimental data. Either 500 or 700 picosecond timescales composed of a linear part (25 fs step) between −1 and 1 ps, and a logarithmically growing scale onward was used. The synthetic datasets were generated using eqs S1−S5 in the Supporting Information (convolution with IRF (80 fs), chirp, and simulated CA). Furthermore, Gaussian noise was added to test how well the initial parameters can be recovered in noisy experiments. The noise was simulated using the MATLAB function randn, by generating random values with mean 0 and variance 1 and multiplying the result by the desired standard deviation. Coherent Artifact Module. Typically, a pure solvent measurement is used to determine the CA contribution in transient absorption. In OPTIMUS this measurement is analyzed using the CA module to determine the cross-correlation width and the parameters for the chirp fitting and correction. The parameters are useful later on in the fitting of the sample dataset. Figure S1 in the Supporting Information presents the results from the fitting of a solvent measurement. The fitting was performed using variable widths for the Gaussian functions at each detection wavelength and a wavelength dispersion order of 2. The starting values were selected based on the duration and position of the signals. The boundaries were set to ∼60% of the starting value to allow the algorithm to find the optimal values. The results demonstrate excellent approximation of the CA and the chirp in the data. Global Lifetime Analysis Module. The synthetic dataset for testing this module was generated based on an artificial DAS composed of 5 lifetimes (Figure 2A). The spectral shapes of the DAS were given by Gaussian functions. The calculation of the data was as described above. The parameters for generating the chirp are indicated in Table 1. In Figure 2B, the complete synthetic data surface (Gaussian noise, std dev 8.5%) is presented. The CA can be seen in the figure as a sharp positive−negative oscillation at the very beginning of the signal. The synthetic dataset was fed into the GLA module, and values for the fit parameters were guessed. Typically, some fit parameters can be guessed closely since their potential values are either obvious from the data or were previously estimated. For example, the IRF width can be approximated by the widths of the pump and the probe pulses, or it can be estimated from a solvent measurement, the central wavelength and its time point, necessary for the chirp fitting, can also be guessed closely as they are obvious from the data. In this respect, for the fitting tests, those parameters were given good guess values (Table 1). However, for the other parameters, e.g., lifetimes, dispersion order, dispersion factors, number of CA components, the guesses were based on common sense, i.e., overall look of the chirp curvature and the kinetics; in general, following the steps when dealing with an unknown sample. The results from the fitting of the synthetic dataset are shown in

Figure 2. Evaluation of OPTIMUS-GLA. (A) Simulated DAS used to generate a synthetic dataset; (B) synthetic dataset with added CA, chirp; (C) result from the fitting of the dataset (cf. Table 1 for details); (D) DAS recovered from the synthetic dataset; (E) comparison of synthetic data and fit at a single decay channels; (F) residuals plot. The optimization of the 9 fit parameters on the time-wavelength data matrix (56 × 381) takes ∼12 s.

Table 1. Testing of OPTIMUS-GLA on a Synthetic Dataset τ1 τ2 τ3 τ4 τ5 IRF λC c0 c1, c2

original

guess

boundaries

fit

0.16 ps 6.0 ps 40.0 ps 150.0 ps 1 ns 0.08 ps 650 nm −0.1 ps 0.44, 0.1

0.5 ps 10.0 ps 50.0 ps 100.0 ps 1 ns 0.05 ps 649 nm −0.1 ps 0, 0

0.08−inf ps 1.0−inf ps 10.0−inf ps 50.0−inf ps fixed 0.03−0.2 ps 600−700 nm −0.2 − 0 ps −1, −1−1, 1

0.16 ps 6.49 ps 41.5 ps 142.9 ps 1 ns 0.079 ps 647 nm −0.115 ps 0.43, 0.1

Figure 2C−F and Table 1 (cf. also Supporting Information, Figure S2). Notably, even at these high noise levels, the algorithm implemented here is capable of accurately recovering the initial parameters (Table 1) and the DAS shapes (Figure 2D) (Please check the Supporting Information, Table S1 for the 95% confidence intervals for lifetime determination.). The last lifetime (1000 ps) was kept fixed for the presented simulation as it is longer than the used timescale. The purpose of including it was to evaluate the behavior when shorter timescales are used, as is often the case in real experiments, and to get an idea on the reliability of such estimates. As anticipated, the determination uncertainty of this parameter is relatively large (any value >900 ps can yield a reasonable fit), especially when the amplitude of the component is considered with respect to the noise level. The information loss with the decrease of the S/N ratio reminds that under real experimental conditions excessive noise could obscure important lifetime components. Global Target Analysis Module. The synthetic dataset for the GTA testing was generated using the kinetic model in Figure 2332

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry

enlarge the solution space. Hence, the solution space needs to be explored thoroughly, by starting the fitting from a variety of initial points. The implementation of global optimization and parallel processing in OPTIMUS is very helpful in this respect. Finally, the results of the kinetic modeling need to be judged based on previous knowledge to select the most reasonable model. The GTA was performed on the synthetic data using global optimization, initiated from 60 different starting values sets within the given boundaries. The parameters for the data simulation, the starting conditions for its fitting and the resulting fit values are shown in Table 2.

3A and includes chirp and CA. The dataset preparation (Gaussian noise, std dev 8.5%) was as explained above; however,

Table 2. Testing of OPTIMUS-GTA on a Synthetic Dataset

Figure 3. OPTIMUS-GTA testing: (A) Kinetic model used for simulating a dataset; (B) assignment of the kinetic rate matrix; (C) synthetic dataset based on the kinetic model in part A and the SAS in Figure 4C; (D) result from the fitting (cf. Figure S4 in the Supporting Information for the fit quality); (E) residuals; (F) fit plotted without chirp and CA. The optimization of the 10 fit parameters on the timewavelength data matrix (61 × 281) takes ∼10 s per CPU per starting value set.

the DAS was formed in a complex manner: (i) the kinetic model was converted into a transfer matrix; (ii) the transfer matrix was subjected to eigendecomposition to get the lifetimes and the eigenvector matrix; (iii) the eigenvector matrix was weighted by an excitation vector; (iv) SAS were generated using Gaussian shapes (Figure 4C); (v) DAS (Figure 4A) was calculated from the SAS and the weighted eigenvector matrix. Naturally, the fitting of the synthetic data (and of any experimental data) with a kinetic model represents a challenge due to the great number of model degrees of freedom which

The results in Figure S3 in the Supporting Information are from the unconstrained fit, while the results in Figure 3D−F and Figure 4B,D are from a fit where, based on the knowledge about the initial SAS of C3 (Figure 4C), part of the C3 SAS (610−630 nm) was set to 0 during the fitting. The result is included to exemplify the importance of incorporating previous knowledge in kinetic model testing. Although the two best fits (Figure 4B,D and Figure S3A,B in the Supporting Information) are identical and the parameter values match well (given the S/N ratio), the values used for simulating the test dataset (Table 2), constraining the C3 SAS reduces significantly the number of local minima in the solution space. For example, the local minimum solution in Figure S3C,D in the Supporting Information with a wrongly determined SAS of C3 is no longer possible in the constrained case. The results illustrate not only the importance of introducing previous knowledge in the kinetic model but also the way kinetic models are discriminated based on such knowledge. Under real experimental conditions, models that deliver physically unreasonable rates or spectra that do not agree with previous results are typically discarded. Considering the S/N level in the simulated dataset (8.5%), the rate constants recovered in the GTA that give reasonable SAS are typically within ∼10% of the original value for the unconstrained case and within ∼5% for the constrained case. Similar precision should be expected under real experimental conditions given the right model is selected.

Figure 4. Evaluation of the capabilities of OPTIMUS-GTA: (A) simulated DAS used to generate a synthetic dataset; (B) DAS resulting from the fitting of the synthetic dataset with constraint SAS of C3 in the 610−630 nm range; (C) simulated SAS; (D) SAS resulting from the fitting. The unconstrained solution is shown in Figure S3 in the Supporting Information. 2333

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry

Figure 5. Evaluation of OPTIMUS-LDA: (A) simulated LDM used for generating the synthetic datasets; (B−D) LDMs recovered from the different synthetic datasets (Figure S5A−D in the Supporting Information) using the ID regularization term. The optimal α was selected as described in the text. The optimization of the fitting parameters takes ∼100 s per CPU depending on the regularization type and the S/N level.

fit parameters is the same as for the GLA and GTA, except for the IRF width when the FD and SD regularization terms are used. Effect of the Noise Level on the LDM Recovery. Figure 5B−E shows the LDMs recovered from datasets with different noise levels (Figure S5A−E in the Supporting Information) simulated using the “ideal” LDM in Figure 5A. Intentionally, the LDM used for the simulation of the datasets was formed by lifetime distributions with different, but relatively narrow widths, and with different spacing and amplitude. Such a construction of the LDM allows detailed evaluation of the capabilities and the limitations of the LDA method. The in-depth analysis of the noise effect on the LDA shows that the increase in the data uncertainty leads to broadening of the lifetime distributions, which, depending on their distance, could eventually merge into a broader distribution. Nevertheless, even at very low S/N levels (Figure 5D,E), the LDA delivers an invaluable view of the studied kinetics, which cannot be obtained by the other methods. The results presented here compare very well with the resolution capabilities demonstrated in previous works,20−23,26,44−48 with the important difference that here the tests were performed taking into account the additional complexity present in datasets from ultrafast spectroscopy, e.g., need for IRF deconvolution, chirp, and CA correction. Evaluation of the Effect of the Regularization Operator, L. The purpose of the regularization operator is to improve the conditioning of the problem and thus help with its numerical solution. However, the choice of the regularization term is often not obvious. Henceforth, here the effect of the most common regularization operators (ID, FD, SD) on the recovery of the LDM is examined. In Figure S10 in the Supporting Information, the LDMs obtained with the different regularization terms at two different noise levels (0.14 and 2.5%) are shown. The best solution for the current application is obtained using the ID term, while FD and SD introduce uncertainties in the amplitude determination of the short (1 ns) lifetimes (0.14% and 2.5% noise) and also lead to an increase of the oscillatory effects (e.g., Figure S10 in the Supporting Information, 0.14% noise, distributions at ∼685 nm). The latter is most likely due to the tendency of FD and SD terms to sharpen the distributions. However, often the sharpening is associated with overshooting/ringing effects, which make the recovered LDMs less reliable. The ID term appears most suitable for data produced in ultrafast spectroscopy, and the LDMs obtained with this term have the highest reliability. The LDMs shown in Figure 5 are obtained using the ID term.

Lifetime Density/Distribution Analysis Module. Extensive evaluation of the OPTIMUS-LDA performance, and hence of the TR method, for application in ultrafast optical spectroscopy is conducted in this work. The generation of the synthetic dataset for LDA testing differs from the procedure for the GLA module testing only by the generation of the pre-exponential amplitudes. Since LDA is performed using a quasi-continuous lifetime distribution, both spectral and lifetime distribution shapes need to be defined. In both cases, Gaussian shapes were used to simulate a three-dimensional lifetime distribution map (LDM).27 The LDM is a quasi-continuous analog of the DAS and accordingly is used for calculating the synthetic datasets. For the tests, five Gaussian lifetime distributions were simulated, each with a characteristic spectral shape. To test the resolution capabilities of the LDA module, the simulated LDMs (Figure 5A) contain both closely and more distantly spaced distributions.

Figure 6. (A) Dataset simulated using the LDM in Figure 5A with CA, chirp, and noise (1.2%); (B) fitting result using the LDM in Figure 5C; for detailed illustration of the fit quality, see the Supporting Information, Figures S6−S9.

Furthermore, four levels of Gaussian noise, std dev 0.14, 1.2, 2.5 and 4.6%, were added to the datasets constructed from the simulated LDM (Figure 6A and Figure S5B−E in the Supporting Information, cf. Figure S6−S9 in the Supporting Information for a better perception of the noise levels). The complexity of the datasets was extended by adding chirp and CA. The testing of the LDA algorithm was further extended to testing the performance of three different regularization terms, L, (ID, FD, and SD), as well as to testing the above-mentioned methods for determination of the regularization parameter, α, (L-curve, MPM, and GCV). The dimensions of the dataset matrices are the same as in the GLA testing. In the LDA module, the lifetimes are fixed (100 lifetimes spread on a decimal log scale) and thus they are not part of the fitting parameters. In this respect, only the IRF width, the chirp, and the CA are being fitted. The estimation accuracy of the 2334

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry Selection of the Regularization Factor α. The optimal selection of the regularization factor is a tedious task and needs to be performed carefully as α defines the balance between the regularization (smoothing) and the perturbation (noise) errors. In Figure 7 is shown an example for the determination of the

spectroscopy and shown the high efficiency and accuracy of the algorithms even for very noisy datasets. The focus of the work was placed on the implementation and the demonstration of the capabilities of the LDA technique for ultrafast spectroscopy. LDA provides a unique insight into the kinetic information contained in the experimental data without assuming a particular model a priori and thus can be used as a first analysis step especially in the case of complex kinetics. The methods is particularly advantageous in the case of nonlinear kinetics caused by, e.g., dynamic Stokes shift or cooling in ultrafast optical spectroscopy or by bimolecular reaction on longer timescales. Another valuable analytical property of the LDMs is that they actually can be used as kinetic footprints due to their high specificity for a given sample and experimental conditions. We have examined carefully the effect of the different regularization terms on the recovery of the LDMs under a number of scenarios and determined that the ID is most suitable for the discussed application. Furthermore, we have provided a selection procedure for the regularization factor to obtain the most reliable reconstruction of the LDM. The data analysis approaches and algorithms presented here can be readily implemented and applied not only for ultrafast spectroscopy but also for basically any time-resoled experimental technique.

Figure 7. Selection of the optimal regularization factor, α, based on the L-curve, MPM, and the GCV methods. The example shown here is for the synthetic dataset in Figure 6A (i.e., 1.2% noise).

optimal α for the 1.2% noise dataset using the three methods implemented here, L-curve, MPM, and GCV. In LDAOPTIMUS, such a figure is always generated to aid the user with the selection of α. Careful examination of the results from the three methods shows that for the type of data analyzed here the MPM typically suggests an α value that favors an oversmoothed solution. The values suggested by the L-curve and the GCV methods are relatively close although their difference changes depending on the regularization term that is used and the S/N level (Figure 8).

■

ASSOCIATED CONTENT

S Supporting Information *

Additional information as noted in text. This material is available free of charge via the Internet at http://pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected] Fax: +49 (0) 69/798 29709. Present Address †

H.H.: Frankfurt Institute for Advanced Studies, Frankfurt, Germany. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS The authors thank Dr. Markus Brown, Peter Trojanowski, and Peter Eberhardt (Goethe University) for the insightful discussions and the contribution to program development. C.S. and J.W. acknowledge DFG (Grant WA 1850/4-1).

Figure 8. Dependence of the optimal α, estimated by the different methods (L-curve, MPM, and the GCV), on the regularization term, L, and the S/N level; 1, ID; 2, FD; 3, SD.

■

REFERENCES

(1) Zewail, A. H. Angew. Chem. 2000, 39, 2586−2631. (2) Sundström, V. Annu. Rev. Phys. Chem. 2008, 59, 53−77. (3) Tamai, N.; Miyasaka, H. Chem. Rev. 2000, 100, 1875−1890. (4) Zinth, W.; Wachtveitl, J. ChemPhysChem 2005, 6, 871−880. (5) Kennis, J. T. M.; Groot, M.-L. Curr. Opin. Struct. Biol. 2007, 17, 623−630. (6) Middleton, C. T.; de La Harpe, K.; Su, C.; Law, Y. K.; CrespoHernández, C. E.; Kohler, B. Annu. Rev. Phys. Chem. 2009, 60, 217−239. (7) Wand, A.; Gdor, I.; Zhu, J.; Sheves, M.; Ruhman, S. Annu. Rev. Phys. Chem. 2013, 64, 437−458. (8) Bamann, C.; Bamberg, E.; Wachtveitl, J.; Glaubitz, C. Biochim. Biophys. Acta, Bioenerg. 2014, 1837, 614−625. (9) Klimov, V. I. Annu. Rev. Phys. Chem. 2007, 58, 635−673. (10) Snellenburg, J. J.; Laptenok, S. P.; Seger, R.; Mullen, K. M.; van Stokkum, I. H. M. J. Stat. Softw. 2012, 49, 1−22. (11) van Stokkum, I. H. M.; Larsen, D. S.; van Grondelle, R. Biochim. Biophys. Acta, Bioenerg. 2004, 1657, 82−104.

As previously observed,23,37 the GCV method gives a slightly lower value for the regularization factor than the L-curve, thus favoring a less smooth solution when ID is used; the behavior changes for FD and SD (Figure 8). Nevertheless, the combination of the two methods defines a range where the optimal α is located. Within this range, the solutions yield very similar LDMs, which do not contain strongly expressed over- or under-smoothing effects. The final LDM can be selected from within this range without any significant loss of information. The results in Figure 5 are obtained using α from within the range defined by the GCV and the L-curve.

■

CONCLUSIONS We have presented here the implementation strategies for a number of data analysis techniques for ultrafast optical 2335

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336

Article

Analytical Chemistry (12) van Wilderen, L. J. G. W.; Lincoln, C. N.; van Thor, J. J. PLoS One 2011, 6, e17373. (13) Jaumot, J.; Gargallo, R.; de Juan, A.; Tauler, R. Chemom. Intell. Lab. Syst. 2005, 76, 101−110. (14) Ruckebusch, C.; Sliwa, M.; Pernot, P.; de Juan, A.; Tauler, R. J. Photochem. Photobiol., C 2012, 13, 1−27. (15) Smilde, A.; Bro, R.; Geladi, P. Multi-way Analysis with Applications in the Chemical Sciences; John Wiley and Sons: Hoboken, NJ, 2004. (16) Holzwarth, A. R. In Biophysical Techniques in Photosynthesis. Advances in Photosynthesis Research; Amesz, J., Hoff, A. J., Eds.; Kluwer Academic Publishers: Dordrecht, The Netherlands, 1996; pp 75−92. (17) Arnaut, L.; Formosinho, S.; Burrows, H. Chemical Kinetics: From Molecular Structure to Chemical Reactivity; Elsevier: Amsterdam, The Netherlands, 2007. (18) Chen, C.-T. Linear System Theory and Design, 3rd ed.; Oxford University Press: New York, 1999. (19) Plaza del Pino, I. M.; Parody-Morreale, A.; Sanchez-Ruiz, J. M. Anal. Biochem. 1997, 244, 239−255. (20) Steinbach, P. J.; Ionescu, R.; Matthews, C. R. Biophys. J. 2002, 82, 2244−2255. (21) Voelz, V. A.; Pande, V. S. Proteins: Struct., Funct., Bioinf. 2012, 80, 342−351. (22) Mulligan, V. K.; Hadley, K. C.; Chakrabartty, A. Anal. Biochem. 2012, 421, 181−190. (23) Lorenz-Fonfria, V. A.; Kandori, H. Appl. Spectrosc. 2007, 61, 428− 443. (24) James, D. R.; Ware, W. R. Chem. Phys. Lett. 1986, 126, 7−11. (25) Livesey, A. K.; Brochon, J. C. Biophys. J. 1987, 52, 693−706. (26) Landl, G.; Langthaler, T.; Englt, H. W.; Kauffmann, H. F. J. Comput. Phys. 1991, 95, 1−28. (27) Croce, R.; Müller, M. G.; Bassi, R.; Holzwarth, A. R. Biophys. J. 2001, 80, 901−915. (28) Müller, M. G.; Slavov, C.; Luthra, R.; Redding, K. E.; Holzwarth, A. R. Proc. Natl. Acad. Sci. U.S.A. 2010, 107, 4123−4128. (29) James, D. R.; Ware, W. R. Chem. Phys. Lett. 1985, 120, 455−459. (30) Trojanowski, P.; Plotner, J.; Grunewald, C.; Graupner, F. F.; Slavov, C.; Reuss, A. J.; Braun, M.; Engels, J. W.; Wachtveitl, J. Phys. Chem. Chem. Phys. 2014, 16, 13875−13888. (31) Kovalenko, S. A.; Dobryakov, A. L.; Ruthmann, J.; Ernsting, N. P. Phys. Rev. A 1999, 59, 2369−2384. (32) http://www.theochem.uni-frankfurt.de/~chslavov/optimusfit/ and http://www.optimusfit.org/. (33) Dobryakov, A. L.; Kovalenko, S. A.; Weigel, A.; Perez-Lustres, J. L.; Lange, J.; Müller, A.; Ernsting, N. P. Rev. Sci. Instrum. 2010, 81, 113106. (34) Moore, E. H. Bull. Am. Math. Soc. 1920, 26, 394. (35) Penrose, R. Math. Proc. Cambridge Philos. Soc. 1955, 51, 406−413. (36) Beechem, J. M.; Gratton, E.; Ameloot, M.; Knutson, J. R.; Brand, L. In Topics in Fluorescence Spectroscopy; Lakowicz, J. R., Ed.; Plenum Press: New York, 1991; pp 241−305. (37) Hansen, P. C. Rank-Deficient and Discrete III-Posed Problems: Numerical Aspects of Linear Inversion; SIAM: Philadelphia, PA, 1998. (38) Tikhonov, A. N. Dokl. Akad. Nauk SSSR 1963, 151, 501−504. (39) Miller, K. SIAM J. Math. Anal. 1970, 1, 52−74. (40) Hansen, P. BIT 1990, 30, 658−672. (41) Lian, J.; He, B. Brain Topogr. 2001, 13, 209−217. (42) Golub, G. H.; Heath, M.; Wahba, G. Technometrics 1979, 21, 215−223. (43) Wahba, G. Spline Models for Observational Data; Society for Industrial and Applied Mathematics: Philadelphia, PA, 1990. (44) Steinbach, P. J.; Chu, K.; Frauenfelder, H.; Johnson, J. B.; Lamb, D. C.; Nienhaus, G. U.; Sauke, T. B.; Young, R. D. Biophys. J. 1992, 61, 235−245. (45) Ablonczy, Z.; Lukács, A.; Papp, E. Biophys. Chem. 2003, 104, 249− 258. (46) Vecer, J.; Herman, P. J. Fluoresc. 2011, 21, 873−881. (47) Manzo, A. J.; Goushcha, A. O.; Berezetska, N. M.; Kharkyanen, V. N.; Scott, G. W. J. Phys. Chem. B 2011, 115, 8534−8544.

(48) Esposito, R.; Altucci, C.; Velotta, R. J. Fluoresc. 2013, 23, 203− 211.

2336

DOI: 10.1021/ac504348h Anal. Chem. 2015, 87, 2328−2336