10713
J. Phys. Chem. 1991,95, 10713-10718
+
constants for the reaction of F CH30H. MP2 level calculations show no classical barrier to hydroxyl-side abstraction and only a small barrier to methyl-side abstraction. This correlation between reaction exothermicity and classical barrier heights for R l versus R2 is opposite to the usual trend for hydrogen abstraction reactions. Our calculations yield methoxy branching fractions that reproduce both qualitatively and quantitatively the experimental results in paper 1 in this series. The improved quantitative agreement between theory and experiment when the barrier for R2 is removed suggests that the calculated barrier for R2 may be artifactual. Our analysis indicates that entropic rather than enthalpic considerations are responsible for the anomalously large methoxy yields observed experimentally. Although the calculated rate constants k, and k2 may entail appreciable errors in an absolute sense, we have argued that the calculated branching fractions should remain qualitatively reliable. In view of the uncertainty in the location of TS 1 and to ensure that we did not
obtain the right answers for the wrong reasons, a more detailed look at the potential surface using multireference wave functions and variational TST calculations is needed. The results of this investigation should provide a useful starting point for calculations a t higher levels of theory. Acknowledgment. We are indebted to Dr. Joseph L. Durant, Jr., for suggesting and inspiring this theoretical adjunct to his experimental studies. We thank Dr. Celeste M. Rohlfing, Dr. Carl F. Melius, and Dr. Richard S. Judson for providing valuable suggestions in the preparation of the manuscript. We must also acknowledge the important suggestionsmade by the reviewers that served to markedly enhance the clarity and accuracy of the presentation. The work was performed under the auspices of the Division of Chemical Sciences, Office of Basic Energy Sciences, US.Department of Energy. Registry No. F, 14762-94-8; CHIOH, 67-56-1; H, 1333-74-0.
Differential Method in Chemical Kinetics
c. Steel* Department of Chemistry, Brandeis University, Waltham, Massachusetts 02254-91 10
and K. Razi Naqvi Department of Physics, University of Trondheim, N-7055 Dragvoll, Norway (Received: June 10, 1991)
The following problem is addressed: given a digitally stored data stream and a reaction scheme, evaluate the rate constants without undertaking analytic or numerical integration of the relevant rate equations. Although the basic principles have long been known, and it is common knowledge that smoothing or differentiation according to the cannon of least squares is tantamount to a convolution, there appears to be no criterion for choosing W,the total number of the convoluting integers, i.e., the number of data points within a sampling window. We have found that the best course is to employ a whole range of W, the wider this range, the greater the opportunity for minimizing statistical as well as systematic errors. Using simulated data, the following schemes, the first of which is also amenable to analytic integration, have been successfully analyzed: (i) C k P,,2 C 4 P2; (ii) A C, followed by the previous two steps. The procedure is also applied to a problem of biophysical interest, viz., the determination, by means of a photoselection experiment, of the axial ratio of a randomly labeled prolate ellipsoid of revolution.
Introduction The extraction of rate constants from experimental data is a fundamental step in chemical kinetics. Most reaction schemes can be expressed in terms of coupled ordinary differential equations, but only a handful of these rate equations can be solved analytically and the integrated forms then utilized for the evaluation of rate constants. The use of the differential rate equation itself for the determination of rate constants and reaction orders was pioneered over 100 years ago by van’t Hoff,’ but the method has not found wide acceptance because of the difficulty of obtaining slopes from noisy or sparse data. It has been used to study the mixed first- and second-order decay of excited states and radicals,2-5 although-as we shall see below-other routes, including direct integration of the rate equation, are Now that large data sets can be gathered, stored, and rapidly analyzed, it is worthwhile to (a) ask if the differential method has general utility, (b) explore ways of optimizing the extraction
of rate parameters and their attendant uncertainties, and (c) devise a procedure for taking account of systematic errors inherent in the very process of digital differentiation and smoothing. So as to be able to make comparisons with a case which can be handled analytically, we shall first consider the above-mentioned mixed first- and second-order consumption of a species. Thereafter, we shall use the principles developed en route to analyze a scheme not amenable to integration. Finally, we present a biophysical example and conclude by suggesting that the method can be profitably brought to bear on other kinetic problems.
Background The reaction scheme and rate equation for the mixed decay are
c
-
2c(1) van’t Hoff, J. See, for example: Laidler, K. J. Chemical Kinetics, 2nd ed.; McGraw-Hill: New York, 1965; p 15. (2) Linschitz, H.; Steel, C.; Bell, J. A. J . Phys. Chem. 1962, 66, 2574. (3) Heppel, G. E. Photochem. Photobiol. 1965, 4, 7. (4) Grellmann, K. H.; Scholz, H.-G. Chem. Phys. Lett. 1979, 62, 64. (5) El-Sayed, F. E.; MacCallum, J. R.; Pomery, P. J.; Shepherd, T. M. J . Chem. SOC.,Faraday Trans. 2 1979, 75, 79. (6) Porter, G.; Wright, M. R. J . Chim. Phys. 1958, 55, 705. (7) Burkhart, R. D.; Boileau, S.;Bovin, S . J. Phys. Chem. 1987, 91, 2189.
0022-3654/91/2095-10713$02.50/0
kl
PI
ki P2
d@(t)/dt = &(t) = -k,@(t) - k2e2(t) where @(t), the concentration of species C at time expressed, with Bo e(O),as k,@0 exp(-k,t)
e(?)= k l + k26,[1 - exp(-klt)] 0 1991 American Chemical Society
(3) 1,
can be
(4)
10714 Thp Journal of Physical Chemistry, Vol. 95, No. 26, 1991
Steel and Naqvi
Various techniques have been applied to the extraction of kl and k2.2-7 One637is based on rearranging eq 4 to exp(k,r) = A
+ B/@(t)
[A/B = k2/kl]
and determining kl from the long-time portion of the decay curve where the contribution of k2 is mininum. However, as has been pointed this way of determining k , is precarious. Another m e t h ~ d ~is. ~to. rewrite ~ the rate equation as -d[ln @(t)]/dr = k ,
+ k2@(t)
0.
m=3 true slope
O.0. “0. I
-3
Originally,2 the analogue data were replotted and the slopes determined by graphical means, a tedious and error-prone exercise. Since the advent of digital data acquisition, the discretized form of the left-hand side, A[ln @(t)]/At, has been useds with a fixed window (vide infra).
0
8
0’0.03
W
0-0. “0-0. W
0
3
6
12
9
15
18
Sample number ( i )
Methodology Convolution: Smoothing and Differentiation by Least Squares (LS). Consider a data set in tabular form with each row listing concentrations and times: [@(i),r(i)], where i, the sample number, takes integral values from 0 to N . Our first task is to produce from these ( N 1) primary data points a secondary set of slopes and concentrations [ e 0 ) , e 0 ) ] , where the underlining serves to emphasize that the original data have been subjected to a convolution or smoothing operation. The mathematical basis for least-squares polynomial (LSP) differentiation and smoothing of equally spaced data is well established and amply dicussed.8-’1 By far the commonest application of LSP convolution is in replacing a noisy experimental curve such as a spectrum or a chromatogram-having, more often than not, several peaks and valleys-by a smoothed version; to this end the parameters that may be adjusted include W (the number of convoluting intege:s or the number of points in the smoothing window), the degree of the smoothing polynomial, and the number of repetitions of the process. The influence of these parameters on the extent of noise reduction and signal distortion was investigated empirically by Enke and Nieman,Io who reported the results of their study in a review paper dealing with the capabilities and the limitations of LSP smoothing. Some pitfalls of LSP smoothing, particularly the distortion resulting from the application of long LSP sequences with a low-degree polynomial to rapidly varying data, have been pointed out by Marchand and Marmet.” A kineticist does not wish to doctor a noisy data stream and end up with a smoother curve that is simply more pleasing to the eye; such curve fitting would, in the words of Enke and Nieman, be “of cosmetic value only”. The motive behind analyzing noisy rate data by means of the differential method is to procure rate constants suffering as little as possible from the distortion that afflicts LSP slopes and concentrations; this purpose can be met only by developing a new strategy. The spirit, if not the scope, of this paper is similar to that of the Enke-Nieman review: we wish to examine empirically the potential and drawbacks of LSP fitting for the purpose of extracting rate constants from kinetic data. The essentials of LSP smoothing are recapitulated here with the aid of Figure 1, in which W(=2m + I), the number of points in a window, is 7. In general, m can have integral values 1, 2, 3, .... The first step is to calculate e(i = m) and = m),the concentration and slope at the midpoint of the window. This is most easily done by defining a new index w, which runs from -m to +m, in terms of which the data points in the window are fitted to some smoothing polynomial
+
c(i
+ pw + yw2 + ... (5) Thus, at the midpoint of the window one has e ( w = 0 or i = m) @ ( w ) = CY
= a and @(w = 0 or i = m ) = @. (8) Kerawala, S. M. Indian J . Phys. 1941, 15, 241. (9) Savitzky, A.; Golay, M. J. E. Anal. Chem. 1964, 36, 1627. (IO) Enke, C. G.;Nieman, T. A. Anal. Chem. 1976, 48, 705A. ( 1 1 ) Marchand, P.; Marmet, L. Reu. Sci. Instrum. 1983, 54, 1034.
Figure 1. Concentration of species C a s a function of sample number, i. The index w is used for labeling within a window. Inset a shows how smoothing, in this case first-order, results in systematic errors in both the concentration and slope associated with the midpoint of the window. W is the number of data points in a window.
The subsequent steps in smoothing usually consist in advancing the window one point at a time and determining e(m),e(m l), e ( m + 2), and so o ~ . ~ - I ’ We shall see later that, for kinetic purposes, it is inadvisable to ,use such overlapping windows. Rather, the second window should cover, as shown in Figure 1, the next set of Windependent data points, in this case i = 7-13. In this way the original data set [@(i),f(i)],0 Ii IN , is reduced to a smaller set [ e ( j ) 9 0 ’ ) ] , j = m, ( m + W), (m +2W), .... To be more precise, the new set should be designated [&;U’),& (j;W)] to indicate that it has been generated by using a specific window size W. The number of pairs, q, in the secondary set will satisfy the relation N + 1 = qW + r, where r is the number of remanent primary data points, lying at the end of the data stream, that cannot be utilized. Determination of Rate Consk@ and practical Considerations. Next, the secondary data set [eG;W),e0;W)] is used as input for fitting the modified rate equation
+
-@G;W) = k i ( W ( j ; W ) + k z ( W 2 0 ; W )
(6)
to obtain rate constants and statistical errors (standard deviations) designated ki(W) and ui(W), with i = 1 or 2, to stress that they too refer to a specific window width W. To take advantage of digital smoothing, it is clearly desirable to employ a wide window; the larger the window, the greater will be the noise reduction. On the other hand, it is well-knownlOJl that the LS approach itself will introduce-unless the particular polynomial chosen happens to coincide.with the functional form of @systematic errors in both and which will increase with the size of the window and influence the values of ki(W);the inset in Figure 1 provides a palpable visual example of these systematic errors. Our task, then, is to devise a strategy that combats noise without succumbing to large systematic deviations. It turns out that, if one determines k , and k2 by using not a single window but a range of window sizes, one can employ a ki( W) vs W plot to find the best estimates of ki consistent with the data on hand. At fmt sight, it may appear most reasonable to calculate the slopes by using some fixed “optimum” window. However, it was found that by using the largest possible set of windows, one can optimize the information about the desired parameters and extract the values pertaining to the ideal, infinitesimal window. This concept constitutes a departure from earlier practice and is central to the method described below.
e
Examples ( 1 ) Mixed-Order Decay: (a) Determination of k , and kZ. Suppose that we wish to determine k , in a system in which the bimolecular component dominates. Specifically, let us choose k2 = 5 X 1O8 M-I s-’, eo= 1.O X 10-5 M, and k, = 500 s-l so that k 2 e o / k l= 10. Experimental data sets, [@(i),t(i)], were simulated
The Journal of Physical Chemistry, Vol. 95, No. 26, 1991 10715
Differential Method in Chemical Kinetics
*
J 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 t (ms)
0 0.2
0.1
II
Smoothing and differentiotion Fig. I
Kinetic onalyaia Eq. 6
and determinotion of brat ki(0)
7001.
II
3001 ,
Fig. 3 and Eq.7
-k,(W)
Figure 4. Flow chart summarizing the steps used in obtaining rate constants using the differential method.
I 300
1
:
“:I
- u,( W )
200
. 100 550
-
:
-
k,(W) = a
0 C
450
A 0
100
0.5
300
200
I .o
400
I .5
’w LOW‘
Figure 3. Mixed first- and second-order decay. (a) First-order rate constant k , (SI) as a function of W,the number of data points in a 0% noise. Windows are nonosmoothing window: (-) 1% noise; verlapping. W,is the relative window defined by twindow/t0.5, where 10.5 is the half-life of the reaction. (b) Standard deviation in kl as a function of W. (c) Rate constant k , as a function of Wfor nonindependent or sliding window analysis. (e-)
by using the analytical expression for e, eq 4, and adding Gaussian noise. By n% noise we mean n = (2u/e,,,) X 100, where u is the standard deviation in the noise and emax is the maximum value of e, in this case eo. In the more general case, when there is no convenient analytical expression for the time variation of the concentration, data can be generated by numerical integration of the appropriate rate equation(s). Each set consisted of 1250 data points taken every 0.8 p (=k) over a period of 1 ms in which time the concentration drops to about 10% of its intial value (Figure 2). First-order smoothing, e ( w ) = a pw, was used to obtain the [&;W),e(j;W)] data sets, one for each value of W,which were then fitted to eq 6 to obtain kl( W) and ul( W). Plots of these quantities are shown as solid curves in Figure 3a,b. Since at least three points are required to extract kl(W)and k2(W) from eq 6, the maximum value of Win this case is W,,, = 41 5, the largest odd integer 11250/3. In addition to a W scale for the ordinate, it is also instructive to graduate the axis according to a new parameter, the relative where twindow = WAt and aperture, defined by W, = twindow/t0,5, to,sis the time for the concentration to drop to half of its initial value. This scale gives an immediate measure of the window width in terms of an experimentally measurable parameter ( 2 0 . 5 ) . Returning to Figure 3a, one sees that the k l (W)values appear to be evenly distributed along a straight line with a small positive slope. Clearly, with a single narrow window (say W C 25), the statistical error would be very large indeed, while with a wide window the systematic error might be serious. In fact, in this case, the weak aperture dependence of k,(W) implies that, although
+
only first-order smoothing is being carried out, one need not incur serious systematic error in the rate constants even at large apertures. Reasons for the good performance of the linear approximation and a remedy for dealing with the systematic error will be presented later. The [k,(W),WI set displayed in Figure 3a is fitted to an extrapolation polynomial
+ b W + cw2 + ...
(7)
each k l ( W) value being assigned a weight wl( W) = [ul( W)]-*, where ul(W), the statistical error in kl(W) decreases-as expected from the previous discussion-with increasing aperture (Figure 3b). The best estimate of kl, in which systematic error has been removed, is then kl = kl(0) = a. Once again, the underlining serves to emphasize that the kl(W)values have been subjected to regressional smoothing using eq 7. The degree of the extrapolation polynomial is chosen to be the lowest beyond which the sum of the residuals, C[kl(W)- kl(W)I2,shows no further decrease. In this case, the scatter in the data warrants no more than first-order fitting, kl(W) = a + bW, and k,(O) = 501 s-I is obtained. The uncertainty associated with this value is estimated by forming the weighted average of all the ul(W), (al)= Cwl(W)uI(W)/Cwl(W), and comes out to be 21 s-l. The soundness of this method of error analysis for a single run may be judged from the fact that similar results were obtained when a set of nine independent runs were analyzed: The set [k,(’)(O), kl(2)(0),..., k1(9)(0)]gave an average value of 494 SI with a standard deviation of 20 s-I. The same procedure, when applied to the evaluation of k2, yielded an average value of 5.04 X lo8 M-’ s-l with a standard deviation of 0.03 X lo8 M-I s-l. Notice that the relative statistical error in kl ( = ( u , ) / k , )is about 10 times larger than that in kS,reflecting the dominance of the second-order component (k2eo/kl= 10). By simulating data with noise in the 0.34% range, we also confirmed that the statistical uncertainties in the rate constants were linear with respect to noise. Figure 4 outlines a flow chart summarizing the key steps described above; these are (I) smoothing and differentiation within a window, (11) rate constant and error determination for a range of fmed windows, and (111) removal of systematic error to obtain ubest” kiand their attendant uncertainties. Even when a simple personal computer is used, the whole procedure for obtaining ( u i ) and ki can be carried out in 2 min. (b) Systematic Errors. If it is desired to examine more closely the systematic error introduced by using finite windows, the values of b ( 0 ) obtained from the noisy signal can be used, in conjunction with the reaction scheme, for generating noise-free model data, [ e o ( i ) , t ( i ) ] .Such a set, when the procedures outlined in steps I and I1 above are carried out, provides rate constants kp( W)over a range of Ws, as shown by the dotted curve (Figure 3a). The systematic error in k,(W), brought about by the finite window size,
10716 The Journal of Physical Chemistry, Vol, 95, No. 26, 1991
TABLE I: Relative Percent Systematic Errors in and -&for Two Reactions i C k P (i = 1, 2), under Different Conditions relative 7% systematic error for W, = 1.0 reaction smoothing order, i order e ki
Steel and Naqvi h
3
10
.A
-e
2
1
1 2 3
7.40 7.40 -0.34
4.05 -0.13 -0.13
-0.80 7.70 -0.10
1
1.26 1.26 -0.004
2.11 -0.006 -0.006
-0.83 1.27 0.003
2 3
is then e,( W) = [ky( W) - k?(O)];when subtracted from k,( W), it provides an alternative method of obtaining k,. For example, at W,,, = 415, t l = 25.3 s-I, so the value of kl is then k,(Wmax) - e l ( W,,,,,) = 529.5 - 25.3 = 504 s-I with a standard deviation, ol(Wma,),of 15 s-'. (c) Compensation of Errors. In the example considered, W,,, corresponds to nearly two half-lives, yet first-order smoothing results in only a 5% systematic error in k l . At first sight this may appear surprisingly small. However, for a decay in which -8 decreases with increasing time, it can be shown that the systematic errors in the values of E and at the midpoint of the smoothing window (Figure la) compensate each other when k, is determined. Table I exemplifies this point by considering two simple reactions, viz., -8 = k , P , (i = 1, 2), and comparing, at a relative window aperture of W, = 1, A#, AsE, and Ask,, the relative percentage errors in 0,and k,, where s is the degree of the smoothing polynomial in eq 5. The definitions A@ = 100 X [(E- e)/@ = 100 X [(a,- e)/@] and As& = 100 X - &)/&I imply that one must bear in mind the relations between the expansion coefficients pertaining to different orders: The least-squares expression for a changes in going from first to second order ( a l # a2)but that for p remains unchanged (PI = p2) and conversely in going from second to third order (a2= a3,p2 # 8J. In general, and as = = os, s = 2, 4, 6, .... Consequently, # A2E and A l e = A& and so on. These relations find an application in Table I, where it can be seen that although the error in decreases, as one goes from first-order to second-order remains unchanged, and, consequently, smoothing, the error in the error in k, becomes larger. Only by proceeding to third order can one lower the systematic error in k, below its first-order value. However, this improvement in systematic error is accompanied by a deleterious effect on the statistical error. When a given set of noisy data within a window are fitted to a smoothing polynomial (eq 5), the standard deviations in the concentration and the slope at the midpoint get larger with increasing order, corresponding to greater flexibility in the choice of the fitting parameters a,0, 7,etc. These deviations are reflected in the statistical error associated with k,(W). For the example previously considered, although (ul) = 21 s-l for first-order smoothing, the value doubles in going to third order. Clearly, since the systematic error can be allowed for so well, nothing is to be gained in this case by going beyond first order. (d) Overlapping Windows. It was stated above that successive windows were not allowed to overlap in obtaining a secondary [S(j;W),C?(j;W)] set. Since smoothing is generally carried out by sliding the window along the horizontal axis in steps of one data point,g it may be asked why this procedure is not followed here. Clearly, when a set is formed using overlapping windows, the adjacent points in the set, having been formed from blocks of data which differ only in two [C?(i),t(i)] points, cannot be statistically independent and will, therefore, not provide reliable estimates of statistical uncertainties. For example, if the run analyzed above is subjected to sliding window analysis, the secondary data set of (1 250 - 1 + 415) points generated for W = 415, when fitted by eq 6, yields ~ ~ ( 4 1 5= )0.9 s-l, whereas the three independent points obtained from nonoverlapping windows give ~ ~ ( 4 1 5=) 28 s-'. As it happens, the former process is also more time-consuming. The lack of statistical independence
I
0 0.1
0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 t (ms)
e
-e,
[(os
Ale
e
-e
[a(j;W),e(j;W)]
- @ t ) = -koAoexp(-kot)
+ k 1 8 ( t )+ k 2 g 2 ( t )
(9)
where A. is the initial (known) concentration of species A at time t = 0. In this case eq 9 cannot be integrated to give an expression for 8(t). A typical simulated data set for the indicated parameters is shown in Figure 5. As before, the problem is to extract the rate constant kl for the minor decay channel (k2emax/kl 10). If, as in the first example, one used the total data set, the Occurrence of the peak in the present case would vitiate the LS concentrations and slopes at large windows." To overcome this problem, we were compelled to split the analysis into two parts (shown in the figure) and to deal separately with the regions on either side of the peak. With the occurrence of two or more peaks, as in oscillating reactions or in quantum beats, each region in which concentration varies monotonically with time would become a small fragment of the total data stream. Under these conditions, the reservations expressed by Marchand and Marmet" would apply also to kinetic (12)Press, W.H.; Flannery, B. P.; Vetterling, W. T. Numerical Recipes: The Art of Scientific Computing Cambridge University Press: Cambridge, U.K., 1986;Chapter 10. (1 3) Noggle, J. H. Physical Chemistry on a Microcomputer; Little Brown: Boston, 1985;Chapter 8. (14)Shoefeld, P. S.;deVoe, J. R.Anal. Chem. 1976,48, 403R.
0.8 1
-
The Journal of Physical Chemistry, Vol. 95, No. 26, 1991 10717
Differential Method in Chemical Kinetics I
600
100
-
-
3 w
_______________-----‘,tort
.
-- 300 - 200
***.
**
4
-400.
--- .....
c
Y
-900 .
100
-
-
‘i’229
193
287
q-4
’
231
1
289,
- 1400
J Figure 6. Slopes as a function of window aperture for nine independent runs having the same parameters as for the run displayed in Figure 5. J is the index of the aperture center and ( 2 J t 1) is the aperture width.
problems, and our strategy of employing a wide range of window sizes would become unfeasible. (a) Determination of Initial Slope. In the first region, successively larger windows, all starting at the origin and centered on i = 1, 2, 3, ..., are taken, and as before (see inset to Figure I), LS slopes associated with the midpoint of these windows are determined. Figure 6 shows such plots for nine independent runs with the same parameters as those used for generating the curve shown in Figure 5 . Notwithstanding the inevitable noise at small windows, the initial slope (=koAo)can be ascertained by gradually opening the window up to the maximum permissible width (which is determined by the position of the peak) and polynomial fitting of the data set after discarding the error-prone low-J values. Finally, extrapolation to J = 0 gives the desired initial slope and illustrates once again that-contrary to what one might have expected-ideal, “zero-window” values are best estimated by using wide rather than narrow windows. (b) Determination of k l and k2. The later portion of the second region, dominated by decay and hence by k , and k2, is used as before to generate secondary sets for a range of windows. In analogy to eq 6, these q pairs of secondary data points are then fitted to = -koAo exp[-kotU)l + kl(W)W;W)+ kdW20‘;W) (10)
[e(j),e0’)]
-m;W)
However, in this case, since the right-hand side of eq 10 is not a polynomial, a nonlinear least-squares procedure (e.g., Simplex) has to be used to obtain k l (W)and k2(W). To remove systematic error, ki( W) is then plotted against W to obtain ki(0), the best values of ki. The plots look quite like that shown in Figure 3a. Nine independent runs yielded k , = 514 f 60 s-l and k2 = (4.99 f 0.14) X 108 M-’ s-l, so that the standard errors increase by about a factor of 3 as compared to the simple mixed-order decay. The starting point, istart,of the segment used for estimating kl(W) is, it will be noticed, quite arbitrary; the numbers quoted above were obtained with is,,,, = 300. One would be justified in thinking that, by not using a larger portion of the second region, one is being profligate with one’s data. To investigate the (dis)advantages of enlarging the segment by moving closer to the maximum, we analyzed noise-free data, choosing ish,, = 100, 200, 300. Two noteworthy points emerge from this study, the results of which appear in Figure 7. One notices, in the first place, that the “dependence of the systematic error, moderate when is,,,, = 300 or 200, becomes steep for istart= 100. This behavior can be traced to the fact that, for istart= 300, the first term on the right-hand side of eq 1 0 has become much smaller than the two terms representing the decay of the product; the aforementioned compensation of errors therefore comes into play here again and leads to a significant reduction of systematic error. The second prominent feature of Figure 7 is the stepwise change in k , ( W ) for is,,,, = 100. To grasp this behavior, let us note first that one now has 1250 - 100 = 1150 primary data points for determining k , ( W). On referring to Figure 7 and recalling that q is now given by the relation of 1150 = q W + r, one sees that for 193 IW I 229, q = 5 ; since r = 185 and 5 for W = 193 and 229, respectively,
P
3,
3,
33
30
1 2 3 4
1 0.954 0.967 0.978 0.984 0.988 0.991 0.993 0.994 0.995
1 1.315 1.728 2.098 2.406 2.655 2.855 3.016 3.147 3.255
1 1SO5 2.341 3.396 4.641 6.062 7.651 9.401 1 1.308 13.368
1 1.168 1.369 1.518 1.624 1.700 1.755 1.796 1.828 1.853
5 6
7 8 9 10
and since the terminal points are those for which compensation works best, the value of k1(193) will be subject to a larger systematic error than kl(229). This accounts for the slightly upward trend in the nearly horizontal part of the data for which q = 5 . When W takes the value 231, q drops to 4 and r jumps to 226; since r is even larger than that for W = 193 and q smaller, the systematic error in kl(W) increases abruptly at W = 231. The slight improvement thereafter (until W = 287) is again due to the steady reduction in r. The staircase behavior of kl(W), which becomes noticeable only when systematic error becomes large, can be understood along these lines. It should be added here that, in the presence of noise, these steps will be partially or wholly submerged in statistical fluctuations; moreover, extrapolation to W = 0 still leads to a reliable value for k l ( 0 ) . (3) Biophysical Example. Proteins are often idealized as ellipsoids of revolution and their rotatory diffusion coefficients determined by a combination of optical, spectroscopic, and hydrodynamic For a protein, randomly labeled with a suitable spectroscopic probe, the decay of flash-induced anisotropy, r ( t ) , can be expressed as17 i ( t ) r(t)/r(O)= 0.4 exp(-t/&) + 0.4 exp(-t/&) + 0.2 exp(-t/&) (1 1) where the rotational correlation times 4i ( i = 1, 2,3) are functions of D , and D,,, the perpendicular and parallel rotational diffusion coefficients, respectively, which in turn are related to the axial ratio p = a/b. Determination of reliable diffusion coefficients by fitting the observed time course of anisotropy to a triexponential curve is by no means an easy mission.I6 However, we have found that for dealing with a prolate ellipsoid of revolution the following differential technique is rapid as well as reliable. Although absolute magnitudes of 4idepend upon the volume (V) of the ellipsoid, which at this stage is now known, relative values q$ = 4i/4s(where 4sis the rotational correlation time of a sphere with the same volume as the ellipsoid), which depend only on p , can still be calculated.I6 For convenience, some values of d j for a prolate ellipsoid are shown in Table 11. Let -do-] be (15) Scheraga, i-I. A.; Mandelkern, L. J. Am. Chem. Soc. 1953.75, 179. (16) Small, E. W.; Isenberg, 1. Biopolymers 1977, 16, 1907. (17) Acuiia, A. U.;Gonzalez-Rodrigues, J.; Lillo, M. P.; Razi Naqvi, K. Biophys. Chem. 1987, 26, 63.
Steel and Naqvi
10718 The Journal of Physical Chemistry, Vol. 95, No. 26, 1991
”O
L 9
0 I.I
3
2
1
0
Time (in units of
4
$0)
Figure 8. Anisotropy decay for a prolate ellipsoid with axial ratio p = 1, 3, 4, 6, 8, 10. The time scale is in units of &, where -& (s-I) is the slope of the initial decay.
400
800
1200
1600
2000
Sample number ( i )
Figure 10. Decay of flash-induced anisotropy in a randomly labeled prolate ellipsoid of revolution (of axial ratio p = 3) with an admixture of 3% Gaussian noise.
Once 4o becomes known, the effective volume V of the protein can be calculated from the equation”
0.28
c
>.
0.23
Here k, T, and 1 are the Boltzmann constant, absolute temperature, and solvent viscosity, respectively. F , and FIIare the rotational Perrin factors1*J9which can also be expressed in terms of the now known p .
a
e
c 0
.-
2
0.18
0.13
p 1
3
5
7
9
Axial ratio, p
Figure 9. Anisotropy at time 2& for prolate ellipsoids as a function of axial ratio p .
identified with the initial slope in the anisotropy decay: 4o = 441 4& 2&. It is convenient to plot model data for 3 using a dimensionless time scale in which the unit of time is do. When this is done for several values of p , all curves have the same initial slope, -1, but differ significantly for t > 4o (see Figure 8). A plot of P(t = 2 4 ~ against ~) p appears in Figure 9. When one is dealing with actual experimental data, simulated in Figure 10, the initial slope is obtained by the method exemplified above. Once again, it proved beneficial to use windows encompassing the total data range, precluding only the low aperture ‘noisy” outputs. After &, has been obtained in this fashion, the P(i) vs i curve is used to obtain the best value of P at some time that is a multiple of 40, say 240. This is done by considering a set of P(i) values around 240 and carrying out regression interpolation. The axial ratio can then be found by referring to Figure 9. Nine independent data sets like those shown in Figure 10 gave a 2% standard deviation in 4o while the value of the axial ratio was 2.9 f 0.2. It is noted that this procedure makes use of two constraints (40and p ) to define the curves, which stands to reason, for one has two independent parameters, D , and Dl,,at one’s disposal. The final step, redundant here since we are working with simulated data, would consist of computing the theoretical anisotropy curve a n d verifying t h a t t h e residuals a r e distributed evenly; systematic deviation(s) would evidently imply departure from the assumed model.
+
+
Concluding Remarks Unless one is in possession of an instrument or a transducer with a delta-function response, convolution-let us remembercannot but distort the input signal. The experimenter’s aim is, therefore, to defeat this distortion by an appropriate choice of the transducer or to recover the original input by carrying out a deconvolution of the output. Now, LSP smoothing or differentiation of equidistant data points amounts essentially to a convolution;”’ so far as we can see, the resulting distortion has been accepted as a fait accompli, and attempts are generally made to minimize it by using narrow windows. After conducting a thorough examination of the applicability of LSP smoothing and differentiation to three kinetic problems, we have come to the counterintuitive conclusion that there arise many instances when the use of wide windows is to be recommended rather than deprecated. We have found that-by an apt choice of window sizes, careful consideration of the attendant systematic errors, and judicious extrapolation-one can simulate an ideal, infinitely narrow window without sacrificing the advantages associated with wide windows. In our approach, extrapolation to a delta-function window performs a role analogous to deconvolution. The principles propounded here, though they have been illustrated by considering three kinetic examples, should prove serviceable in a wider context. Acknowledgment. C.S. acknowledges support from the Research Corporation; K.R.N. is grateful to the University of Trondheim for leave of absence and to the Norwegian Research Council for Science and Humanities (NAVF) for munificent monetary emoluments. (18) Perrin, F. J . Phys. Radium 1934, 5 , 497. (19) Koenig, S. H. Biopolymers 1975, 14, 2421
