Semianalytical Methods To Determine First-Order Rate Constant

Semianalytical methods for the determination of first-order rate constant distributions are discussed. In general the expression for the overall decay...
0 downloads 0 Views 327KB Size
Download PDF
Langmuir 1997, 13, 961-969

961

Semianalytical Methods To Determine First-Order Rate Constant Distributions† Luuk K. Koopal,*,‡ Maarten M. Nederlof,§ Willem H. Van Riemsdijk,| and Peter A. Barneveld‡ Department of Physical and Colloid Chemistry and Department of Soil Science and Plant Nutrition, Wageningen Agricultural University, NL-6703 HB Wageningen, The Netherlands Received November 20, 1995. In Final Form: July 23, 1996X Semianalytical methods for the determination of first-order rate constant distributions are discussed. In general the expression for the overall decay function is an integral equation based on the local decay function and the distribution function. The analytical expressions of the distribution function obtained by inversion of the integral equation are a series of derivatives of the overall decay function. The higher the order of the approximation, the more derivatives are required. The most common method is based on the Laplace transform technique; newly derived are the coefficients for the third-order method. The second method is due to Schwarzl and Staverman, it provides a general scheme to find a distribution function from an integral equation and has been used here to obtain an expresssion for the rate constant distribution. The advantage of this method is that it provides a weighting function that maps the true distribution into its approximation; i.e., it visualizes the quality of the approximation. The third method is newly developed. It uses an approximation of the local decay function to solve the integral equation for the distribution function. The local decay function approximation or LODA method provides a physical interpretation of the approximation involved in obtaining the distribution by showing the function that approximates the true local decay function. The three methods are compared on the basis of two synthetic data sets, one composed of exact data and one of nonexact data. For nonexact data presmoothing of the data is required. For this purpose a smoothing spline technique is applied in which the smoothing parameter is obtained objectively by generalized cross validation in combination with physical constraints. It is shown that the newly developed LODA-G2 method, which needs the first and second derivative of the overall decay function, combines a good resolution with a relatively low sensitivity to experimental error.

Introduction Many processes both in industrial practice and in natural environments are affected by the heterogeneity of the substrate. In order to describe the effect of the heterogeneity on a process, information on the heterogeneity is required. In the case of adsorption or binding studies, such information can be obtained by studying the sorption process with special probes followed by an analysis of the results for heterogeneity. Methods to investigate the heterogeneity of binding reactions have been developed for both equilibrium1-8 and kinetic conditions.9-14 Very recently Cerefolini and Re15 have indicated the similarities between the two types of analyses. The present paper focuses on the heterogeneity analysis of kinetic data. * Corresponding author: e-mail address, [email protected]. † Presented at the Second International Symposium on Effects of Surface Heterogeneity in Adsorption and Catalysis on Solids, held in Poland/Slovankia, September 4-10, 1995. ‡ Department of Physical and Colloid Chemistry. § Present address: KIWA NV Research and Consultancy, P.O. Box 1072, 3430 BB Nieuwegein, The Netherlands. | Department of Soil Science and Plant Nutrition. X Abstract published in Advance ACS Abstracts, February 15, 1997. (1) Koopal, L. K.; Vos, C. H. W. Langmuir 1993, 9, 2593-2605. (2) LumWan, J. A.; White, L. R. J. Chem. Soc., Faraday Trans. 1991, 87, 3051-3062. (3) Von Szombathely, M.; Brau¨er, P.; Jaroniec, M. J. Comput. Chem. 1992, 13, 17-32. (4) Jagiello, J.; Schwarz, J. A. J. Colloid Interface Sci. 1991, 146, 415-424. (5) Jagiello, J. Langmuir 1994, 10, 2778-2785. (6) Nederlof, M. M.; Van Riemsdijk, W. H.; Koopal, L. K. J. Colloid Interface Sci. 1990, 135, 410-426. (7) Nederlof, M. M.; Van Riemsdijk, W. H.; Koopal, L. K. Environ. Sci. Technol. 1992, 26, 763-771. (8) Nederlof, M. M.; Van Riemsdijk, W. H.; Koopal, L. K. Environ. Sci. Technol. 1994, 28, 1037-1047. (9) Stanley, B. J.; Topper, K.; Marshall, D. B. Anal. Chim. Acta 1994, 287, 25-34.

S0743-7463(95)01054-7 CCC: $14.00

For complicated systems the overall reaction rate is determined by a series reactions that each have their own order and rate constant. To obtain the rate constant distribution from the overall reaction rate, the order of the reactions has to be known. With complicated reactions the reaction order is often unknown, and an assumption has to be made with respect to the order; generally a firstor pseudo-first order reaction is assumed. Under the assumption of a series of first-order reactions, the overall decay function is essentially a multiexponential decay. Such a decay function appears in a large variety of studies.16,17 In principle, the rate constant distribution function can now be obtained without further assumptions about the rate constant distribution itself. Yet the determination of the distribution is not trivial because small variations in the overall reaction rate due to experimental error may lead to large spurious variations in the distribution function.17,18 In order to avoid these problems, special numerical techniques have to be used. References 9, 16, 17, and 19 can be consulted for discussion and examples. Three advanced numerical techniques, based on regularized least squares to find an optimal (10) Shuman, M. S.; Collins, B. J.; Fitzgerald, P. J.; Olson, D. J. In Aquatic and terrestrial humic materials; Christman, R. F., Gjessing, E. T., Eds.; Ann Arbor Science Publishers: Ann Arbor, MI, 1983; pp 349370. (11) Olson, D. L.; Shuman, M. S. Anal. Chem. 1983, 55, 1103-1107. (12) Olson, D. L.; Shuman, M. S. Geochim. Cosmochim. Acta 1985, 49, 1371. (13) Nederlof, M. M.; Van Riemsdijk, W. H.; Koopal, L. K. Environ. Sci. Technol. 1994, 28, 1048-1053. (14) Lavigne, J. A.; Langford, C. H.; Mak, M. K. S. Anal. Chem. 1987, 59, 2616-2620. (15) Cerefolini, G. F.; Re, N. J. Colloid Interface Sci. 1995, 174, 428440. (16) Zimmermann, K.; Delaye, M.; Licinio, P. J. Chem. Phys. 1985, 82, 2228-2235. (17) Provencher, S. W. Comput. Phys. Commun. 1982, 27, 213-227. (18) Noble, B. The state of the art in numerical analysis; Academic Press: London, 1979; pp 915-966. (19) Provencher, S. W. Comput. Phys. Commun. 1982, 27, 229-242.

© 1997 American Chemical Society

962

Langmuir, Vol. 13, No. 5, 1997

solution for the rate constant distribution, have been discussed recently.9 The methods based on the maximum entropy principle minimize artifactual information optimally,9 but the procedures are fairly complicated. For more simple solutions to the problem, further approximations have to be made. A much more simple analysis results if it is assumed, a priori, that the distribution function is of a certain type,20,21 but this is more a fitting exercise than a true distribution analysis. Another way to reduce the problem is to assume that the number of reactions is very limited. Examples of this approach can be found in studies of Schechter22,23 and Mark and Rechnitz.24 This procedure is, however, only satisfying if the overall reaction rate is indeed determined by only a limited number of reactions.25 In the present paper we will discuss a fourth group of methods that provide analytical expressions for the distribution function in the form of a series of derivatives of the relation describing the overall decay function (e.g., the amount of probe left at time t). Partly the origin of these methods goes back to the Laplace transform technique.26,27 These methods have the advantage that they show in a direct way the relation between the overall decay function and the heterogeneity without making assumptions about the rate constant distribution function. Another advantage is that they provide easy access to that part of the distribution function that corresponds to the “window” of measured data points. Connors28 used the first derivative of the overall decay function to approximate the rate constant distribution, whereas Olson et al.10-12 used the first two derivatives. In order to obtain reliable derivatives, presmoothing of the data is required.8,12,13 The method of Olson has been used in several studies regarding metal ion binding.10-14,29,30 A matter of confusion in these studies is the amount of smoothing of the data. In the present paper we will (I) outline the basic assumptions and equations involved in the kinetic heterogeneity analysis, (II) briefly discuss the approximations based on the Laplace transform technique (LT), (III) show that the general method of Schwarzl and Staverman31,32 leads to equivalent results and provides a way to visualize the quality of the various approximations, and (IV) present a new approach in which the the local decay function (the decay function for a homogeneous reaction) is replaced by an approximation that allows for a direct inversion of the integral equation.33 The latter methods are closely related to the LIA methods for the determination of equilibrium affinity constant distributions6,7 and will be called the LODA (local decay approximation) methods. One of these methods has been used to reanalyze the data of Olson.13 (20) Rate, A. W.; McLaren, R. G.; Swift, R. S. Environ. Sci. Technol. 1992, 26, 2477-2483. (21) Rate, A. W.; McLaren, R. G.; Swift, R. S. Environ. Sci. Technol 1993, 27, 1408-1414. (22) Schechter, I. Anal. Chem. 1992, 64, 729-737. (23) Schechter, I.; Schro¨der, H. Anal. Chem. 1992, 64, 325-329. (24) Mark, H. B.; Rechnitz, G. A. Kinetics in analytical chemistry; Interscience: New York, 1968; pp Chapter 7. (25) Marshall, D. B. Anal. Chem. 1989, 61, 660. (26) Doetsch, G. Theorie und Anwendung der Laplace Transformation; Springer: Berlin, 1937. (27) Widder, D. V. The Laplace Transform; Princeton University Press: Princeton, NJ, 1946. (28) Connors, K. A. Anal. Chem. 1975, 47, 2066-2067. (29) Cabaniss, S. E. Environ. Sci. Technol. 1990, 24, 583-588. (30) Gutzman, D. W.; Langford, C. H. Environ. Sci. Technol. 1993, 27, 1388-1393. (31) Schwarzl, F.; Staverman, A. J. Physica 1952, 18, 791-798. (32) Schwarzl, F.; Staverman, A. J. Appl. Sci. Res., Sect. A 1953, 4, 127-141. (33) Nederlof, M. M. Analysis of binding heterogeneity. Ph.D. Thesis, Wageningen Agricultural University, Wageningen, The Netherlands, 1992; p 243.

Koopal et al.

Last but not least, special attention will be given to an objective way of presmoothing the data. In order to illustrate the methods and problems involved in the kinetic heterogeneity analysis, we will focus the treatment on the example of the dissociation of metal ionsubstrate complexes. The study of metal ion complexes is of special relevance for the mobility of metal ions in the natural environment. Normally the available data sets for these systems are limited, and for such data sets analytical methods are very well suited. Most of the cited studies regarding the kinetic analysis, refs 9-14, 20, 21, 29, and 30, also concentrate on this type of system. The methods are however general and can be applied to obtain the rate constant distribution as long as the reactions are (pseudo) first order and the reverse reaction is so slow that it can be neglected. Basic Assumptions and Equations The kinetics of the dissociation of ligand-metal complexes in solution can be studied by monitoring the amount of metal ions that is released as a function of time according to the dissociation reaction

MLi f Li + M

(1)

The amount of released metal can be followed by adding a sink (S) that instantaneously takes away the ion M that is released

M + S f MS

(2)

The process described by eq 1 can be followed indirectly by measuring the amount of MS formed, when the dissociation of the MLi complex is the rate-limiting step. The amount of MS can, for instance, be measured accurately by a photometric technique.12,14 For a correct result the rate constant for dissociation of the MS complex has to be orders of magnitude lower than that of the dissociation of MLi. When it is assumed that the dissociation of the complex MLi, into M and Li, is a first-order reaction, the amount of complex MLi, left at time t (Ci(t)), can be written as t

Ci(t) ) Ci(0) e-ki

(3)

where ki is the rate constant of the dissociation of MLi and Ci(0) is the amount of MLi at t ) 0. For a heterogeneous ligand mixture Ci(t) is called the local decay function because each ligand type i has its own rate constant ki. The overall dissociation of a mixture of ligands is given by the summation of the decay contributions of all the individual site types, resulting in an overall decay function n

C(t) )

Ci(0) e-k ∑ i)1

t

i

(4)

For a continuous range of dissociation constants eq 4 can be rewritten as an integral equation

C(t) ) C(0)

∫0+∞e-kt g(k) dk

(5)

where g(k) is the normalized rate constant distribution function and C(0) is the normalization constant that expresses that at t ) 0 all ligands have to be occupied equally with M. Practically speaking, C(0) is the total amount of ML at t ) 0. Since the value of the rate constant may vary several orders of magnitude, it is appropriate to switch to a normalized rate constant distribution function f that is a function of log k instead of k

Rate Constant Distributions

Langmuir, Vol. 13, No. 5, 1997 963

∫-∞+∞e-kt f(log k) d log k

C(t) ) C(0)

(6)

Equalizing the respective cumulative rate constant distribution functions leads to f(log k) ) 2.3k g(k). Note that when C(0) is not known, as is often the case in practical situations, the corresponding non-normalized distribution functions

G(k) ) C(0) g(k)

(7a)

F(log k) ) C(0) f(log k)

(7b)

or

have to be used. Absence of knowledge of C(0) is for the present discussion no obstacle, and for the moment we assume that C(0) is unity. After C(t) is measured, the rate constant distribution function f(log k) can in principle be obtained by inversion of eq 6. As mentioned in the Introduction, several ways exist to solve eq 6 for the distribution function. Below we will discuss a series of methods that result in analytical expressions for the distribution function. Before doing so we first present the data set used to test the methods. To illustrate the methods, a synthetic data set based on a known distribution function is used. This has the advantage that the calculated distribution can be compared with the true distribution. A synthetic data set C(t) is calculated using eq 6 and a multi-Gaussian distribution n

f(log k) )

wi N(µi,σi) ∑ i)1

(8)

where N(µi,σi) is the Gaussian distribution with mean µi and standard deviation σi and wi is the fraction of the individual Gaussians. We constructed a distribution with four peaks (n ) 4) as shown in Figure 1a. The parameters are given in Table 1. The synthetic decay function was obtained by performing the integration of eq 6 numerically using eq 8 as distribution function. The result, expressed as C(t)/C(0), is shown in Figure 1b. A smooth monotonically decreasing function is obtained without distinct steps, whereas in the distribution, sharp peaks are present. This illustrates how difficult it will be to retrieve the distribution function. To be able to study the effect of a small random error in the data, a second data set is generated which is identical to the first one except for the presence of a small error of σ ) 1 × 10-3 in the overall decay function. As a consequence, the relative error varies from 0.5% for high values of the decay function to 20% for low values of C(t). Secondly, the number of points was limited to 120 (20 points per decade) to illustrate that in practice it is often difficult to obtain a very large data set. Methods To Determine Rate Constant Distributions Laplace Transform Approach. The overall decay function C(t) is given in eq 5 as the Laplace transform of the rate constant distribution function G(k). The rate constant distribution function can thus be obtained by taking the inverse Laplace transform of C(t). Following Doetsch26 and Widder,27 the inverse Laplace transform can be written as

G(k) ) lim

mf∞

( )( ) ( ) (-1)m m!

m k

m+1

∂mC(t) ∂tm

(9) t)m/k

Figure 1. Rate constant distribution function (a) and the corresponding overall decay function (b) used to illustrate the methods. See Table 1 for the parameters. Table 1. Synthetic Data Set, Parameters of the Distribution Function mean, µi width, σi fraction, wi

peak 1

peak 2

peak 3

peak 4

-2 0.0686 0.1

-1.25 0.0686 0.1

-0.5 0.0686 0.1

1.0 0.6142 0.7

In practice, however, arbitrary high-order derivatives of C(t) are not available since C(t) generally is a measured function consisting of only a limited set of data points. Still eq 9 can be used to obtain approximations for G(k). These approximations involve linear combinations of derivatives of C(t) up to the mth derivative and can be referred to as mth order approximatons, Gm(k). Changing to the distribution function F(log k) instead of G(k) and from t to log t does not change the order of the approximation, and the general approximation formula can then be written as m

Fm(log k) )

∑

v)1

( )

av(0.43)v-1

∂vC(t)

∂ log tv

(10)

t)γ/k

where Fm denotes the mth order approximation of the true F. Note that eq 10 allows us to evaluate the derivatives at t ) γ/κ rather than only at t ) m/k as required in eq 9. The purpose of this generalization becomes clear in the next section. The coefficients av are found by substitution of finite values of m in eq 9 and writing out the derivatives in terms of d log t for γ ) m and are given up to the third-order approximation in Table 2. Briefly these approximations will be called LT1, LT2, and LT3, where the number indicates the order of the approximation. The LT1 and LT2 methods have been used before,10-12,28 the LT3 method not. The quality of the different approximations is shown in Figure 2 where a comparison is made between the calculated distribution functions and the true distribution function. As expected

964

Langmuir, Vol. 13, No. 5, 1997

Koopal et al.

Table 2. Summary of the Coefficients av and γ of Eq 10 for the Distribution Function Fm(log k) According to the Various Approximations, with m (e3) the Order of the Approximation, LT ) Laplace Transform Approximation, SS ) Approximation According to Schwarzl and Staverman, LODA ) Local Decay Function Approximation m)1 a1 γ1

LT1

LT1

LODA-1a

-1 +1

-1 +1

-1 +0.7

m)2 a1 a2 γ2

LT2

LT2

LODA-L2a

LODA-G2

-1 +1 +2

-1 +1 +2

-1 +1 +1.6

-1 +1.48 +2.2

m)3 a1 a2 a3 γ3

LT3

SS3a

LODA-G3b

-1 +1.5 -0.5 +3

-1 +1.56 -0.66 +3.3

-1 +0.69 +0.79 +1.1

a Optimal solution, it is also possible to arrive at the LT solution. Optimal solution, it is only possible to arrive at the LT solution by using a combination of LODA functions.

b

the best result is obtained with the third-order approximation. In general the distributions are flattened compared to the true distribution. The three sharp peaks appear as a shoulder in the first-order approximation and as strongly smoothed peaks in the other two approximations, compared to the true distribution. The top of the wide peak is slightly displaced in all three approximations. Schwarzl and Staverman Approach. In their study of viscoelastic behavior, Schwarzl and Staverman31,32 have discussed, in general, the inversion of integrals like eq 6 that contain a kernel (in the present case exp(-kt)) in combination with a distribution function. They have shown that for a measured function C(t) and any type of kernel, the distribution function can be approximated by eq 10. The approximations obtained with the inverse Laplace transform, for the case that the kernel equals exp(-kt), can be considered as special solutions of the general eq 10 for γ ) m and the av values as shown in Table 2. According to Schwarzl and Staverman a general method to obtain the values of the coefficients av and γ is by substituting the equation for C(t), our eq 6, into eq 10. For the present situation this results in

∫-∞+∞K(k′,k) f(log k′) d log k′

Fm(log k) ) C(0)

(11)

where k′ is a dummy variable and where the kernel K(k′,k) equals m

K(k′,k) )

∑

v)1

( )

av(0.43)v-1

∂ve-k′t

∂ log tv

(12)

t)γ/k

The function K(k′,k) is a combination of derivatives of the local decay function with respect to log t. The better this function resembles a Dirac-Delta function, the closer Fm(log k) will approximate F(log k). For any m the coefficients av and γ are determined by demanding K(k′,k) to be as sharp as possible and to have one and only one maximum at kt ) γ. The kernel K(k′,k) should preferably be positive, otherwise the distribution may contain negative parts and at extreme values of k′, K(k′,k) should

Figure 2. Comparison of the distribution functions obtained with the Laplace transform technique with the true distribution function: (a) first-order method; (b) second-order method; (c) third-order method.

become zero. The higher the value of m, the sharper the kernel can be made. For the first- and second-order approximation the demands to K(k′,k) lead to the result as obtained with the first- and second-order LT approximation.33 For the thirdorder approach several possibilities exist. When the position of the maximum of the kernel is fixed at kt ) 3, the results are equivalent to the LT3 method. A sharper kernel can be obtained by enlarging a2 and lowering a3. The sharpest kernel is obtained for a2 ) 1.56 and a3 ) -0.66. The position of the maximum is then found at kt ) 3.3. Figure 3 shows the kernels for the four approximations. The higher the order of the approximation is, the sharper the kernel is. The dashed line in Figure 3 represents the optimized solution that gives a slightly better result than the third-order LT approximation. The kernels are given for the case k ) 1. Different values for k just shift the curve along the log k′-axis. The advantage of the work of Schwarzl and Staverman is that K(k′,k) provides a measure of the quality of the approximation involved. Moreover, K(k′,k) can be seen as a weighting function that maps the true distribution

Rate Constant Distributions

Langmuir, Vol. 13, No. 5, 1997 965 Table 3. Relation between the Coefficients (ri and βi) in the LODA Eq 13 and the Coefficients av of Eq 10 for Fm(log k) R1 ) 0 a1 ) -1 a2 ) 0 a3 ) 0

Figure 3. The function K(k′,1), see eq 12, for the LT1, LT2, LT3, and SS3 approximation. The dashed line corresponds with the third-order method for which the coefficients were optimized numerically following Schwarzl and Staverman (SS3).

function into its approximation. In principle it is therefore possible to use the kernel for a deconvolution of the calculated distribution function Fm(log k) and to obtain a more accurate representation of F(log k). Local Decay Function Approximation (LODA). A totally different solution to the problem of finding F(log k) has been suggested by Nederlof.33 The strategy in this approach is to replace the kernel in eq 6, i.e., the local decay function exp(-kt), by a suitable mathematical function that makes it possible to solve eq 6 analytically for f(log k). This type of approach was applied successfully for the analysis of equilibrium binding data, and it was called in that case local isotherm approximation, LIA.6,7,13 Simple functions that can be used to replace the local decay function are, for instance, a step function and a linear function. A general and flexible approximation of the local decay function, which comprises also the step and the linear function and has the desired property that it leads to an analytical solution for f(log k), is the socalled LODA function. The LODA function is composed of two equations

(ktγ ) kt )R ( ) γ

β1

ξ1,LODA ) 1 - R1 ξ2,LODA

β2

2

for kt e γ for kt > γ

(13a) (13b)

that replace the true local decay function, exp(-kt). In eq 13, γ is the “breakpoint” which divides the approximation in its two parts and Ri and βi are constants. It will appear that γ is also equivalent to γ in eq 10. For R1 ) 0 and R2 ) 0 the LODA approximation is a simple step function with the step situated at kt ) γ. For R1 ) 1 and R2 ) 0, ξ1 is either a simply linear (β1 ) 1) or a slightly curved (β1 * 1) function that reaches its zero value for kt ) γ. In the more general case that 0 < R1 < 1, the full eq 13 have to be used. In order to make sure that at kt ) γ a continuous function is obtained, the following constraints should be obeyed

1 - R1 - R2 ) 0 R1β1 - R2β2 ) 0 i.e., for a given set of R1 and β1, R2 and β2 are fixed. The general expression for the approximated distribution function that can be derived by using the LODA eq 13 in eq 6 (rather than exp(-kt)), is again provided by eq 10. However, the LODA approach gives only solutions for m e 3. For the derivation of this result, ref 6 should be consulted. The coefficients av and γ are in this case

0 < R1 < 1 a1 ) -1

a2 )

β 2 - β1

R1β12 + R2β22 1 a3 ) 2 R1β1 + R2β22

R1 ) 1 a1 ) -1 a2 ) 1/β1 a3 ) 0

determined by the coefficients Ri and βi of the LODA function eq 13. The relations between av and Ri and βi are shown in Table 3 for R1 ) 0, for 0 < R1 < 1, and for R1 ) 1. The values of Ri, βi, and γ either are chosen or can be obtained by optimization. Relevant cases are discussed below. Approximating the local decay function by a step function (R1 ) R2 ) 0) leads to R1 ) -1 and corresponds for γ ) 1 to the first-order LT approximation. The advantage of the present method is that it shows the physical meaning of the approximation involved. The LODA method also allows us to locate the step by a best fit criterion. By applying a best fit criterion on a log kt scale results in γ ) 0.7. The two approximations are shown in Figure 4a on a log kt scale. The effect of the shift in γ on the calculated distribution function is that the wide peak shifts to its true position. The approximation of the exponential decay function for R1 ) β1 ) 1 is a linear function that goes from ξ ) 1 at kt ) 0 to ξ ) 0 at kt ) γ. The coefficients in eq 10 become R1 ) -1 and R2 ) +1. For γ ) 2 the expression for the distribution function corresponds to the secondorder LT approximation. The optimal solution is determined by the least sum of squares method applied on a log kt scale. The best fit of the local decay function was found for γ ) 1.6. In Figure 4b the two approximations are illustrated on a log kt scale. The calculated distribution function for γ ) 1.6 shows the right position of the broad peak and slightly too low values for the peak positions of the small peaks. Instead of approximating the local decay function with a linear function, it is also possible to use R1 ) 1 and a finite value for β1. By applying the least sum of squares method on a log kt scale, the best fit of the local decay function was obtained for β1 ) 0.68 and γ ) 2.2; see Figure 4c. The corresponding expression for the distribution function is of the second order with coefficients R1 ) -1 and R2 ) 1.48. The coefficient R2 is considerably larger than that obtained for the second-order LT approximation but about equal to that observed for the third-order LT approximation. As shown in Figure 5a, in the calculated distribution the broad peak is resolved rather well and the small peaks are in the right place. We will call this method LODA-G2, where the G stands for general and the 2 emphasizes that it is a second-order method. Using the general approximation of the local decay function (eq 13 with 0 < R1 < 1) and optimization of all parameters (Ri,βi,γ) for a log kt scale leads to γ ) 1.1, R1 ) 0.68, R2 ) 0.32, β1 ) 0.77, and β2 ) 1.64. The approximation is shown in Figure 4d and indicated as LODA-G3. For these parameter values the resulting distribution function is of third order and the coefficients are a1 ) -1, a2 ) 0.69, and a3 ) 0.79 (see Tables 2 and 3). This result relies less strongly on the second and more strongly on the third derivative of the overall decay function than the third-order LT approximation. In the calculated distribution function the peaks are slightly better resolved than with the LT3 approximation; see Figure 5b.

966

Langmuir, Vol. 13, No. 5, 1997

Koopal et al.

Figure 4. Approximations of the local decay function for the LODA methods and the true local decay function on a logarithmic kt scale: (a) step functions on a kt scale (R1 ) 0, R2 ) 0), γ ) 0.7 (LODA1) and γ ) 1 (local decay function corresponding with the LT1 method); (b) linear functions on a kt scale, R1 ) 0, β1 ) 1, γ ) 1.6 (LODA-L2), and γ ) 2 (local decay function corresponding with the LT2 method); (c) LODA-G2 function, R1 ) 1, β1 ) 0.68, γ ) 2.2; (d) LODA-G3 function R1 ) 0.68, R2 ) 0.32, β1 ) 0.77, β2 ) 1.64, and γ ) 1.1.

Figure 6. Weighting functions K(k′,1) see eq 12, for the LODAG2 and LODA-G3 approximations.

Figure 5. Comparison of the true distribution with the distribution functions obtained with (a) the LODA-G2 and (b) the LODA-G3 methods.

Knowledge of the coefficients av of the LODA-G2 and the LODA-G3 distribution functions also allows us to calculate the kernel K(k′,k), according to eq 12, that

corresponds with each of these two distribution functions. For both the LODA-G2 and -G3 method the weighting functions that map the true distribution into its approximation are shown in Figure 6 for the case k ) 1. Comparison of these functions with those obtained for the LT2 and LT3 method (see Figure 3) shows that especially the LODA-G2 weighting function is sharper than the LT2 function. A disadvantage of the two methods is that the weighting functions become negative for some values of log k outside the peak region. These negative values may lead to irrealistic slightly negative values in the calculated distribution function. In general we may say that the LODA methods are conceptually simple: the physical approximation involved can be seen directly from the agreement of the LODA eq 13 with the true local decay function. However, results show that even a close agreement may lead to considerable deviations from the true distribution if narrow peaks are present. The LODA-G2 method that does not need a third

Rate Constant Distributions

Langmuir, Vol. 13, No. 5, 1997 967

derivative is the most promising. Although the peaks are somewhat less distinctly resolved than with the thirdorder LODA method, the narrow peaks are in the right position and slightly better resolved than with the LT2 method. In addition the broad peak is considerably better resolved than with the LT2 method. The LODA-G2 result is comparable to the LT3 result; however, the LODA-G2 method will be less sensitive for experimental error, since it does not need a third derivative. Data Treatment Derivatives of experimental data, which are prone to errors, are hard to obtain. To deal with this problem, a smoothing spline can be used to smooth the data.34,35 The smoothing spline method can be applied provided that the measured data can be represented as

yi(xi) ) g(xi) + i

(14)

where g(xi) is the sought value of the unknown function g(x), yi(xi) is the measured value at x ) ξi, and  is assumed to be a normally distributed random error in y only. The smoothing spline procedure previously applied to study affinity constant distributions,13 can be adapted for use with the present problem. A quintic spline was used because in that case derivatives up to the third, needed for the calculation of the distribution functions, are continuous. In the smoothing spline method the amount of smoothing is determined by the smoothing parameter; this quantity determines the trade-off between the goodness of fit and the smoothness of the curve. The higher the smoothing parameter p, the more weight is given to the smoothness of the curve and the less weight to the goodness of fit. The larger the experimental error, the larger the smoothing parameter will be. The determination of the correct smoothing parameter is extremely important. However, often an arbitrary value of the smoothing parameter is used (refs 11, 29, and 30). This may lead to large differences in the obtained distribution functions when the analysis is performed by different users. To prevent confusion, an objective way to establish the smoothing parameter should be used. We propose to use generalized cross validation (GCV)36 in combination with physical constraints. The GCV criterion provides a first estimate of the smoothing parameter on the basis of the random error implicitely present in a series of data points. Although the GCV technique has been proven to result in the correct value for the smoothing parameter for very large data sets, this estimate of the smoothing parameter is sometimes too small. This is due to the facts that the error may not be distributed purely random and that the data set is of limited size. To correct for this, knowledge of the physical characteristics of both the overall and the local decay function can be used to improve the estimate of the smoothing parameter; see also ref 13. The following physical constraints apply. Since the first derivative of the local decay function with respect to t is negative and the distribution function is always positive, it follows that the first derivative of the overall decay function with respect to t should always be negative or equal to zero. From the same reasoning it follows that the second derivative of C(t) should always be positive or equal to zero and that its third derivative should be negative or zero. These constraints can be adjusted to (34) Silverman, B. W. J. R. Statist. Soc. B 1985, 47, 1-52. (35) Woltring, H. J. Adv. Eng. Software 1986, 8, 104-107. (36) Craven, P.; Wahbla, G. Numer. Math. 1979, 31, 377-403.

Figure 7. The log(GCV) function for the synthetic data set on which a random absolute error of σ ) 1 × 10-3 on C(t) is imposed (panel a). log(GCV) is plotted as a function of the smoothing parameter, log p. The minimum occurs at log p ) -4.2. In panel b, the smoothing spline representation of the overall decay function is shown.

the heterogeneity analysis method to be used; for instance, for a second-order approximation the constraints on the first and second derivative are used. In practice the procedure to obtain p is as follows. First the GCV minimum and the corresponding value of p are calculated. Subsequently the spline corresponding to the GCV smoothing parameter is tested for the constraints. When one or more constraints are not obeyed, the smoothing parameter is enlarged until all the constraints are fulfilled. Confidence intervals of the distribution functions can be calculated on the basis of the variance in the spline result using Monte Carlo simulations.34 For the results presented below we calculated the variance of the distribution on the basis of 100 realizations of the spline. Calculated Distributions To illustrate the method in the case of a small error in the data, the data set with a known error of σ ) 1 × 10-3 in C(t) is used. First the GCV minimum is calculated. The GCV function is shown in Figure 7a for a range of values of the smoothing parameter p. A minimum is found at log p ) -4.2. In the present case the value of log p ) -4.2 obeys all the constraints. The smoothing spline function for log p ) -4.2 representing the overall decay function corresponding with the data set is shown in Figure 7b. The bars denote the confidence interval. It follows that only for high values of log t, where the relative error is the largest, is the uncertainty significant. The distribution functions calculated using the various approximations are shown in Figures 8 and 9. The uncertainty in the overall decay function for high values of log t is reflected in a large uncertainty in the distribution

968

Langmuir, Vol. 13, No. 5, 1997

Koopal et al.

Figure 9. LODA distribution functions for the synthetic data set containing a small absolute error in C(t): (a) LODA-G2 method; (b) LODA-G3 method. The error bars are derived from the spline function through C(t).

see Figure 9b. The LODA-G2 method (Figure 9a), which does not need a third derivative, results in much smaller error bars and a slightly better result than that obtained with the second-order Laplace method. Also here, the two small peaks with the low rate constants are lost and appear as a “foot” of the distribution function. Discussion

Figure 8. LT distribution functions calculated for the synthetic data set containing a small absolute error in C(t): (a) LT1; (b) LT2; (c) LT3. The distributions are calculated on the basis of the spline function shown in Figure 7b. The error bars indicate 66% confidence intervals.

for low values of log k. For the first-order LT approach (LT1, Figure 8a) a reliable result is obtained over almost the entire log k range, since only the first derivative is needed. For the second-order LT approximation (LT2, Figure 8b) the range where the resulting distribution is quite uncertain shifts to higher values of log k. The lower two small peaks are not resolved. The three small peaks appear as a “foot” of the broad peak, which is recovered very well. The result of the LT3 approximation is shown in Figure 8c. Despite the fact that the third-order Laplace approach has in principle the best resolving power, the distribution function obtained in practice is not better than that shown in Figure 8b; details in the distribution function cannot be interpreted due to the large error bars. For methods relying more strongly on the third derivative (e.g., the LODA-G3 method) the situation will be worse;

We are now able to compare the results obtained for exact synthetic data with those obtained for data on which a random error was imposed. It follows that the firstorder methods are not capable of resolving narrow peaks in a distribution. Even for exact data, narrow peaks are only resolved as very weak shoulders. On the other hand for the third-order methods that are capable of resolving the narrow peaks more distinctly, it applies that the error bars are too large to allow for the interpretation of small peaks, even for a small experimental error. Thus, only for highly accurate data can higher order methods be used. For the second-order methods the peak with the smallest rate constant is “covered” by large error bars. On comparison of the LT2 method with the newly developed LODA-G2 method, it is found that the resolving power of the latter method is better than that of the LT2 method. The broad peak is reproduced very well with the LODAG2 method and is of much better quality than that obtained with the LT2 method. Compared to the result obtained with the LT3 method, the broad peak is recovered only slightly worse with the LODA-G2 method, but the error bars are much smaller. In conclusion we might say that the LODA-G2 method has the best potentials for use in practice.

Rate Constant Distributions

An illustration of results obtained in practice with the LODA-G2 method can be found in ref 13. The results of the whole procedure, that is, a spline representation of the overall decay function and the rate constant distribution function calculated from the spline, can be used for further analysis; see e.g., Lavigne and Langford.14 On the basis of the rate constant distribution it is also possible to look for analytical expressions for the overall decay rate. An analytical expression for the overall decay function can be obtained, for instance, by using a γ distribution function for log k.37 Cerefolini and Re15 have shown that Elovich type equations can also describe the overall decay function. In this case the underlying distribution function is a series of normalized exponential functions. Conclusions The overall decay function of a set of first-order reactions can be analyzed to obtain the rate constant distribution. A method used previously is the inverse Laplace transform (LT) technique. In analogy with the affinity spectrum method and the local isotherm approximation (LIA) method for the analysis of equilibrium data, two alternative methods of analysis are presented in this paper. The analogue of the LIA method is called the local decay function approximation (LODA) method. The resulting expressions for the distribution function are in all cases a series of derivatives of the overall decay function. As a consequence of the fact that approximations have to be made to arrive at an analytical expression for the distribution function, the calculated distribution function is an approximation of the true distribution function. For wide and smooth distribution functions the methods approximate the true distribution function very closely. However, a series of narrow peaks is difficult to recover with the methods. In principle the higher order methods are best suited to reproduce details in the distribution function. The two alternative methods both provide a means of visualizing the quality of the approximation involved. Schwarzl and Staverman have shown that for each approximation a weighting function can be obtained that (37) Shkilev, V. P.; Bogillo, V. I. React. Kinet. Catal. Lett. 1993, 50, 355-361.

Langmuir, Vol. 13, No. 5, 1997 969

maps the true distribution to its approximation. In principle this weighting function makes it possible to further deconvolute the calculated kinetic spectrum and to obtain a better estimation of the true spectrum. The LODA methods provide a physical interpretation to the approximation involved; the LODA equation (13) shows how the local decay function is approximated for each method. In order to apply any of the three methods for kinetic heterogeneity analysis to experimental data which include an experimental error, a treatment of the data is necessary. For this purpose an adapted smoothing spline technique is suggested. An objective estimate of the smoothing parameter is obtained by applying the generalized cross validation (GCV) technique in combination with a number of physical constraints related to the overall decay function. This procedure is essentially an extension of the method by Olson et al.,11 who also used a smoothing spline but took an arbitrary value of the smoothing parameter. In order to be able to obtain details in the distribution functions, very accurate data are needed. For ordinary quality data wide and relatively smooth distributions of rate constants can be reproduced relatively well with all methods. Narrow distributions or narrow peaks superposed on a distribution cannot be found unless the quality of the data is extremely good. The secondorder methods, although in principle less good than the third-order methods, are most suited to reveal narrow peaks, as they are less sensitive to experimental error than the third-order methods. The newly developed LODA-G2 method is the most promising both with respect to its resolving power and with respect to the sensitivity to experimental error. Acknowledgment. This work was partially funded by the Netherlands Integrated Soil Research Programme under Contract Number PCBB 8948 and partially by the EU Environmental Research Programme on Soil Quality under Contract Number EV4V-0100-NL (GDF). The GCV software of H. J. Woltring, called GCVSPL, which was used for the calculations presented in this paper, is available via [email protected] and [email protected]. LA9510543