Biophysical Journal 84:2112-2129 (2003)
© 2003 The Biophysical Society
Application of Singular Value Decomposition to the Analysis of Time-Resolved Macromolecular X-Ray Data
Marius Schmidt*,
,
Sudarshan Rajagopal,
Zhong Ren,
,
and
Keith Moffat,
* Physik-Department E17, Technische Universitaet Muenchen, 85747 Garching, Germany; Department of Biochemistry and Molecular Biology, University of Chicago, Chicago, Illinois 60637 USA; BioCARS, Argonne National Laboratory, Argonne, Illinois 60439 USA; and Renz Research, Des Plaines, Illinois 60018 USA
Correspondence: Address reprint requests to Marius Schmidt, E-mail: marius{at}hexa.e17.physik.tu-muenchen.de.
 |
ABSTRACT
|
|---|
Singular value decomposition (SVD) is a technique commonly used in the analysis of spectroscopic data that both acts as a noise filter and reduces the dimensionality of subsequent least-squares fits. To establish the applicability of SVD to crystallographic data, we applied SVD to calculated difference Fourier maps simulating those to be obtained in a time-resolved crystallographic study of photoactive yellow protein. The atomic structures of one dark state and three intermediates were used in qualitatively different kinetic mechanisms to generate time-dependent difference maps at specific time points. Random noise of varying levels in the difference structure factor amplitudes, different extents of reaction initiation, and different numbers of time points were all employed to simulate a range of realistic experimental conditions. Our results show that SVD allows for an unbiased differentiation between signal and noise; a small subset of singular values and vectors represents the signal well, reducing the random noise in the data. Due to this, phase information of the difference structure factors can be obtained. After identifying and fitting a kinetic mechanism, the time-independent structures of the intermediates could be recovered. This demonstrates that SVD will be a powerful tool in the analysis of experimental time-resolved crystallographic data.
|
INTRODUCTION
|
|---|
Characterization of reaction intermediates is critical to understanding the pathways and mechanism by which a protein performs its biological reaction. Direct structural information on these intermediates is difficult to obtain because they are often unstable and have short lifetimes. Thus, methods used to probe them must either be capable of fast time resolution or increase the lifetime of the species to be studied. Chemical or physical trapping methods have been used to increase the lifetime (Moffat and Henderson, 1995
; Schlichting and Chu, 2000), always with the caveat that the nature of the trapping might disturb the true structure of the intermediates. In time-resolved crystallography, in contrast, no perturbation of the structure of the intermediates is required (except for those perturbations associated with the crystalline state, which is also true for trapping methods). The relative simplicity of trapping methods is given up in favor of a method that is technically challenging (Ren et al., 1999).
Excellent time resolution as low as 100 ps at third-generation synchrotron x-ray sources is possible using polychromatic Laue crystallography (Bourgeois et al., 1996), which allows the visualization of extremely short-lived intermediates. A number of time-resolved Laue studies have been performed with time resolutions varying from nanoseconds to milliseconds (reviewed in Ren et al., 1999). The most detailed of these are nanosecond pump-probe time-resolved studies on the photolysis of the CO-myoglobin complex (Srajer et al., 1996, 2001) and on the photocycle of the blue-light photoreceptor known as photoactive yellow protein (PYP) (Perman et al., 1998; Ren et al., 2001). In both studies, time-dependent difference Fourier maps are generated in real space from measured structure factor amplitudes as the reaction proceeds. Interpretation of such difference Fourier maps is not trivial. It is hindered by a low signal-to-noise ratio arising from error in the difference structure factor amplitudes and in the phase of the parent structure, and from the difference Fourier approximation itself (Henderson and Moffat, 1971). Signal may be difficult to differentiate from noise by simple visual inspection of the map (Moffat, 2001; Srajer et al., 2001). Furthermore, a difference map that corresponds to a single time point will consist of an admixture of difference features arising from all time-independent, intermediate structures that are significantly populated at that time. "Deconvolution", or separation of this mixture into pure, time-independent intermediates, is essential to determine the chemical reaction mechanism and the structure of each intermediate (Moffat, 1989).
Both issues, the differentiation of signal from noise and the separation of intermediates, may be addressed by a mathematical procedure commonly used in the analysis of time-resolved data, singular value decomposition (SVD) (Golub and Reinsch, 1970). SVD takes data—e.g., a set of optical absorption spectra or electron density obtained under different conditions, such as time, pH, or voltage—and represents it by two sets of vectors, which are weighted by their corresponding singular values. In time-resolved spectroscopy, for example, the "left" set of singular vectors (lSVs) constitute a time-independent orthonormal basis set from which all time-dependent difference spectra in the data matrix are constructed. The "right" singular vectors (rSVs) describe the time-dependent variations of the corresponding lSVs. The singular values correspond to the degree to which their respective lSVs and rSVs contribute, in a least-squares sense, to the data matrix. Because the vectors that model the data matrix are weighted by singular values, the data matrix can be approximated by a subset of singular values and vectors that contains primarily signal, thus reducing the noise present in the data. This procedure acts as a mechanism-independent filter of noise that is objective (up to the point at which the particular subset of significant singular values and vectors is chosen). The reduced representation of the data facilitates the interpretation of the rSVs with a chemical kinetic mechanism by means of a least-squares fit. This then allows the condition-independent (here, time-independent) intermediates to be obtained. The only requirement for the application of SVD is that the observable varies linearly with the concentration of the intermediates, which is the case with difference spectra and electron density.
SVD has been successfully used in a number of areas such as the analysis of spectroscopic data (Henry, 1997; Henry and Hofrichter, 1992, and references therein), of molecular dynamics simulations (Romo et al., 1995; Doruker et al., 2000) and the temporal variation of genome-wide expression (Alter et al., 2000). However, time-resolved crystallographic data differs substantially from, for example, time-resolved spectroscopic data. The signal-to-noise is much lower at a typical grid point in a difference electron density map than at a typical wavelength in a difference absorption spectrum; the signal tends to be concentrated at a subset of grid points, usually near the active site or chromophore of the protein, rather than distributing over a wide wavelength range, with the consequence that the majority of grid points exhibit only noise; and the difference electron density calculated from thousands of reflections is distributed over 105–106 grid points, rather than at 102 wavelengths. Further, time-resolved crystallographic data exhibit systematic as well as random errors, arising, for example, from the difference Fourier approximation (Henderson and Moffat, 1971) and from variation in the extent of reaction initiation from time point to time point.
It was therefore not obvious that SVD would be immediately applicable to time-resolved crystallographic data, nor how it should be implemented, and what the limitations on its successful use might be. We evaluate these issues here by preparing mock data that simulate a time-resolved crystallographic experiment on PYP. These data were generated for several chemical kinetic mechanisms of varying complexity in which the level of noise on the structure amplitudes, the number of measured time points, and the extent of reaction initiation were varied to simulate a realistic time-resolved experiment as closely as possible. Analysis of these mock data by SVD reveals that the method is powerful and can be successfully applied, within certain limits. These limits in turn suggest which types of noise are most important and how experiments should be conducted, if SVD is to be subsequently applied. The noise-reduced data can be further analyzed to identify the chemical kinetic mechanism and determine the time-independent structures of the intermediates in that mechanism.
|
A GENERAL GUIDE TO THE SVD ANALYSIS
|
|---|
The following paragraphs explore the applicability of SVD to time-resolved x-ray data by using mock data. In Section 1, Application of of SVD to time-dependent difference electron density, we present the principles. This paragraph mainly describes how the time-dependent difference electron density maps can be related to the columns of a single large matrix, which can be decomposed by SVD. In Section 2, Generation of the mock data, the generation of realistic mock data is outlined in detail. Special attention is given to the incorporation of noise of different origins. Section 3, Application of SVD to mock data, details the outcome of the application of SVD to the mock data. It shows the influence of noise on the magnitude and shape of the singular values and vectors, respectively. In Section 4, SVD as a noise filter, we demonstrate how SVD can contribute to the phasing of, and to a substantial reduction of noise in, the crystallographic difference maps. We call this method SVD flattening. Section 5, Extraction of mechanism, is the key section of this paper. Time-independent difference electron density maps are extracted by modeling or fitting the rSVs by several candidate, chemical kinetic mechanisms. Here, we pinpoint the limits of the SVD-driven analysis arising from different sources of noise. We demonstrate that an SVD analysis of experimental time-resolved crystallographic data is indeed feasible. In a last step, we elaborate a method to discriminate between compatible and incompatible candidate, chemical kinetic mechanisms. This method depends on representation of the time-independent difference maps on an absolute scale. This scale is not available in the maps extracted from the singular vectors. However, absolute scale may be restored after intermediate structures are determined. Therefore, we call this scheme "posterior analysis".
Fig. 1 presents a road map to analyze crystallographic difference electron density maps by means of SVD. It is intended as a reminder to a reader who is already familiar with the general concepts of this paper. The upper panel (for details read Sections 1–4) describes the largely objective mathematical decomposition of the data into singular vectors and values followed by noise reduction by means of SVD flattening. A successful fit of the right singular vectors by a sum of exponentials determines the minimum number of intermediates (Section 5). The middle panel describes the fit of the rSVs by chemical kinetic mechanisms and the assessment of the extracted time-independent maps by the criteria specified in Section 5. The scheme of "posterior analysis" in the lower panel discriminates further among the remaining mechanisms until the best mechanism(s) compatible with the data is (are) found. We will refer to this figure throughout the paper as different steps in the analysis are reached.
View larger version (32K):
[in this window]
[in a new window]
|
FIGURE 1 Road map for analysis of time-dependent crystallographic difference maps by SVD. (Upper panel) The data matrix A is decomposed by SVD and reconstituted with the most significant rSVs and lSVs. The resultant matrix A' can be SVD-flattened and fit by a sum of exponentials. If outliers are apparent, the original data matrix must be corrected; if not, determine the number of relaxation times. (Middle panel) List all chemical kinetic mechanism that are consistent with the number of relaxation times, fit a candidate mechanism, extract time-independent difference electron density maps, and assess the quality of such density using the three criteria mentioned in the text. If any of the maps fails a criterion, select another mechanism. (Bottom panel) If all criteria are passed, refine atomic structures of the intermediates. Calculate time-dependent difference maps using concentrations from the mechanism and the refined structures of the intermediates and the ground state. Compare these maps with the experimental difference maps by the mean-square deviation in the difference electron density. Choose the mechanism that exhibits the lowest deviation.
|
|
|
APPLICATION OF SVD TO TIME-DEPENDENT DIFFERENCE ELECTRON DENSITY: PRINCIPLES
|
|---|
Most time-resolved experiments are of the pump-probe type (Moffat, 2002). After reaction initiation by a brief laser pulse, the structure is probed at a time delay t with an intense polychromatic x-ray beam. The results of such experiments are N data sets of time-dependent x-ray structure amplitudes. By combining these with the structure factors of the dark/unperturbed state, a time-dependent difference electron density map 
(ti) is calculated for each time point ti (I = 1..N) at M equidistant grid points in the asymmetric unit. By assigning grid point m (m = 1..M) of map ti to the element a (a = 1..M) of vector ai, the data matrix A is generated. The assignment of grid point m on to element a must be consistent throughout the entire set of maps and the maps must be sorted in temporal order. Each column vector ai of data matrix A contains an entire difference electron density map for time point ti. By tracing the content of equivalent grid points in consecutive maps along the time axis, a time course of difference electron density is generated. The time courses are represented in the row vectors of the data matrix A. The SVD procedure then decomposes data matrix A according to Eq. 1 into an 
matrix U, each of whose N columns are called lSVs; the 
diagonal matrix S whose diagonal elements are the singular values; and the transpose of the square matrix V, VT, each of whose rows is an rSV:
 | (1) |
The data matrix A can be best approximated in a least-squares sense by matrix A' obtained by using the S < N most significant columns of U (denoted U'), rows of VT (denoted V'T), and elements of S (denoted S') (Eq. 2) for the reconstruction:
 | (2) |
In this way SVD acts as an effective noise filter by separating the signal contained in the first few lSVs used for reconstruction from the noise in the remaining insignificant singular vectors, omitted from the reconstruction. The number S of significant singular values/vectors can be judged by examination of the magnitude of the singular values and the magnitude of the autocorrelation of each row of VT. For data that contain systematic noise across the time axis, due to, for example, fluctuating laser power and extent of reaction initiation, the effectiveness of this noise filter can be enhanced by improving the autocorrelation, a measure of the function's smoothness, through the "rotation" algorithm described in detail by Henry and Hofrichter (1992), who also provide a comprehensive review of the SVD procedure.
|
GENERATION OF THE MOCK DATA
|
|---|
PYP is a small (14 kDa), water soluble blue-light receptor of the photoautotrophic purple bacterium Ectothiorhodospira halophila. It undergoes a fully reversible photocycle characterized by several spectroscopically distinct intermediates (Hoff et al., 1994; Ujj et al., 1998; Ren et al., 2001). The ultimate goal of a time-resolved crystallographic experiment on PYP is to identify the mechanism of this photocycle and to determine the structure of each intermediate. The number of spectroscopically distinct intermediates may not necessarily match the number of structurally distinct intermediates (Ng et al., 1995). Nevertheless, for the generation of realistic mock data, we use a dark state, I0, and three intermediates, I1, I2, and I3 (see Fig. 2). The photocycle is completed by the recovery of the dark state. For the atomic structure of the dark state, I0, we used the Protein Data Bank (Berman et al., 2000) entry 2phy (Borgstahl et al., 1995). The first intermediate, I1, was constructed from Protein Data Bank entry 3pyp (Genick et al., 1998) by depleting this structure of all hydrogens and superimposing it on I0. The entry 2pyr (Perman et al., 1998) and a modified PYP structure from entry 2pyp (Genick et al., 1997) were used as the second and third intermediate, I2 and I3, respectively. After removing water, all structures were transferred to the room temperature unit cell of wild-type PYP, which was used throughout the simulation (see Tables 1 and 2). From the dark state structure and the intermediate structures, we calculated structure factors FI0(h) and FI1(h)..FI3(h) to a resolution of 1.9 Å.
We selected two mechanisms to generate the mock data from a generic chemical mechanism for interconversion of three intermediates in a branched reaction (Fig. 3 A). These are the irreversible sequential mechanism S (Fig. 3 B) and the dead-end mechanism DE with I3 in a side path (Fig. 3 C). For each mechanism, we chose rate coefficients to create a simple and a more complicated case, for a total of four kinetic mechanisms. Tables 3 and 4 (below) list the rate coefficients used to integrate the coupled differential equations to determine the time-dependent fractional concentrations for each of the four mechanisms, as shown in Fig. 4. In sequential mechanism S1 (Fig. 4 A), the time at which the peak concentration of each intermediate is reached is very well separated, whereas in S2 (Fig. 4 B), the rate coefficient for the decay of I2 is chosen to be larger than that of I1, such that I2 is completely buried within I1 and attains only a low maximum concentration. If noise is present, the recovery of this intermediate from the data will be a particularly challenging test for the analysis. The mechanism DE employs two irreversible steps with rate coefficients k+1 and k+4 and a reversible reaction with the rate coefficients k+3 and k-3. This leads to the formation of two transients for I2. In mechanism DE1 (Fig. 4 C), the second transient for I2 starting at 
4 ms is barely visible, whereas in mechanism DE2 (Fig. 4 D), the second transient decays in parallel with I3. The relaxation time for the decay of the first transient of I2 is close to 1/k-3, whereas the second transient approximately parallels the decay of the equilibrium between I2 and I3, whose relaxation time is close to 1/k+4.
View larger version (13K):
[in this window]
[in a new window]
|
FIGURE 3 (A) General chemical kinetic mechanism for a branched reaction with three intermediate states I1, I2, and I3. The dark state I0 is excited by light. The first intermediate I1 relaxes back to the dark state via two other intermediates, I2 and I3. The direct path from I1 to I0 is not considered. (B) Irreversible sequential kinetic mechanism S; (C) Dead-end kinetic mechanism DE. An equilibrium between I2 and I3 is allowed as a side path or dead end; I2 relaxes back to the dark state.
|
|
View this table:
[in this window]
[in a new window]
|
TABLE 3 Extraction of time-independent difference maps and reaction coefficients from time-dependent difference maps generated from mechanism S1 under a variety of experimental conditions
|
|
View this table:
[in this window]
[in a new window]
|
TABLE 4 Extraction of time-independent difference maps and reaction coefficients from time-dependent difference maps generated from mechanisms S2, DE1 and De2 under a variety of experimental conditions
|
|
View larger version (24K):
[in this window]
[in a new window]
|
FIGURE 4 Time-dependent concentrations of intermediates calculated from the four different chemical kinetic mechanisms used in the simulations. Reaction coefficients are those given in Tables 3 and 4. Solid line: I1; dotted line: I2; dashed line: I3. (A) Irreversible sequential mechanism S1 in which peak concentrations of intermediates are well separated in time. (B) Irreversible sequential mechanism S2 in which concentrations of I2 are low and remain below those of I1 at almost all times. (C) Mechanism DE1 in which the second intermediate is in a side path and the second transient of I2 is negligible. (D) Mechanism DE2 in which there is a pronounced second transient of I2, and the concentrations of I3 and I2 are comparable at longer times.
|
|
For each of the four mechanisms, time-dependent structure factors F(h,t) were calculated by the vector addition of time-independent structure factors of the intermediates FIj(h) and the dark-state structure factor weighted by their individual time-dependent concentrations cj(t) according to Eq. 3:
 | (3) |
Note that Eq. 3 is general for any kinetic mechanism. J is the total number of intermediates and equals 3 in this instance. Due to stoichiometric constraints, the concentration of the dark state can be calculated from the concentrations of the intermediates. In typical simulations, 20% reaction initiation with cI1(0) = 0.2 was employed, which is a realistic experimental value for PYP entering the photocycle in the crystal. However, this value can vary from crystal to crystal, and hence in some simulations the initial concentration of I1 was randomly selected from time point to time point to be between 14% and 26% or between 5% and 17%. Values of t in Eq. 3 were chosen equidistant in logarithmic time to cover the entire time course.
Noise on the structure amplitudes |F(h,t)| was based on experimental standard deviations (
-values) measured in one particular reference Laue data set collected from a PYP crystal in the dark. To preserve the correlation between the magnitude of the structure amplitude and the magnitude of the -value, we defined 10 equally spaced bins between the minimum and maximum value of the structure amplitudes found in this data set. We then assigned the -values to these 10 bins. For each mock |F|, one of the bins was chosen based on the magnitude of the structure amplitude and a -value was picked randomly from the set of -values present in this bin. A multiple of the experimental -value was used as the width of a Gaussian distribution of error values. An error value was picked randomly from this distribution using the Box-Muller algorithm (Box and Muller, 1958) and added to |F(h,t)|. If the experimental data set did not contain an entry for a particular h, this reflection was discarded throughout the mock data sets. If the picked error value was negative, the structure amplitude could also become negative. In this case the amplitude was reset to 0.5 of the error value. This procedure creates "erroneous" or "noisy" time-dependent data sets of amplitudes, |F#(h,t)|, with standard deviation (|F#(h,t)|), with completeness identical to the experimental data sets (Tables 1 and 2). Dark-state structure amplitudes |FI0#(h)| with standard deviation (|FI0#(h)|) were derived from |FI0(h)| in an identical manner. In all cases, resolution was limited from 15.0 to 1.9 Å. Time-dependent difference structure amplitudes were then calculated from a "noisy" time-dependent data set and the "noisy" dark-state data set. Excitation by an intense laser pulse may generate strain in the crystal, from which streaky reflections result; and the integrated intensities of reflections having a lower peak photon count are more difficult to measure accurately. Therefore, we assigned higher noise values to the time-dependent amplitudes and lower values to the dark structure amplitudes. Each data set of difference structure amplitudes is labeled according to the amount of noise present in the time-dependent amplitudes and the dark-state amplitudes, respectively. For example, a 2 s/1 s data set corresponds to error values from Gaussians with width 2 for |F#(h,t)| and 1 for |FI0#(h)|. The crystallographic statistics for the experimental reference Laue data sets and the simulated data sets with various levels of noise are shown in Tables 1 and 2.
Each time-dependent data set and the dark data set were scaled together in 20 resolution bins using XMerge (McRee, 1999). Weighted difference amplitudes |Fw(h,t)| were calculated from the "noisy" differences |F#(h,t)| according to Eq (4):
 | (4) |
in which the weighting factor w(h,t) was calculated according to Ren et al. (2001):
 | (5) |
F(h,t) is the standard deviation of the difference amplitude and was calculated as:
 | (6) |
For each time point, time-dependent difference electron density maps (t) were calculated using |Fw(h,t)| and phases calculated from the dark-state PYP atomic structure I0 on a 0.75 Å x 0.75 Å x 0.65 Å grid to a resolution of 1.9 Å. Pure, time-independent, noise- and error-free intermediate difference maps, Rj (j = 1..3), were also calculated as a reference (see also Fig. 2) by using vector subtraction of structure factors of the dark structure FI0(h) from those of the three intermediates FIj(h).
|
APPLICATION OF SVD TO MOCK DATA
|
|---|
For each time-dependent difference map, one entire PYP molecule plus a margin of 1 Å was masked out. The difference electron density at the grid points within this mask at all time points constitutes the data matrix A and was subjected to SVD (Fig. 1, Objective SVD). Typically, the mask contains 86,000 grid points, but larger masks containing 200,000 grid points could also be successfully used. To execute the necessary steps for an SVD of our time-dependent difference maps, we developed a program, Singular Value Decomposition for Time-Resolved Crystallography, written in Fortran 77. A data matrix having 21 column vectors with 86,000 grid points each can be generated from difference maps and decomposed in 0.5 min. As an example, Fig. 5 shows the first six lSVs obtained from SVD of the 2 s/1 s data of the mechanism S1 with three intermediate states. Signal should be present in the three most significant lSVs only. In panels lSV1–lSV3, the signal is very well defined and clusters where the intermediate atomic structures differ from the dark state. However, due to the noise, some signal has spread to the fourth lSV: there is some difference density on the tail of the chromophore and at other locations, such as on the phenolate oxygen of the chromophore ring. lSV5 and lSV6 contain only noise; the electron density features are scattered randomly in space and do not cluster. As an additional test, difference map SV5-21 was reconstructed from the 17 least significant singular vectors. Here, too, difference electron density features are scattered and do not occupy chemically sensible positions. In contrast, if the difference electron density maps were reconstituted from the first four singular vectors, strong spatially contiguous signal can be observed. SV1-4 in Fig. 5 shows the resultant difference map for time point 15 of the time course. For comparison, the corresponding mock difference map TP15 is also shown. The agreement between SV1-4 and TP15 is excellent.
View larger version (89K):
[in this window]
[in a new window]
|
FIGURE 5 The lSVs derived from SVD analysis of mock data generated with kinetic mechanisms S1 and DE2 on the 2 s/1 s noise level. The sign of the lSVs is arbitrary; that of the reconstruction is correct. (A, B, and C) First three most significant lSVs, contour level red/white -3 /-4 , blue/cyan 3 /4 ; (D, E, and F) at a lower contour level, red/white -2 /-3 ; blue/cyan 2 /3 . (G, H, and I) Singular vectors lSV4–lSV6, contour level red/white -2 /-3 ; blue/cyan 2 /3 . Arrows in lSV4 indicate signal. (J) Difference electron density SV5-21 at time point 15 is reconstructed from 17 least significant singular values and singular vectors lSV5–lSV21; contour level red/white -2 /-3 ; blue/cyan 2 /3 . (K) Difference electron density SV1-4 at time point 15 is reconstructed from the four most significant singular values and singular vectors lSV1–lSV4; contour level: red/white -3 /-4 ; blue/cyan 3 /4 . (L) Original difference electron density TP15 at time point 15; contour level: red/white -3 /-4 ; blue/cyan 3 /4 .
|
|
Fig. 6 shows how the singular values and the rSVs behave as noise is progressively added. On the left side, the magnitude of the singular values (circles) is shown together with the autocorrelation of rSVs (squares). On the right side the first few rSVs are plotted. Panels A and B show results from noise-free data. Only three significant singular values can be detected unambiguously, as expected. The same conclusion holds for the amplitudes of the rSV (Fig. 6 B). The high autocorrelation of rSV4 to rSV7 shows that they still contain smoothly varying amplitudes. However, the amplitudes of both the rSVs and lSVs are negligible compared to the true signal. At the 1 s/0.5 s noise level (Fig. 6, C and D), the system is still well-behaved. No significant amplitudes except those of the first three rSVs can be observed. When the noise is further increased to the 2 s/1 s level (Fig. 6, E and F), some signal can be observed in an additional, fourth lSV, and a fourth rSV begins to rise. When the noise is further increased, there might be a significant fifth singular value (Fig. 6 G). However, inspection of the lSVs after application of the rotation method confirms that only four singular values and vectors are needed in the reconstitution. At high noise, the autocorrelation becomes a sensible additional criterion to select the number of significant basis vectors. At the highest noise level of 10 s/5 s (Fig. 6, I and J), the fourth lSV contributes even more to the signal and a fifth significant lSV emerges. However, reconstruction of maps with more than five basis sets is still unnecessary; the map reconstructed with all the least significant vector sets, SV6-21, does not show any significant signal (data not shown).
View larger version (44K):
[in this window]
[in a new window]
|
FIGURE 6 (Left column) Dependence of the magnitude of the singular values S (•) and the autocorrelation of the corresponding rSV ( ) on the ordinal number of S, for mock data generated with kinetic mechanism S1.  , magnitude of the singular value after rotation;  , autocorrelation of the corresponding rSV after rotation. (Right column) Dependence of the magnitude of the corresponding rSVs on time. •, first rSV; , second rSV;  , third rSV;  , fourth rSV; x, fifth rSV; *, sixth rSV;  , seventh rSV. (A and B) 0 s/0 s, no noise present. (C and D) Noise level 1 s/0.5 s; (E and F) noise level 2 s/1 s; (G and H) noise level 5 s/3 s; (I and J) noise level 10 s/5 s.
|
|
At noise levels at and higher than 5 s/3 s, the identification of significant singular values is hindered by an offset. The ratio of the magnitude of the significant singular values to this offset becomes increasingly unfavorable. Fig. 6 clearly shows this offset at the 5 s/3 s and 10 s/5 s noise levels (Fig. 6, G and I): the rSVs begin to deviate strongly from zero at the end of the reaction. Less significant rSVs become larger in amplitude and vary rapidly in time; however, the significant rSVs still vary smoothly. As shown by Henry and Hofrichter (1992) the 5 s/3 s and higher noise levels may be too high for a successful application of SVD, because the offset becomes comparable to the amplitude of the next significant rSVs. Nevertheless, even at the 10 s/5 s noise level, a reconstruction of the time-dependent difference maps with a sensible selection of a set of singular values and vectors is feasible because it is possible to identify the signal easily by inspecting the lSVs.
The number of significant singular values and vectors can be selected by examining the amplitude of the singular values and the autocorrelation of the rSVs and by assessing the quality of the lSVs. The first two criteria must be used if no other evidence exists to identify signal in the lSVs. The third criterion is usually not available for spectroscopic data. For example, the shape of an lSV difference absorption spectrum might obey general constraints but it also may vary within large boundaries. Therefore, a difference absorption spectrum observed in the lSVs may or may not contribute significantly to the signal. In marked contrast, the quality of the lSVs can readily be assessed in crystallographic data. Signal tends to cluster on or near to atoms of the chromophore or of the active site. Thus, the decision on which singular values and vectors are significant is facilitated.
|
SVD AS A NOISE FILTER: PHASE IMPROVEMENT
|
|---|
The two difference maps shown in Fig. 7 demonstrate one of the most important properties of the SVD, the ability to separate signal from noise. These difference maps were calculated from kinetic mechanism S1 on the highest, 10 s/5 s level of noise. Fig. 7 A shows the original mock difference map; the map in Fig. 7 B was reconstituted from the first five lSVs derived from SVD. The dark-state atomic structure and the structure of I3 are shown as a guide to the eye. In the lower left corner, one can easily identify the authentic difference electron density on and close to the chromophore. In Fig. 7 A, numerous additional electron density features arising largely from noise are distributed randomly in space. In contrast, Fig. 7 B is nearly free of these features but preserves all of the signal on the chromophore. SVD is able both to reduce the noise in difference maps and to fully preserve the signal.
There are three sources of noise in a difference Fourier map: noise in the difference structure amplitudes, errors in the native (dark) phase, and the phase error introduced by the difference approximation itself (Henderson and Moffat, 1971). The first two are random but the third is systematic. Fig. 8 A illustrates the effect of the difference Fourier approximation on attempts to measure the true difference structure factor FT. FT is derived from the vector sum of the dark-state and intermediate-state structure factors weighted by their concentrations, and hence is time dependent. The difference Fourier approximation changes both the phase 
T and amplitude |FT| of the true difference structure factor to DA and |F|, respectively. To examine the effects of noise due to the difference Fourier approximation alone, we examined difference maps from noise-free data generated from kinetic mechanism S1. Difference structure factors were generated by the difference approximation. The mean absolute phase difference 
|T - DA|
was found to be 45° for acentric and around 43° for all reflections in all data sets, respectively. The values of T - DA are uniformly distributed between -90° and 90° (Fig. 9 A), which explains the above result. All noise-free difference maps can therefore be considered as conventional Fourier maps calculated with the relatively small phase error of 45° for acentric reflections and essentially correct phases for centric reflections. Fig. 9 B shows the effect on the phase difference T - DA if noise is added. Obviously, the box function in Fig. 9 A is convoluted with a Gaussian; the phase differences are distributed to larger values and the mean absolute phase error is increased to 56° (Fig. 9 B) and 66° (Fig. 9 C). When the signal-to-noise ratio is very low, for example in the final time points, the mean absolute phase difference approaches 90°, as expected for a pair of randomly varying phases on the unit circle.
View larger version (11K):
[in this window]
[in a new window]
|
FIGURE 8 (A) Argand diagram for a reflection h to illustrate the construction of the time-dependent structure factor F from the dark-state structure factor FI0 and the intermediate structure factors FI1, FI2, and FI3, each weighted by their corresponding time-dependent concentrations cIj, resulting in vectors 0, 1, 2, and 3, respectively; FT, true difference structure factor; F, difference structure factor after difference Fourier approximation; DA, phase from difference approximation; T, true phase of the difference structure factor. (B) Influence of noise: FI0 changes to FI0# and F to F#; the true difference structure factor FT becomes the noisy difference structure factor FN and, after the difference Fourier approximation, F#.
|
|
If noise introduced by the difference approximation were independent of time, SVD would be the right tool to minimize this noise and recover the magnitude and the true phase of the difference structure amplitude. To test this idea, we first analyzed difference maps calculated from 0 s/0 s data by SVD. SVD was applied to the time-dependent difference maps and the maps were reconstituted with the first three basis vectors. The mean phase difference |T - DA| decreased by 5° to |T - SVD|; phases SVD were derived from a Fourier back-transformation of the reconstituted maps. Subjecting the maps to SVD thus has recovered part of the phase information necessary to correct for the difference approximation. We were able to recover more phase information and lower the average phase difference further by applying a phase recombination scheme (Fig. 10) employed by Ren et al. (2001). Here, additional information is supplied by using the dark-state structure factor FI0 and the time-dependent structure amplitude |FL|. The value of FSVD, whose amplitude and phase are calculated from a Fourier back-transformation of the reconstituted maps, is added to FI0. The resultant vector (dotted line in Fig. 10) determines the phase FL of the time-dependent structure amplitude |FL|. Due to the phase error in FSVD, the triangle characterized by FL, FI0, and FI0 + FSVD does not close. We adjust FSVD so that the resultant F* points to FL with phase *. New time-dependent difference maps can be calculated from the difference structure factors F*. If time-independent noise features still persist, the phase error can be reduced further by iteration. Fig. 11, A and F, shows the result after a second cycle: the phase improvements remain rather small. We conclude that a partial correction of error introduced by the difference approximation is completed within one SVD iteration cycle. Additional cycles do not contribute to a significant further reduction of this error.
View larger version (41K):
[in this window]
[in a new window]
|
FIGURE 11 Phase improvement by SVD flattening. (A–E) Dependence of the mean absolute difference in phase angle between the true value T and the current value in cycle j, cycle j, on time, i.e., on the ordinal number of the data set. (A) 0 s/0 s, no noise present; (B) noise level 1 s/0.5 s; (C) noise level 2 s/1 s; (D) noise level 5 s/3 s; (E) noise level 10 s/5 s; (F) dependence of this mean absolute difference in phase angle on cycle number for time point/data set 15. •, starting value; , value after first cycle; , value after second cycle; , value after third cycle; x, value after fourth cycle.
|
|
We proceed further by adding realistic noise to the data. A time course of difference maps was generated from kinetic mechanism S1 with various levels of noise. From Fig. 8 B, one can see that FT will change to F# as noise is added. In Fig. 11, the mean absolute phase difference |T - DA| is shown as |T - cyclej|. If noise is added, all values of this phase difference lie significantly above the value of 43° observed in the noise-free data. Further, the values vary in a time-dependent fashion. For example, |T - cyclej=0| observed in the 10 s/5 s data set (Fig. 11 E) is 78° at the first time point, rises to 86° at the seventh, decreases to 74° at the 15th, and finally reaches 90° in the last time points. This time dependence can be understood in light of the underlying structural changes. When an intermediate characterized by a small structural difference from I0 is populated, the overall signal-to-noise ratio in the difference structure amplitudes decreases and a larger phase error results. In kinetic mechanism S, intermediate I1 is followed by I2, which decays to I3. The structure of I2 deviates from I0 only at and close to the chromophore, whereas the structures of I1 and I3 differ much more extensively from I0. Therefore, at time points 5–10, when I2 is most populated, the phase of the difference structure factors is more erroneous, but, in the last time points 19–21, the structural perturbation vanishes completely. Consequently, the mean phase error approaches 90° as expected from completely randomized phases.
However, SVD can contribute to the determination of the true difference structure factor FT from F# even in the presence of a large amount of noise. After SVD, phase recombination can be done similar to the procedure demonstrated for the 0 s/0 s maps in Fig. 10. In addition, the noise in the time-dependent structure amplitude |F#| can be subjected to appropriate weighting, as outlined by Ren et al. (2001). The result are improved difference structure amplitudes |Fcyclej*| and phases cyclei, which can be used to calculate new time-dependent difference maps, and the procedure iterated. The phase improvements are in the order of 10°–15° as demonstrated for time point 15 in Fig. 11 F. Two to three cycles are usually sufficient for acceptable convergence, depending on the noise level.
We call this method SVD flattening, by analogy with solvent flattening (Wang, 1985). Difference electron density features originating either from phase error or from random noise in the difference structure amplitudes are partially excluded into less significant singular vectors. In contrast, those features persisting over several time points are most likely signal and retained in the maps. In SVD-flattened maps, the difference density is better defined on a higher contour level, the connectivity of the difference electron density is improved, and noise features reduced by comparison with the original maps. The necessary information is gathered automatically from the time axis in a purely analytical fashion.
The final result of the SVD procedure, with or without SVD flattening, is a set of noise-reduced (or signal-to-noise enhanced) phased difference maps, which can be used for further analysis of the mechanism (Fig. 1, Prepare data matrix A'').
|
EXTRACTION OF MECHANISM FROM DIFFERENCE FOURIER MAPS: INTRODUCTION
|
|---|
The goals of any time-resolved spectroscopic or crystallographic experiment are the characterization of spectral or structural intermediates and the identification of their associated chemical reaction mechanism. Intermediates can be extracted from time-resolved data if their time-dependent concentrations can be described accurately, but in most instances it is not possible to extract them from complicated time traces without prior or additional knowledge (Lozier et al., 1992). In spectroscopic experiments, additional knowledge such as the nonnegativity and shape of absorption spectra, amplitude constraints, and temperature dependence may be sufficient for an unambiguous analysis of the data (Zimanyi and Lanyi, 1993; Nagle et al., 1995, Van Brederode et al., 1996). To test whether such an approach is feasible with time-resolved crystallographic data, we analyzed our simulated time-dependent difference electron density. The best time-dependent difference maps are represented by a noise-reduced or SVD-flattened data matrix A' or A'', which can be decomposed into matrices U', S', and V'T (Eq. 2). The columns of U' are the lSVs, the rows of V'T are the rSVs, and the singular values comprise the diagonal matrix S'. The N rSVs, vn, represent the temporal variation of the corresponding lSVs. Least-squares fitting of correctly weighted rSVs by some function is mathematically equivalent to fitting the entire data matrix using global analysis, with the advantage that the number of parameters necessary to describe the fit has decreased dramatically (Henry and Hofrichter, 1992). If a simple chemical kinetic mechanism holds (Moffat, 2001), the fit function must obey a sum of exponentials describing simple relaxation processes. The minimum number Q of intermediates in the mechanism then depends on the number of relaxation times observed by analysis of the rSVs. For example, three relaxation times are predicted for all kinetic mechanisms outlined in Fig. 3. The elements vn,i of the S significant rSVs, vn (n = 1..S), are fit by a sum of exponentials whose amplitudes Dn,q (q = 1..Q) plus an offset Dn,0, which differ for each rSV, and a unique set of relaxation rates, kq, simultaneously used to fit all significant vn. The fit is weighted by the square of the corresponding singular value sn according to Eq. 7:
 | (7) |
It is crucial to subsequent analysis of the time-resolved crystallographic data that the number of distinct relaxation times, common to all rSVs, can accurately be determined, as well as the amplitudes of each relaxation process for each rSV. We therefore explore how the magnitude of experimental errors influences the accuracy of such a determination.
Determination of relaxation times in the presence of noise
We identified the variation in the extent of reaction initiation (RI) as a major source of error in the time domain. Data are typically collected on different crystals for each time point. The laser pulse energy and the energy density (mJ/cm2) are prone to significant fluctuations, as data collection is lengthy and data often are combined from different experimental trials. To analyze the consequences of these realistic experimental variations, we allowed the number of time points in the data set to vary from five to 21 as the noise present in the structure amplitudes varied from the 1 s/0.5 s level to the 10 s/5 s level. In addition, the extent of RI for each time point was varied between 14% and 26% and in another case from 5% to 17%. After applying SVD to the mock data sets and preparing data matrix A' with the best-difference maps (Fig. 1, Prepare data matrix A'), the number of significant singular vectors was chosen based on the magnitude of the singular values, on the autocorrelation of their corresponding rSVs, and on visual inspection of the lSVs, as explained above. The relaxation times were initially determined from the time courses in the rSVs, and the sum of exponentials was fitted to the rSVs, as shown in Fig. 12, A–D, for data on the 1 s/0.5 s–10 s/5 s noise level. Even at the 10 s/5 s level, the three predicted relaxation times are clearly observable despite the large offset at the end of the reaction. The fit is satisfactory with a mean-square deviation (MSD) of 385 (see Tables 3 and 4). If only two relaxation times were used, an MSD of 2155 was obtained; four relaxation times yielded an MSD of 322, quite similar to the one obtained from three relaxation times. We conclude that the minimal number of relaxation times and intermediates can be accurately selected even at high experimental noise.
When the RI was varied, the data points in the rSVs scatter around the fitted curve. Fig. 12, E and F, illustrate the effect when both a large variation of the RI from 5% to 17% and two time points as outliers (time point 4 with 0.1% RI and time point 11 with 40% RI) are included in simulation on the moderate 2 s/1 s noise level. Both outliers are clearly observable in the rSV (arrows). However, a fit with three exponentials is now difficult to distinguish from one with two exponentials; the MSD of both is comparable (see Tables 3 and 4). Because these two data points are clear outliers, we can adjust the initial data matrix A to account for the systematic error. The difference electron density at time point ti is multiplied by a factor determined by the ratio of the fit value from the sum of exponentials and the magnitude of the element of the first rSV at this time point (Fig. 1, Correct for extent of reaction initiation). The result is data matrix A', whose rSV after a second cycle of SVD exhibit much smaller variations (Fig. 12 F). We conclude that variation in the extent of RI propagates through the SVD analysis into the rSVs. If the variation is sufficiently extreme, outliers can be identified, corrected in the original data matrix, and the SVD analysis repeated. That is, used judiciously, SVD can reduce certain sources of noise arising from systematic as well as from random error. If the random noise increases further as shown for the 5 s/3 s level, SVD flattening can be used in addition to increase the reliability of the relaxation times (Fig. 12, G and H; Tables 3 and 4).
Fitting a chemical kinetic mechanism to the rSV
The ability to accurately fit the rSVs by a sum of exponentials already limits the number of possible mechanisms that can account for the data. Only those reaction mechanisms that produce the observed number of relaxation times/rates can be used for further analysis (see Fleck, 1971). One general reaction mechanism for a reaction involving four states, I0, I1, I2, and I3, and three relaxation rates is shown in Fig. 3 A. All possible mechanisms for this case can be obtained from this general mechanism by setting some of the rate coefficients to zero. Two examples are given in Fig. 3, B and C. If a candidate reaction mechanism with J intermediates and rate coefficients kr (r = 1..R) is chosen, the time-dependent concentration of the jth intermediate, cIj(kr,t), can be calculated at time point ti by solving the coupled differential equations describing the mechanism, which is given in the most general form by:
 | (8) |
where the amplitudes Ap,j and exponents Bp,j are dependent on the set of rate coefficients kr (Matsen and Franklin, 1950). By varying the magnitude of the rate coefficients, the concentrations are fit to the elements of the most significant rSVs using the scale factors En,j for each concentration (Eq. 9):
 | (9) |
The time-independent difference electron density Ij of a candidate intermediate Ij, represented by the vector bIj, is then calculated according to Eq. 10 (see also Henry and Hofrichter, 1992):
 | (10) |
Each of the finite number of candidate mechanisms that exhibit J intermediates may be subjected to this process (Fig. 1, List all possible mechanisms). Each candidate mechanism then yields its own set of candidate intermediates Ij, represented by Ij. How shall these candidate intermediates be assessed? If all intermediates for a particular mechanism pass the assessment, the candidate mechanism and intermediates are consistent with the data; but if one or more intermediates fail the assessment, the candidate mechanism and the intermediates are rejected. This type of assessment lies at the heart of all investigations of mechanism based on kinetic data.
After fitting the rSV with candidate mechanisms (Fig. 13, A–F), time-independent difference electron density was extracted and assessed using three criteria.
View larger version (42K):
[in this window]
[in a new window]
|
FIGURE 13 Fit of the candidate chemical kinetic mechanisms S or DE to the three or four most significant rSVs with various levels of noise. •, first rSV; , second rSV; , third rSV; , fourth rSV. The solid lines are a fit to a sum of three exponentials. (A) Noise level 1 s/0.5 s, candidate mechanism S; (B) noise level 1 s/0.5 s, mechanism DE; (C) noise level 2 s/1 s, candidate mechanism S; (D) noise level 5 s/3 s, candidate mechanism S; (E) noise level 10 s/5 s, candidate mechanism S; (F) noise level 5 s/3 s, variable level of reaction initiation between 5 and 17%, corrected for outliers, after SVD flattening, candidate mechanism S.
|
|
As a first criterion, the extracted time-independent difference maps Ij were compared at all grid
TOP
ABSTRACT
INTRODUCTION
