)V1ESI.J RETURNING MATERIALS: PIace in book drop to LIBRARIES remove this checkout from All-[3!!!L. your record. FINES wiII ——— be charged if book is returned after the date stamped below. APPLICATIONS OF COMPUTER-BASED MODELING METHODS TO ANALYTICAL CHEMISTRY By Peter Dale Wentzell A DISSERTATION Submitted to Michigan State University in partial fulfillment of the requirements for the degree of DOCTOR OF PHILOSOPHY Department of Chemistry 1987 ABSTRACT APPLICATIONS OF COMPUTER—BASED MODELING METHODS TO ANALYTICAL CHEMISTRY By Peter Dale Wentzell Several applications of computer-based modeling methods to analytical chemistry are described. Areas of application include reaction-rate methods of analysis, response surface modeling, and liquid chromatography. Two reaction-rate methods of analysis are introduced which attempt to minimize the effect of between-sample variations in the pseudo-first-order rate constant. The first method uses two rate measurements during the reaction and is therefore called the two-rate method. A model for the propagation of random errors for this method is developed and shown to be valid. The second method introduced is based on the extended Kalman filter, a recursive analog to nonlinear least squares fitting. Chemical systems are used to show that both methods are less susceptible to between—run variations in the rate constant than traditional methods. The Kalman filter is also shown to enhance the capabilities of differential kinetic methods, those reaction- rate methods which take advantage of differences in the rate constant to determine two or more analytes simultaneously. A comparison of traditional and optimized reaction-rate methods of analysis, including the two—rate method, is presented for reactions following first-order kinetics. The comparison focuses on the susceptibility of seven reaction-rate methods to rate constant variations and their precision in the absence of those variations. To aid in this comparison, models are developed which illustrate the dependence of each method on systematic errors in the rate constant. Response surface modeling is an important area of many scientific disciplines and a new approach to this technique, based on the Kalman filter, is described. The new method has numerous advantages over traditional approaches to linear least squares modeling, including speed, simplicity, and adaptability. The new method is applied to three chemical systems to demonstrate its utility. A theoretical model for recycle heartcut chromatography (RHC) is developed to evaluate its potential. The model shows that RHC is only useful for obtaining resolution improvements if long heartcut loops are used and/or dispersion in the heartcut loop can be reduced. The utility of RHC for enhancing chromatographic selectivity is demonstrated, however. Also presented in this work are a program for nonlinear curve fitting and a summary of methods for generating normally distributed random numbers for simulation. To my parents, Carl and Rochelle, who planted the seed and nurtured its growth. iv ACKNOWLEDGEMENTS No man is an island, and that is particularly true for‘ the work presented here. In my five years at Michigan State University many people have contributed to my progress and I will take this opportunity to acknowledge as many as possible. My deepest appreciation goes to Dr. Stan Crouch, of course, who provided guidance and direction while allowing me to pursue my own avenues of inquiry. His constant encouragement and support were invaluable. I would also like to acknowledge Tom Doherty, whose enthusiasm for all fields of science inspired much of the research which is presented here, and some which is not. Special thanks are also due Kim Ratanathanawongs, who was always there for a walk around campus and an encouraging word when things got hectic. There are numerous others worthy of mention (in no particular order): Helen Archontaki, for her constant good humor; Max Hineman, who will always be a physical chemist and existentialist at heart; Dean Peterson, for his friendship over the years and many lunchtime discussions; Bob Harfmann, whose stopped—flow wizardry is unequalled; Mark Victor, John Stanley, and Keith Trischan, three of the finest programmers I have ever met; Paul Kraus, who is one of the mOIt competent researchers I know; Cheryl Stults, who managed to have a baby and still finish ahead of the rest of uI; Nelson Herron, who will be V remembered for his perseverance and brownies; Kristine Kurtz, for breakfast at Theio’s and tolerating late night experiments with remarkable good humor; Stephen Medlin, who bakes a hell of a cake even if he does prefer to program in Pascal; Mayda Lopez-Nieves, for the Puerto-Rican avocado among other things; Larry Bowman, who is willing to tackle any task from lasers to plumbing to carpet—laying; Jimmy Wu, who shares my love of teaching; Pat and Marguerite Wiegand, who inspired us all; Dr. Chris Enke, for his vision; Drs. Miltos and Stella Karayannis, who I’ve just begun to know; and Dr. Adrian Wade, with whom a research dynasty is hopefully just beginning. I would also like to thank the departmental support staff who have been so helpful during my stay at MSU, especially Dr. Tom Atkinson for his computer genius, Ellen the graduate secretary, Bill in graphics, and Ron, Tom, and Scott in the electronics shop. In addition, there are those outside the immediate MSU community who deserve credit. These include my family (especially Kendra for all her letters), Margaret-Anne Vaughan for her support over the years, and Dr. Louis Ramaley for his guidance while I was an undergraduate. To all of these individuals, and any I may have overlooked, I offer my sincerest gratitude. "You read the book, you turn the page, You change your life in a thousand ways, The dawn of reason lights your eyes, With the key, you realize, To the kingdom of the wise. " - Alan Parsons Project "Nothing Left to Lose" vi TABLE OF CONTENTS CHAPTER I. III. LIST OF TABLES LIST OF FIGURES INTRODUCTION A REACTION-RATE METHOD OF ANALYSIS INSENSITIVE TO BETWEEN-RUN CHANGES IN RATE CONSTANT GENERAL CONSIDERATIONS A. Two-Rate Method B. Propagation of Random Errors ,,,,, EXPERIMENTAL SECTION A. The Exchange Reaction B. The Molybdophosphate Reaction C. Instrumentation D. Data Analysis RESULTS AND DISCUSSION A. pH Effects“ B. Thermal Effects“, C. Propagation of Random Errors CONCLUSIONS LIST OF REFERENCES COMPARISON REACTION-RATE METHODS OF ANALYSIS FOR SYSTEMS FOLLOWING FIRST—ORDER UNETICS EXPERIMENTAL SECTION vii PAGE xi 30 33 CHAPTER PAGE RESULTS AND DISCUSSION 33 A. Susceptibility to Rate Constant Variations 33 B. Precision with Invariant Rate Constant 39 C. Additional Factors 43 CONCLUSIONS 45 LIST OF REFERENCES 46 IV. THE KALMAN FILTER FOR REACTION-RATE METHODS OF ANALYSIS 49 GENERAL CONSIDERATIONS 51 A. The Kalman Filter 51 B. Minimizing the Effect of Rate Constant Variations 57 C. Differential Kinetic Methods 61 EXPERIMENTAL SECTION 65 A. The Molybdophosphate Reaction 65 B. Simulated Data 66 C. Data Processing 66 RESULTS AND DISCUSSION 67 A. Minimizing the Effect of Rats Constant Variations 67 B. Differential Kinetic Methods ‘78 C. Additional Applications 85 CONCLUSIONS 86 LIST OF REFERENCES 88 V. MODELING CHEMICAL RESPONSE SURFACES WITH THE KALMAN FILTER 90 THEORY ....... ' 92 A. Relationship to Deterministic Least Squares 92 B. Model Definition 93 EXPERIMENTAL SECTION~ pg 96 A. Analytical Systems , 96 3- Computational Aspects .. 98 viii CHAPTER PAGE RESULTS AND DISCUSSION 99 A. Seven-Variable FIA System 99 B. Five-Variable FIA System 102 C. Two-Variable Polarographic System 107 D. Advantages of the Kalman Filter Approach 111 E. Relationship to Simplex Optimization M... 113 CONCLUSIONS 114 LIST OF REFERENCES 115 VI. RECYCLE HEARTCUT CHROMATOGRAPHY: THEORY AND EVALUATION 117 EXPERIMENTAL SECTION 122 A. HPLC Apparatus 122 B. Separations 125 RESULTS AND DISCUSSION 125 A. Resolution Improvement 126 B. Limit Recycle Number 127 C. Length Limitations 130 D. Evaluation of Resolution Enhancements 132 E. Selectivity Enhancements 136 CONCLUSIONS 138 LIST OF REFERENCES 140 VII. CONCLUSIONS AND FUTURE PERSPECTIVES 142 APPENDIX 1. XYFIT - A PROGRAM FOR NONLINEAR CURVE FITTING 145 LIST OF REFERENCES 161 APPENDIX II. GENERATING NORMALLY DISTRIBUTED RANDOM NUMBERS FOR SIMULATION 163 LIST OF REFERENCES 171 ix TABLE 2.1 2.2 2.3 3.1 3.2 3.3 LIST OF TABLES Summary of limiting cases of the two-rate method.m__~_.___ Effect of pH variations on the initial rate method (IRM) and two-rate method (TRM). Slope of calibration curves for three reaction-rate methods at three temperatures. Equations for relative errors arising from rate constant variations for selected reaction-rate methods. Comparison of results Obtained by various reaction—rate methods under conditions Of large changes in temperature. Based on calibration curves generated at 26 °C. Precision of various reaction-rate methods under conditions of an invariant rate constant. PAGE 11 21 24 36 38 41 FIGURE 1.1 2.1 2.2 3.1 3.2 4.1 4.2 4.3 4.4 LIST OF FIGURES Classification of modeling methods in science. Contour surface showing theoretical precision of the two-rate method as a function of measurement time. Contours are represented in units of the RSD in the rate measurement at t. Diagonal lines extending from the origin represent actual measurement situations in that the ratio Dt/ t1 remains constant regardless of changes in the rate constant. Comparison of theoretical and experimental random errors for the two-rate method as a function of measurement time. The dashed line represents the theoretical model. Theoretical errors expected in the concentration determined by various reaction-rate methods as a function of the relative deviation of the rate constant from its nominal value. Instrumental elements of systems for reaction-rate methods of analysis. Qualitative representation of factors influencing the update of the of the Kalman filter state vector after a new measurement. Typical evolution of state parameter estimates with the extended Kalman filter in accordance with the model At=Ar(1-e"“)+B. Comparison of calibration curves generated with the extended Kalman filter (EKF) method and the initial rate method at three temperatures. Fits of molybdophosphate reaction curves at three temperatures to a first-order model by nonlinear least squares ([P]=10 mg/mL). xi PAGE 14 26 35 44 56 68 71 72 FIGURE 4.5 4.6 4.7 4.8 4.9 4.10 4.11 5.1 5.2 5.3 Fits of molybdophosphate reaction curves at three temperatures to a first-order model by the extended Kalman filter ([P]=10 mg/mL). Evolution of final absorbence estimates at three temperatures. ' Evolution of the determinant of the information matrix at three temperatures. Typical simulated data for differential kinetic analysis by the Kalman filter in accordance with the model: A; = AMI-ed“) + Ann-6"“) + B. Values used were AI=1.0, A2=0.1, k2/k1=3, 1% noil_...e._m.,._... Results of differential kinetic analysis by the Kalman filter as a function Of rate constant ratio. A1:100’ A3:001. Results of differential kinetic analysis by the Kalman filter as a function of rate constant ratio. Ax=0.1, Az=1.0. Comparison of actual (points) and predicted (curved solid lines) errors for the Kalman filter approach to differential kinetic analysis as a function of rate constant ratio. Comparison of coefficients determined by the Kalman filter and singular value decomposition (SVD) for the seven-variable FIA system. Response curves for the determination of isoprenaline by FIA. (a) Experimental univariate responses as a function of coil length (other variables at four levels, as indicated by symbols). (b) Experimental univariate responses as a function of ferricyanide concentration. (c) Contour map of the two-variable response surface in accordance with the derived second- order model. Contours represent 10 mm differences in response from 100-180 mm. Other variables are held at their optimum values. Three-dimensional display of points Obtained from the simplex Optimization of the polarographic determination of uranium (VI). The labels on the data refer to the order in which the points were obtained in the simplex procedure. xii PAGE 73 74 76 79 80 81 83 101 105 108 FIGURE 5.4 6.1 6.2 6.3 6.4 6.5 Response surface generated by quadratic modeling of the data in Figure 5.3. Response is in mA, concentrations are in mM. Instrumental configurations for recycling methods in HPLC. HPLC system and valve configuration used to implement recycle heartcut chromatography. Valve switching sequence used in recycle heartcut chromatography. Theoretical maximum resolution improvement attainable with RHC as a function of the length and geometry of the heartcut loop. Selectivity enhancement in the separation of napthalene and biphenyl with RHC after two passes through the column Comparison of actual (histogram) and theoretical (smooth curve) distributions Obtained with the inverse cumulative normal distribution approximation. Results are based on 10,000 samples. xiii PAGE 110 119 123 124 133 137 168 CHAPTER I Introduction "Man tries to make for himself in the fashion that suits him best a simplified and intelligible picture of the world; he then tries to some extent to substitute this cosmos of his for the world of experience, and thus to overcome it. " Albert Einstein "Ideas and Opinions" The word "model" has been defined as "a tentative description of a theory or system which accounts for all of its known properties". The idea of a model, as expressed in the opening quotation, is central to chemistry and science as a whole, as it is a human construct which bridges the gap between pure mathematics and physical reality. It is a ubiquitous concept, and one that allows us to understand, simplify, and make predictions about the physical world. The mathematical formulation of the model implicitly contains a description of the system under study and serves as a means to test our view of physical phenomena. Although the subjects of the work presented here may at first appear unrelated, they are connected by a single thread in that they are all, to some degree, modeling methods. The term "modeling method" is by no means well-defined and is used in its broadest sense here. Some clarification is in order, however. Modeling methods may be classified into two major categories: those methods which attempt to derive mathematical models for the description of a physical system (fundamental modeling methods) and those which apply these models to 1 obtain information about a physical system (applied modeling methods). This distinction is indicated in Figure 1.1. Of the fundamental modeling methods, there are essentially three classes, as represented in the figure. Purely mathematical formulations are perhaps the most common and most useful type of model. Integrated rate laws are an example of this type of model. Simulations are another type of fundamental modeling method, and these include Monte Carlo methods and random walk simulations. Finally, there are empirical modeling methods, which are useful for predictive purposes, but are not derived from a fundamental physical model. There are no doubt cases where the distinction between these becomes . blurred, but the segregation is useful nevertheless. The applied modeling methods can also be divided into three categories. The goal of these methods is not to obtain a model, but rather to extract some information about the system under study (or the model itself) by applying the model. An important modeling method here is the testing of models with experimental data. This may be done to (a) test the validity of the model, or (b) determine whether or not the data conform to the physical assumptions underlying the model. Curve fitting is the primary example of such testing methods. A second type of applied modeling method is the extraction of model parameters (eg. fluorescence lifetimes, rotational constants, etc.) from the physical system. Although curve fitting may also be used here, the context of its usage is quite different. In this case we assume the model is valid and obtain information which will allow us to better characterize the system under study. Finally, models can be used for comparison and Prediction. A given model can be used to predict the behavior of a 60:33 5 23.305 mam—ODOR mo :oSooammmSO .HA Gun—NE 223: 3%.. as .> .5 .5 53...... 51¢ .>_ .5 51¢ .> .>_ ..= .5 As .5 .5 2.5 E 23 Execute 958... «5:33.50... 3052‘ 55.23: 58.5553 . BEE—tom .23: Box—.5050: casesEfi Botaem 5 5 «35¢: 55:28.2 3:53. meofimz 3:25: 358853 «85.2 3.255 physical system under a certain set of conditions (eg. when the temperature is raised) and this may be compared with results predicted under other conditions. Likewise, models for different systems may be compared as an alternative to direct comparison of the systems themselves. Computers have played a central role in modeling methods. First, they may be used to provide numerical solutions to models for which there is no simple analytical form, as will be demonstrated in the following chapters. They are essential for simulation methods, where they allow a large number of iterations of the basic process to be performed. Finally, in empirical modeling and and curve fitting computers reduce the burden of calculation tremendously. All of the modeling methods represented in Figure 1.1 are employed in one way or another in the chapters which follow. For each of the classes shown in the figure, the chapters to which the method applies are given in parenthesis. A chapter by chapter synopsis of the of the methods employed' is now given. Chapter II describes a new reaction-rate method of analysis, the two-rate method, which is insensitive to between—run variations in the rate constant for systems following first-order or pseudo-first-order kinetics. This method makes use of the first—order model to extract the rate constant from a two-point measurement. Changes in the rate constant which occur between samples can be accounted for in this manner. Because the method involves two rate measurements at different times, it was important to develop a model for the propagation of random errors in the rate measurements. The resulting expression is relatively complex and depends on the times of the two rate measurements. No analytical solution to the resulting equation was possible, so numerical solutions were Obtained with the aid of a computer. These are used to develop an error surface for the two-rate- method, which is useful in determining the weaknesses of the method and selecting Optimum measurement times. The validity of the error model is examined by comparing predicted errors to those actually Observed for different measurement times. The model was found to perform very well in spite of some general assumptions which were made in its derivation. In Chapter III, several reaction-rate methods of analysis are compared, particularly from the perspective of susceptibility to between- run variations in the rate constant and errors in the absence of those variations. Models are developed relating systematic errors to changes in the rate constant. These provide insights into the weaknesses and strengths of each method. The validity of the models is verified with experimental results, which support the conclusions drawn. Chapter IV considers the application of the Kalman filter to reaction-rate methods of analysis. The Kalman filter may be considered a modeling method since it requires the definition of a measurement model and a model for the propagation of states in a time-variant system. It allows the extraction of model parameters in much the same way as least squares fitting, but does so in a recursive fashion. The Kalman filter is successfully applied to the problem of between-run rate constant variations and its potential for use with differential kinetic methods is evaluated with simulated data. In Chapter V, the Kalman filter is applied to the problem of response surface modeling. This is an empirical modeling method which allows the exploration of response surfaces for which there is no known analytical description. The Kalman filter is shown to produce the same results as traditional least squares methods, but has a number Of advantages. Three chemical systems are modeled in this study. A theoretical model for recycle heartcut chromatography (RHC) is developed in Chapter VI. While RHC has been used in our laboratory for some time, this is the first theoretical investigation into its potential. The model developed allows the resolution enhancements which can be obtained with this method to be examined under realistic physical conditions. It does, however, require numerical solution with the aid of a computer. The limitations of RHC and the experimental barriers which must be overcome become clear through examination of the model. Nonlinear curve fitting is an important part of modeling methods and Appendix I describes a simplex-based curve fitting program, XYFIT, which is referred to in Chapters 11, III, and IV. While the concepts employed in this program are not original, the particular implementation is unique, and the widespread use of the program justifies the brief description. Finally, Appendix II considers the only example of simulation modeling in this work. Methods for modeling experimental noise with a normal distribution are reviewed, and a specific implementation that was used in Chapter IV is described. The ability to model experimental noise accurately is an important consideration in the evaluation of new methods of signal processing or data analysis. It is hoped that the modeling methods described in the work presented here demonstrate the importance of these methods for solving problems in analytical chemistry. CHAPTER II A Reaction-Rate Method of Analysis Insensitive to Between-Run Changes in Rate Constant "In one moment I’ve seen what has hitherto been Enveloped in absolute mystery, And without extra charge I will give you at large A lesson in Natural History." - Lewis Carroll "The Hunting of the Snark" Reaction-rate methods are now widely used in analytical and clinical chemistry, as evidenced by the numerous review articles [1-5], books [6,7] and two international conferences [8,9] devoted to the subject. Despite their popularity, rate methods of analysis are subject to many errors and experimental uncertainties that do not affect equilibrium-based methods. Among these are changes in first- or pseudo-first-order rate constants due to variations in experimental conditions. These changes may occur during the analysis of a single sample or between samples. It is the latter problem, between-run variations in the rate constant, which is addressed here. Carr [10] has shown how relatively small changes in pH or temperature can adversely affect the precision of rate methods. Other experimental parameters, such as ionic strength and reagent concentration, are also important. Thus, traditional, single-component reaction-rate methods, such as the fixed-time [1,10,11], variable-time [1,12-14], and derivative [15-19] methods, require careful control of the experimental conditions. Several 8 methods have been suggested to alleviate this constraint [20524], but most of these, while reducing the rate constant dependence, do not eliminate it. In this chapter, a new approach to reaction-rate methods of analysis is proposed [25]. This method is based upon making two rate measurements during the course of a first-order or pseudo-first-order reaction and will therefore be referred to as the two-rate method. The new method is, in theory, completely independent of between-run variations in the rate constant. It does not compensate for changes in the rate constant within a single reaction and requires that the reaction follow first-order or pseudo—first-order kinetics. The principles of the new method are presented and it is evaluated with two test reactions. A model for the propagation Of random errors is also developed and tested. A comparison of reaction-rate methods, including the two-rate method, appears in the fOllowing chapter and elsewhere [26]. GENERAL CONSIDERATIONS The calculated response parameter for the two-rate method is the two-rate parameter (TRP). This section begins with a derivation of the generalized two-rate parameter, which is independent of the rate constant. A model for the propagation of random errors with the method is then developed. A. The Two-Rate Method Consider a first-order or pseudo-first—order reaction in which analyte A reacts with reagent B to give product P, with rate constant k, as given in Equation 2.1. At any time t, the rate of this reaction, R: is given by, R: = (d[P]/dt): = k[A]o 6"“ (2.2) where [A]o is the initial concentration of the analyte. Rate measurements made at two times during the course Of the reaction, t1 and t2, can be ratioed to obtain k, as shown in Equation 2.3. Rn __ k( t2- t1) _ kAt m - e - e th 1n(Rt1/}?t2) (2.3) 1‘: At Substitution of this result into eq 2.2 and rearrangement yields an expression for the initial concentration of the analyte which is independent of the rate constant. [A]o 1n(£::/A;t2)(RtI/Rt2)tl/At = TRP(t1,t2) (2.4) 10 Note that this derivation assumes that k is the same at times t1 and ta. The expression in eq 2.4 is the generalized two-rate parameter (TRP(t1,ta)). It is equally applicable to the more common case where rates are measured in terms of some linear response function, such as the absorbance of a product. In this case, the rates are expressed in terms of the rate of change of the response function, and a proportionality exists between the two-rate parameter and the initial concentration of the analyte. ‘ There are at least two limiting cases which lead to simplified versions of eq 2.4. The result of the first simplification will be referred to as the simplified two-rate parameter. It requires that an initial rate measurement be made (TRP(0,t)). The second simplification leads to the extended two-rate parameter, which has the constraint that t1=At=t (TRP(t,2t)). The equations for each of the three cases mentioned here (simplified, extended, and generalized) are summarized in Table 2.1. Each has its own merits, and the selection of an appropriate form depends on the situation under consideration. The two-rate method described here should be applicable to any analytical reaction that follows first-order or pseudo—first-order kinetics. These include enzyme reactions which approximate first-order behavior in the determination of a substrate [27]. B. Propagation of Random Errors The two-rate method involves the measurement of two rates rather than one and requires the evaluation of a non-linear function. Therefore, an important consideration is the effect of random errors in 11 um; I had I ad sad I .36 n J 51... .56 I had I 3. I J .3 I :3 I I 3: I .5 883 uses—2:308 38mm... .5: u . so £5. a a. 8...: .2 .353: u 535. u S. L34?» . a mar? I 5.30 + a: Idamiu «W E u Ana—camp. 63:82:.» 3‘: 3 E «a a...» A: \c v r «A: I 5.8» + a: I «:5» fl 3%..“— Q a Emanuaag popes»: lbWulll a . e «as» + a: I 3 AleuuWE I :6va poem—a8? .8519, 2582 he «38.8.. :38 1.3505... 83 95:8: .1930:— 39793 0.5 no condo U535: no 5.3.:an .H.N 039—. 12 the rate measurements. An error model for the two-rate method is now presented. If errors in the measurement of time are assumed to be negligibly small, the variance in the two-rate parameter is given by the expression, 6(TRP) 2 6(TRP) 2 2 = 2 + 2 2.5 “TR? ( 6R“ )nnan" ( OBI: Rack“ ( ) where (nu)3 and (naP are the variances in the rate measurement at times t1 and ta, respectively. To develop a reasonably simple equation from this generalized expression, we assume equal variances for all rate measurements; that is, (¢m)3 = (anP = (003. This assumption should be approximately correct for many situations and serves as a starting point for examining the propagation of errors. For the experimental systems and rate measurement techniques that we employed, the equal variance assumption was found to be a reasonable approximation. Based on this assumption, a general expression for error propagation can be derived. The derivation, which is lengthy but straightforward, leads to a final expression for the relative variance in the two-rate parameter, and hence the measured concentration. This is given as Equation 2.6. Puma/[Aw = amz/(TRW = me. + At - 1)2 + emu - tail "8. _ (2.6) e2(At)2 R.2 13 As eq 2.6 shows, the relative variance in the two-rate parameter is a function of the times of the two rate measurements. It is important to note that all of the time measurements in eq 2.6 have been adjusted to units of 1' (=1/k) in order to remove the explicit dependence of the equation on the magnitude of the rate constant. Also note that the relative variance in the two-rate parameter is expressed in terms of the relative variance of the rate, as measured at time T. This was chosen as a convenient point of reference for comparison with other reaction- rate methods [26]. Equation 2.6 can be simplified somewhat for the two limiting cases of the two-rate method previously discussed. The simplified forms are listed in Table 2.1 with the corresponding two-rate expression. Figure 2.1 is a contour surface showing the propagation of random errors as a function of the measurement times according to the model developed. The contour lines represent loci of equivalent relative standard deviation in the two-rate parameter. These are expressed in units of the relative standard deviation of the rate measurement at t: 1’. The "lowest" contour, for example, is labelled with a "1" and encompasses the region in which the RSD of the two-rate parameter is less than that of the rate measured at 1'. Subsequent contours delineate regions of 2 RSD(R1'), 3 RSD(R1), and so on. These contours were generated through a rearrangement and numerical solution of eq 2.6, since no analytical solution could be found. As shown in Figure 2.1, results are most imprecise when the time interval between measurements is short (At is small) or when measurements are made at extended time intervals (t1+At is large). This is expected, since the relative error in the rate measurements increases 14 2.0 I I T l f; 4 \ .. . .. / 2 o / ,2 1.0 ‘ c 3 1 c C 2 4:- o.o ; l. } : 0.0 1.0 2.0 3.0 4.0 At (In units Of 1/k) Figure 2.1. Contour surface showing theoretical precision of the two- rate method as a function of . measurement time. Contours are represented in units of the RSD in the. rate measurement at T. Diagonal lines extending from the origin represent actual measurement situations in that the ratio At/ t1 remains constant regardless of changes in the rate constant. 15 as the reaction nears completion, and the relative error in the log term increases as the measurements are - made close together. The convergence of the contour lines as At -> 0 cannot be regarded as reliable since the assumption of small relative errors in At becomes invalid at this point. The contour surface in Figure 2.1 allows selection of optimum measurement times to minimize the effect of random errors in the rate measurements. This selection does not reduce to simply locating the minimum on the surface, since we are anticipating changes in the rate constant that will alter t1 and At, which are expressed in units of 1/k on the contour map. Therefore regions of the surface must be sought where the gradients are small so that random errors in the rate measurements are not detrimental to the results, even when changes in the rate constant occur. The diagonal lines through the contours in Figure 2.1 represent actual measurement situations in that the ratio tI/ At remains constant. The simplified case (t1 = 0) corresponds to a line coincident with the abscissa, and the extended case (t1 = At) is represented by the diagonal line shown with a slope of unity. Other measurement situations can be represented by similar straight lines extending from the origin. For a given measurement situation (fixed value of ti/ At), rate measurements will be made at different points along the appropriate line as the rate constant changes. Optimum measurement times can therefore be selected along these lines. Suggested measurement ranges are listed in Table 2.1. Deviations from the assumption of equal variance in the rate measurements will naturally change the error surface presented here. 16 Agreement between the theoretical model and experimental results, at least for the methods employed here, is demonstrated below. EXPERIMENTAL SECTION Two chemical systems were employed for the reaction-rate studies carried out here. The first is a metal/complex exchange reaction which was used by Pausch and Margerum [28] for the determination of alkaline-earth metals. we employed a slightly modified version of their method for the determination of calcium. The second reaction used was the formation of 12—molybdophosphate which is used for the determination of phosphate [29,30]. A. The Exchange Reaction Calcium stock solutions were prepared from reagent grade CaCOa (J.T. Baker 00., Phillipsburg, NJ.) which had been dried overnight. Solutions of the complexing agent, 1,2-diamino-cyclohexane-N,N,N’,N’- tetraacetate (CyDTA), were prepared from its acid monohydrate (98% purity, Aldrich Chemical 00., Milwaukee, Wis.) without further purification. The stock solution of Pb(II) was prepared from reagent grade Pb(NOa): (Matheson, Coleman, and Bell, Norwood, Ohio). From these, calcium and lead working solutions were prepared. The calcium working solution for the reaction consisted Of the calcium solution to be analyzed and CyDTA added in 10% excess of the highest calcium concentration. The lead reagent consisted of Pb(II) at a 17 concentration three times greater than the highest calcium concentration. To control ionic strength, sodium acetate (Fisher Scientific Co., Fair Lawn, N.J.) was added to both solutions to a final concentration of 0.5 M.’ This, with the addition of small amounts of acetic acid, also served as the pH buffer (pH 5.8, 6.0, and 6.2). The final concentrations used were [Cafl] = 0.4, 0.8, 1.0, 1.2, 1.6 and 2.0 uM, [CyDTA] = 2.2 MA, and [Pb‘t] = 6.0 ”M. The reactions were carried out in a thermostatted reaction cell at room temperature (23 °C) and were monitored spectrophotometrically by the absorbance of the Pb-CyDTA complex at a wavelength of 260 nm. B. The Molybdophosphate Reaction Phosphate solutions in the concentration range of 1 ug/mL to 10 ug/mL P were prepared from oven—dried, reagent grade KH2P04 (Mallinckrodt Inc., St. Louis, Missouri). The molybdate reagent consisted of 0.075 M MO(VI), prepared from NarMoOvZHzO (J.T. Baker CO.), in 0.6 M nitric acid. Reactions were carried out at 21.5 °C, 26.0 °C, and 31.2 °C in a thermostatted reaction cell. The formation of 12—molybdophosphate was monitored spectrophotometrically at 400 nm. C. Instrumentation All kinetics data were collected on a computer-controlled stopped- flow system designed in our laboratory [4,31-33]. The spectro- photometric detection system consisted of a light source (Heath EU-701- 50) with tungsten and deuterium lamps, a grating monochromator (Heath EU-700), a quartz observation cell with an optical path of 1.87 $0.02 cm, 18 and a photomultiplier tube (RCA 1P28A). The monochromator (3 nm bandpass) was interfaced to the reaction cell with a quartz fiber optic light guide, and transmitted radiation passed through a quartz rod to the photomultiplier tube. For each sample analyzed, 400 data points were collected at sampling intervals of 10 to 50 ms. Each point in the final set Of data was the result of averaging eight stopped-flow "pushes". The data acquired in this manner were stored temporarily in a microcomputer [34], which was also used to control the stopped-flow instrument. Control and data acquisition routines for the microcomputer were written in the FORTH language. For analysis and mass storage, the data were sent via a direct serial link to a Digital Equipment Corporation L81 11/23 minicomputer. D. Data Analysis Data analysis was performed on the minicomputer with programs written in FORTRAN-77. Numerical differentiation Of the data was carried out by a digital version of the up/d'own integration method [15,16], which is equivalent to the modified Savitzky-Golay derivative filter described by Nevius and Pardue [35]. The number of points used in this filter ranged from 31 to 51, depending on the noise level and sampling frequency for various sets of data. Nonlinear curve fitting for the estimation of kinetic parameters was performed with a program written in our laboratory which uses the modified simplex optimization algorithm of Nelder and Mead [36]. 19 RESULTS AND DISCUSSION Results are presented here to demonstrate that the two-rate method effectively reduces errors arising from between-run variations in the first-order rate constant. Systems are examined in which the rate constant is changed through variations in pH and temperature. Experimental agreement with the error model developed earlier is shown. A. pH Effects As pointed out by Carr [10], small changes in pH can cause large changes in the rate constant for reactions which involve acids or bases. For example, a pseudo-first-order reaction with a first-order dependence on the hydrogen ion concentration will have a relative uncertainty in the rate constant given by, n/k as 2.303.». (2.7) Because of the dramatic effect of pH variations on the rate constant, we decided to use this as a basis for testing the two-rate method. The exchange reaction for the determination of calcium was chosen for this study. The stoichiometry of the analytical reaction is, CaCyDTA" + Pb”(excess) --> PbCyDTA" + Ca” (2.8) Where CyDTA is the complexing agent, 1,2-diaminocyclohexane-N,N,N’,N’- tetraacetate. The rate determining step in the reaction [27] involves 20 hydrogen ion exchange and gives rise to a first-order dependence on [H’]. The two-rate method was used to analyze five standard calcium solutions, ranging in concentration from 0.4 ”M to 2 uM, at a buffered pH of 6.0. The simplified two-rate parameter, TRP(0,1.3T), was used to generate a calibration curve from these standards. The data exhibited good linearity with a slope of 0.1374 absorbance units (AU)/uM and an intercept of 0.0395 AU (this significantly non-zero intercept was found to be due to calcium impurities in the sodium acetate used to maintain a constant ionic strength). The least squares fit gave a correlation coefficient of I‘xy = 0.9980, a standard error of the estimate (or residual root-mean-square) of so = 0.0064 AU, and a relative standard error of the estimate (Salli-yum” of RS. 2 0.039. For comparison, a calibration curve was also generated using the initial rate as the response parameter. The statistics of the resulting fit were: slope = 0.04236 AU/(s uM), intercept = 0.00815 AU/s, I‘xy = 0.9964, Be = 0.0026 AU/s, and RS. = 0.051. Even though the same initial rate measurements were used for both methods, comparison of the correlation coefficients and relative standard errors of the estimate for the two methods indicates a somewhat better fit for the two-rate method. To test the sensitivity of the two-rate method to between-run changes in the rate constant, "unknowns" of 1 UM Ca2+ were analyzed in buffered solutions of pH 6.0 (same as for standards), 5.8, and 6.2. The pH of the Pb(II) reagent was also adjusted accordingly. These small pH changes altered the rate constant significantly, but had no really noticeable effect on the calcium determination, as shown in Table 2.2. The deviations are within the margins of error expected from the 21 Table 2.2. Effect of pH variations on the initial rate method (IRM) and two-rate method (TRM) . [Ca3’] Measured (HM) X Error [08”] Taken (uM) pHa k (s'1)b IRM TRM IRM TRM 1.00 6.0 0.305 1.118 1.044 12% 4.4% 1.00 5.6 0.476 1.681 0.962 68% -3.3% 1.00 6.2 0.191 0.602 0.988 -40% -1.2% 8. Standards were run at pH 6.0. b. Rate constants were determined from a non-linear least squares fit of the data. 22 calibration curve. In contrast, the results of the initial rate method, also shown in Table 2.2, exhibit gross errors when subjected to the same pH changes. This example clearly demonstrates that the two-rate method is a viable technique for minimizing errors arising from between- run variations in the rate constant. B. Thermal Effects The effects of temperature on the reaction rate constant have been thoroughly studied. With the exception of special cases, such as enzyme reactions, the magnitude of this dependence can be directly related to the activation energy, as given by Carr [10]. o' E 7f- .. iii-2L «T (2.9) The range of activation energies normally encountered with reaction-rate methods can give rise to some significant temperature effects for many reactions. Normally, reagents and reaction cells are carefully thermostatted to minimize this problem. In the experiments conducted here, the temperature was intentionally varied in order to effect a change in the rate constant so that the performance of the two-rate method could be evaluated. Calibration curves generated from five standard solutions at three different temperatures, 21.5 °C, 26.0 °C, and 31.2 °C, were compared. The reaction used for this study was the molybdophosphate reaction for the determination of phosphate [29,30]. Although the 23 mechanism for this reaction is quite complex [37,38], the analytical reaction can be represented in simplified form by, P09" + xMo(VI) --> 12-MPA (yellow) (2.10) where lZ-MPA is the heteropolyanion, 12-molybdophosphate. This reaction is pseudo-first-order in acidic solution with excess molybdate and can be followed spectrophotometrically on the stopped-flow time scale. The rate constants measured at the three temperatures given above were, in order of increasing temperature, 0.34 s", 0.67 s", and 1.12 r1, values in agreement with what was expected from the previously measured activation energy [37]. Three different methods were used to generate calibration curves from the data at the three different temperatures. The methods used were the initial rate method, the fixed-time method, and the two-rate method. For variety, the generalized two-rate method was used with measurement times of t1=0.51 and ta=21'. The slopes Of the calibration curves Obtained by each of these methods are reported in Table 2.3. Other least squares statistics have not been reported since all of the methods gave intercepts which were effectively zero and exhibited good linearity, with correlation coefficients greater than 0.999 and relative standard errors of the estimate not exceeding 4%. The results reported in Table 2.3 show that the slope of the calibration curve obtained by the initial rate method has the largest dependence on temperature (rate constant), while the two-rate method exhibits the smallest dependence. This is in accordance with theoretical expectations. An unexpected trend in the data is the positive 24 Table 2.3. Slope of calibration curves for three reaction-rate methods at three temperatures. Slope of Calibration Curve Method T = 21.5 °C T = 26.0 ’0 T = 31.2 °C Initial Rate Method' 0.0150 0.0291 0.0587 Fixed Time Method” 0.0264 0.0336 0.0391 Two-Rate Method" 0.0402 0.0411 0.0451 a. Slope in units of AU mL s‘1 ug‘l. b. Slope in units of AU mL pg”. Used t1 0, t2 = 21 (1 calculated from data at 26.0 °C). c. Slope in units of AU mL ug’l. Used t1 0.51, t: = 21. 25 correlation between the values obtained by the two-rate method and the temperature. This is especially noticeable with the larger than expected deviation at the highest temperature. This deviation from theory is attributed to the method of rate measurement rather than deficiencies in the two-rate method itself. The cutoff frequency of the digital derivative filter we employed began to cause significant distortions in the relatively fast rising high temperature data and introduced the errors. These effects are not readily apparent until It becomes large relative to the cutoff frequency, so the two lower temperature values are in better agreement. C. Propagation of Random Errors To investigate the validity of the model for the propagation of random errors experimentally, ten replicate sets of data were Obtained with the Pb-Ca/CyDTA system. The calcium concentration for all of these runs was 1 uM and all measurements were made at a buffered pH of 6.0 in a thermostatted cell. These conditions minimized variations in the rate constant, which exhibited only a 1% relative standard deviation between runs as determined by non-linear fitting of the data. Thus, any errors in the rate measurements are expected to be primarily due to sources other than variations in the rate constant. The relative standard deviation in the two-rate parameter was calculated from these replicate samples as a function of time. The results, expressed as a multiple of the relative standard deviation in the rate measurement at 1 (as they were in Figure 2.1), are shown in Figure 2.2. The generalized form of the two-rate parameter, with t: = 0.51, was used so that the most general form of the error model could be assessed. The plot shows 26 0| 0 N O l _I O L .° 0 l 1 I . l I U 1 l r I ' I l 1 1.0 2.0 3.0 4.0 5.0 At (in units of 1) RSD In TRP(O.5.t2) (In units of RSD(R.,.)) .0 0 Figure 2.2. Comparison Of theoretical and experimental random errors for the two-rate method as a function of measurement time. The dashed line represents the theoretical model. 27 the RSD in TRP(0.51,ta) vs. At (where t2 = 0.51+At). The experimental data are represented by the solid line, while the dashed line shows the theoretical behavior. Even with the statistically small sample size used for this study, the agreement between theory and experiment is quite good, and this supports the validity of the model. The general trends in the experimental results are as expected and there is good quantitative agreement in the critical region of 0.51-21. Deviations from the theoretical curve can be attributed to the breakdown of assumptions leading to the model, especially that the absolute errors in the rate measurements are independent of time. CONCLUSIONS The results presented here have demonstrated that the two-rate method can be used to minimize errors arising from between-run variations in the rate constant for systems following first-order or pseudo-first-order kinetics. The method assumes that the rate constant does not vary within a run and also that the equilibrium position of the reaction is not affected by changes in the rate constant. The model developed for the propagation of random errors in rate measurements shows that there are limitations in terms of the measurement times chosen. The error model was found to agree well with experimental results, at least for the system used in this study. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 28 LIST OF REFERENCES H.V. Malmstadt, C.J. Delaney and E.A. Cordos, CRC. Crit. Rev. Anal. Chem., 2, 559-619 (1972). H.L Pardue, in "Advances in Analytical Chemistry and Instrumentation", C.N. Reilley and F.W. McLafferty, Eds.; Wiley, New S.R. Crouch, in "Computers in Chemistry and Instrumentation", H.D. Mattson, H.B. Mark, Jr. and H.C. MacDonald, Eds.; Dekker, New S.R. Crouch, F.J. Holler, P.K. Notz and P.M. Beckwith, Appl. Spec. Rev., 13, 165-259 (1977). H.A. Mottola and H.B. Mark, Jr., Anal. Chem., 56, 96R-112R (1984). H.B. Mark, Jr. and G.A. Rechnitz, "Kinetics in Analytical Chemistry"; Wiley-Interscience, New York, 1968. K.B. Yatsimirskii, "Kinetic Methods of Analysis"; Pergammon, Oxford, 1966. ' First International Symposium on Kinetics in Analytical Chemistry; Cordoba, Spain, September 27-30, 1983. Second International Symposium on Kinetics in Analytical Chemistry; Preveza, Greece, September 9-12, 1986. P.W. Carr, Anal. Chem., 50, 1602-1607 (1978). W.J. Blaedel, and C. Olson, Anal. Chem., 36, 343-347 (1963). G.E. James, and H.L. Pardue, Anal. Chem., 40, 796-802 (1968). S.R. Crouch, Anal. Chem., 41, 880-883 (1969). (1970). E.M. Cordos, S.R. Crouch, and H.V. Malmstadt, Anal. Chem., 40, 1812-1818 (1968). J.D. Ingle, Jr. and Crouch, S.R., Anal. Chem., 42, 1055-1060 (1972). R.H. Calicott and P.W. Carr, Anal. Chem., 46, 1840-1842 (1974). E.S. Iracki and H.V Malmstadt, Anal. Chem., 45, 1766-1770 (1973). 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 29 J.G. Atwood and J.L. DiCesare, Clin. Chem. (Winston-Salem, N.C.), 49, 1263-1269 (1973). J.B. Landis, M. Rebec and H.L Pardue, Anal. Chem., 49, 785-788 (1977). F.J. Holler, R.H. Calhoun and S.F. McClanahan, Anal. Chem., 54, 755-761 (1982). J.E. Davis and B. Renoe, Anal. Chem., 51, 526-528 (1979). G.E. Mieling and H.L. Pardue, Anal. Chem., 50, 1611-1618 (1978). P.D. Wentzell and S.R. Crouch, Anal. Chem., 58, 2851-2855 (1986). P.D. Wentzell and S.R. Crouch, Anal. Chem., 58, 2855-2858 (1986). J.D. Ingle, Jr. and S.R. Crouch, Anal. Chem., 43, 697-701 (1971). J.B. Pausch and D.W. Margerum, Anal. Chem., 41, 226-232 (1969). W. Rieman and J. Beukenkamp, in "Treatise on Analytical Chemistry", LM. Kolthoff and P.J. Elving, Eds.; John Wiley and Sons, New York, 1961; Part II, Vol. 5, pp. 348-351. A.C. Javier, S.R. Crouch, and H.V. Malmstadt, Anal. Chem., 41, 239-243 (1969). P.M. Beckwith and S.R. Crouch, Anal. Chem., 44, 221-227 (1972). R. Balciunas, Ph.D. Dissertation, Michigan State University, East Lansing, MI, 1981. R. Putt, M.S. Thesis, Michigan State University, East Lansing, MI, 1983. B. Newcome and C.G. Enke, Rev. Sci. Instrum., 55, 2017-2022 (1984). T.A. Nevius and H.L. Pardue, Anal. Chem., 56, 2251-2253 (1984). J.A. Nelder and R. Mead, Computer J., 7, 308-313 (1965). P.M. Beckwith, A. Scheeline and S.R. Crouch, Anal. Chem., 47, 1930-1936 (1975). 0.0. Kircher and S.R. Crouch, Anal. Chem., 55, 242-248 (1983). CHAPTER III Comparison of Reaction-Rate Methods of Analysis for Systems Following First—Order Kinetics "If it can’t be expressed in figures, it is not science; it is opinion. " - Robert Heinlein "Time Enough for Love" The popularity of reaction-rate methods in chemical analysis is indicated by the wide variety of methodologies, which have been developed over the years [1-6]. Among the most widely used rate methods are the fixed-time [1,7,8], the variable-time [1,9-11], and the derivative [1,12,13] methods. These traditional reaction-rate methods are often susceptible to variations in experimental parameters which affect the rate constant, such as pH, temperature, ionic strength and reagent concentration. Carr [7] has shown how small variations in the first two variables can adversely affect the precision of the results. In recent years, several workers have attempted to alleviate the problem of between-run variations in the rate constant for systems following first- or pseudo-first-order kinetics. Atwood and DiCesare [14] first noted that the effect of rate constant variations could be minimized by appropriate adjustment of enzyme activity in substrate determinations by kinetic methods. They suggested that the enzyme activity should be adjusted so that the reciprocal of the pseudo-first- order rate constant is equal to the time of measurement. Pardue and 30 31 co-workers [15] extended this observation, noting that for any reaction following first-order kinetics, the optimum rate measurement time for minimizing the effect of between-run variations in the rate constant is t = 1/k = 1. This observation was investigated further by Holler et a1. [16], who demonstrated experimentally that improved results could be obtained with the method. For convenience, this method of measuring rates at t = 1 is referred to as the optimized derivative method throughout this chapter. Davis and Renee [17] described an optimized fixed-time approach which reduces the influence of between-run rate constant variations for the traditional fixed-time method. This method Optimizes measurement times for minimal rate constant dependence. Mieling and Pardue [18] have proposed a multiple-linear-regression procedure which evaluates kinetic parameters to compensate for changes which occur in the reaction curve. A method developed by Cornell [19] in 1962 for fitting exponentials through partial sums is also applicable to the problem of between-run variations in the rate constant. While the Cornell method has been extended to other kinetic systems by Kelter and Carr [20-22], it has not appeared extensively in the analytical chemistry literature and in this sense its application to reaction-rate methods of analysis is fairly recent. Wentzell and Crouch [23] have introduced a new method, the two-rate method, which attempts to eliminate the dependence on between-run variations in the rate constant by using rate measurements made at two times during the course of the reaction. In this chapter, a critical comparison of several reaction-rate methods applied to systems following first- or pseudo-first-order kinetics is made [24]. Both traditional methods (fixed-time, variable- 32 time, derivative, and initial rate methods) and the more recently developed techniques (optimized derivative, Optimized fixed-time, Cornell, and two-rate methods) are considered. Some of the more specialized methods, such as the kinetic difference [25,26] and signal-stat [27,28] methods, are not included in this comparison. The regression-kinetic method Of Mieling and Pardue [18] is also excluded, since .its level of sophistication is beyond the scope of this treatment. Otherwise, we have attempted to be as general as possible, evaluating each method in terms of susceptibility to rate constant variations, precision in the absence of those variations, and several other factors. The experimental results presented are limited to spectrophotometric data, but many of the conclusions drawn are generally applicable. Throughout this discussion, reaction-rate methods are classified in two ways for ease of reference. First, adopting the terminology of Pardue et a1. [29], those methods requiring an instantaneous measurement of reaction rate (derivative, optimized derivative, initial rate, and two-rate methods) are defined as differential methods, while the others (fixed-time, variable-time, optimized fixed-time, and Cornell methods) are referred to as integral methods. In addition, "traditional" methods (fixed-time, variable-time, initial rate, and derivative) are distinguished from "optimized" methods (Optimized derivative, optimized fixed-time, Cornell and two-rate), which attempt to minimize susceptibility to rate constant variations. 33 EXPERIMENTAL SECTION Two chemical systems were used to evaluate the reaction-rate methods discussed in this paper. Phosphate determinations were carried out by means of the 12-molybdophosphate reaction [30,31]. Calcium determinations used the metal/complex exchange reaction described by Pausch and Margerum [32]. All reactions were carried out on an automated stopped-flow system [33,34] and data were processed on a DEC LSI 11/23 computer. Further experimental details appear in the preceding chapter and elsewhere [23]. RESULTS AND DISCUSSION Errors in a reaction-rate method can be divided into two broad categories: (1) errors arising from between-run variations in the rate constant, and (2) errors arising from other sources. Each of these is considered for the methods compared. In addition, other factors which may influence the selection of a rate method are discussed. A. Susceptibility to Rate Constant Variations A theoretical comparison Of the effect of between-run rate constant variations on seven different reaction-rate methods was carried out to evaluate the extent of the dependence in each case. The methods examined were the fixed-time, variable-time, initial-rate, optimized derivative, optimized fixed-time, Cornell, and two-rate methods. 34 Figure 3.1 compares the theoretical percent error in the concentration determined by each method as a function Of the deviation Of the rate constant from its nominal value, ko. Equations used to generate the plot are listed in Table 3.1 and their derivation is straightforward. In generating Figure 3.1, measurement times of t1=0 and t2=1/ko were assumed for the fixed-time method. Also, for the variable-time method, the reaction was assumed to be 2% complete at the time of the measurement. Finally, for the Optimized fixed-time method, measurement times of 0.5/k6 and 1.76/ko, values given by Davis and Renoe [17], were chosen. Examination of Figure 3.1 leads to some interesting conclusions. Of the six methods, the variable-time method and the initial rate method exhibit the largest dependence on the rate constant. For these methods, the reaction curve is very nearly linear in the measurement region and, therefore, the measured parameters are directly proportional to the first-order rate constant. The fixed-time method shows a somewhat smaller dependence on 11, but this is highly dependent on the selection of measurement time. Short measurement intervals increase the dependence on the rate constant, approaching that of the initial rate method at the limit. As the interval becomes longer, there is an approach to equilibrium methods, which have virtually no dependence on the rate constant. The reason for improved results with the optimized derivative and optimized fixed-time methods is easily seen from Figure 3.1. Both of these methods exhibit an asymptotic approach to zero error at the point of zero deviation in the rate constant. Because of the flatness of the curves in this region, the effect of variations in the rate constant is 35 200. . l r I ' C 1 - Fixed-Time Method 0 2 - Variable-Time and Initial Rote MethodI :3 " 3 - Optimized Derivative Method " o 4 - Optimized Fixed-Time Method :3 5 - Two-Rate and Cornell Methods c — - o 100. "I O C O o .1- t C L 0 --—k o—J O L L. m -)- N l l 1 L L -100e ' I ' T I -1 .0 0.0 1.0 2.0 Ak/ko Figure 3.1. Theoretical errors expected in the concentration determined by various reaction-rate methods as a function of the relative deviation Of the rate constant from its nominal value. 36 Table 3.1 Equations for relative errors arising from rate constant variations for selected reaction-rate methods“. Method Relative Error in Concentration Initial Rate Alt/kc Variable-Time -0.01003 -I- 0.98997(Ak/ko) as Ak/ko (2% extent of rxn.) Fixed-Time (t = 1/ k0) Optimized Fixed-Time (t1=0.51, t:=1.761) Optimized Derivative Cornell, Two-Rate "Nf/ko 0.58198(1-e ) -0.5Ak/lro -1.76Alr/Iro_ 0.606(e -1)-0.172(e 0.434 1) Air/Ito (1 + Ak/ko)e- -1 No dependence on k (in theory) ako is the nominal rate constant, k is the actual rate constant during the run, Ak = k - ha, 1 = 1/ko. 37 minimized, although it is not zero. The similarity of the two curves is worth noting since it indicates an approximately equivalent dependence on k for the two methods. The optimized fixed-time method is expected to improve in this respect as a wider measurement interval is chosen, but again, this is due to an approach to equilibrium methods. Another point worth noting with regard to these two methods is that negative errors in the concentration result regardless of the direction of the variation Of the rate constant. This is an undesirable feature since it introduces bias. The Cornell and two-rate methods show, in theory, no errors attributable to between-run variations in the rate constant. The two- rate method compensates for these variations by making rate measurements at two times during the reaction and is limited in practice by random errors in the rate measurements [23]. The Cornell method adjusts to between-run changes in the rate constant by employing a relatively simple algorithm to fit the exponential data. The practical limitations Of this method have not been extensively investigated. Table 3.2 illustrates experimentally the performance of several reaction-rate methOds under conditions of gross changes in the rate constant. The results were Obtained using the molybdophosphate reaction for the determination of phosphate. The rate constant was" varied by changing the temperature. Calibration curves were generated for each method at a temperature of 26.0 °C. "Unknowns" of 4 ug/mL and 8 ug/mL P were then analyzed at this temperature and also at 21.5 °C and 31.2 °C. The errors in the measured phosphate concentrations are reported in Table 3.2. As anticipated, the variable- time and initial rate methods exhibit the greatest errors when the 38 Table 3.2. Comparison of results Obtained by various reaction-rate methods under conditions of large changes in temperature. Based on calibration curves generated at 26.0 °C. % Error in Phosphate Measured Phosphate Method Taken (ug/mL P) T:21.5°C T:26.0°C T=31.ZRC Variable-Time 4.00 -43% -4.1% 59% 8.00 -45% -7.8% 52% Initial Rate 4.00 -49% 2.7% 78% 8.00 -49% 0.5% 88% Fixed-Time ( 11) 4.00 -40% 0.4% 42% 8.00 -33% 0.4% 39% Fixed-Time (21) 4.00 -24% -0.5% 16% 8.00 -22% -0.2% 18% Optimized Fixed-Time 4.00 -14% -3.0% -6.2% (0.5t,1.761) 8.00 -16% 1.5% -12% Optimized Derivative 4.00 -14% 1.1% -12% (0.751,1.51) 8.00 -0.6% -l.7% 7.2% Cornell Method 4.00 -0.1% 0.6% 6.0% (Range=2.7 81 ) 8.00 -0.4% 0.16% 5.1% 39 temperature is changed. As predicted by Ingle and Crouch [35], the variable-time method also gives poor results at the temperature at which the calibration curve was generated, which demonstrates that this method is best applied to systems following zero-order kinetics. The fixed-time results (at t = 1) are somewhat better than the variable-time results, but the errors are still large. The errors resulting from the application of the optimized derivative (t = 1) and the Optimized fixed- time (t1 = 0.51,t2 = 1.761) approaches are, as predicted, comparable and more reasonable. ' Also as predicted, there is a consistently negative bias with these two methods. Finally, the two-rate method (t1 = 0.751,ta = 1.51) and the Cornell method (range Of data used = 2.781) show the least dependence on changes in the rate constant, although the errors at the high temperature seem inordinately large. These high-temperature deviations have been discussed in the preceding chapter and elsewhere [23] for the two-rate method. For the Cornell method, the errors are believed to be due to a limitation of the algorithm which becomes more important as the rate constant increases. B. Precision with Invariant Rate Constant We now consider the case where between-run variations in the rate constant are negligible and errors arising from other sources, such as detector noise or photometric source fluctuations, limit precision [29,35]. To examine precision in the absence Of rate constant variations, ten replicate sets Of data were Obtained for the calcium determination. Rate constant fluctuations were minimized by maintaining careful pH and temperature control. The data Obtained were analyzed 40 by several methods, and the relative standard deviation (RSD) calculated in each case is reported in Table 3.3. Among the differential methods, four techniques are compared: the derivative method (at t = 0.51), the Optimized derivative method, and the two-rate method with two sets of measurement times. The initial rate method is not included since the method of rate calculation is significantly different and would not allow a reliable comparison. Rates for the methods compared were computed by employing a digital version of the up/down integration method [36,37], which is equivalent to a modified Savitzky-Golay derivative filter [38,39]. For our system, we found that the absolute errors in the rate measurements were largely independent Of the time Of measurement. This is a consequence of the method of rate calculation and the fact that absolute errors in the absorbance measurements remained effectively constant over the small absorbance range used. In Table 3.3, the two-rate method shows similar or improved precision over the optimized derivative method, as predicted elsewhere [23]. A single rate measurement at an earlier time shows improved precision over both of these due to a larger rate with comparable errors in the rate measurements. However, this advantage is overshadowed by a greater susceptibility to rate constant variations. Among the integral methods considered, Ingle and Crouch [35] have already demonstrated the inferior precision of the variable-time method when applied to first-order systems. In the current study, the errors would be compounded by a relatively long sampling interval (poor time resolution), so the results for this method are not included in Table 3.3. The remaining integral methods exhibit a similar compromise 41 Table 3.3. Precision of various reaction-rate methods under conditions Of an invariant rate constant. Method Parameters RSD Derivative $0.51 0.57% Optimized Derivative t=1 1.1% Two-Rate ’ t1=0.51 , ta=1 1 1.1% Two-Rate t1=O.51 , tz=1.7 61 0.86% Fixed-Time t1=0, t2=11 1.2% Fixed-Time t1=0, tz=21 0.83% Optimized Fixed-Time t1=0.81, ta=1.231 5.0% Optimized Fixed-Time t1=0.51, ta=1.761 1.7% OFT (smoothed) tit-0.51, tar-1.761 0.41% Cornell Range=2.251, 150 points 0.34% Cornell Range=4.51, 300 points 0.38% Cornell Range=61 , 399 points 0.67% 42 between precision and rate constant dependence as was seen for the differential methods. Pardue et a1. [29] have shown that for small absorbance changes, as is the case here, relative errors can be expected to be inversely proportional to the magnitude of the absorbance change. This is supported by the results in Table 3.3. Both the fixed-time and Optimized fixed-time methods show improved precision when the time interval is increased, approaching the precision of equilibrium methods in the limit. Also, since the fixed-time method (with t1 = 0) utilizes larger absorbance changes than the Optimized fixed-time method for equivalent time intervals, the latter method generally exhibits poorer precision in the absence of rate constant variations. The Cornell method exhibits the best precision of the methods examined, primarily because all Of the data are used. The precision Obtained will be a function of the number of points and the range of data employed in determining the partial sums, as shown in the table. The exact nature Of this dependence warrants further investigation. A direct comparison of the precision Observed for differential and integral methods would be unfair since the method Of calculating the parameters is not the same. Differential methods require differentiation of the analytical signal to compute the instantaneous rate. Since this process amplifies high frequency noise [6,39], some smoothing Of the data is normally required. In this work, smoothing was carried out concurrently with differentiation by a modified Savitzky-Golay derivative filter [38,39]. To Obtain a fair comparison of the precision between differential and integral methods, it was necessary to apply a comparable smoothing filter to the data for the latter. The last entry for the Optimized fixed-time method in Table 3.3 gives the relative 43 standard deviation the results when such a filter was applied. The filter used was a simple Savitzky-Golay [40] linear smooth. The result indicates that, when compared on common ground, the integral methods are more precise than the differential methods. This observation is expected to be generally valid due to the lower precision of instantaneous rate measurements. C. Additional Factors Several other factors may be of importance when considering which rate method to use in a given application. For example, the compatibility Of the rate computation technique to continuous flow methods of analysis may be important if determinations are to be carried out in flowing streams. Of the methods discussed here, only the fixed- time method and the optimized fixed-time method are directly compatible with flow injection analysis (FIA) [41] and air-segmented continuous flow analysis [42]. The other methods can be implemented in flow systems only if the flow is stopped. Another important factor is the instrumental and computational complexity of the rate methods. All of the methods discussed have the same fundamental components in the instrumental design, as shown in Figure 3.2, so the complexity involved in their implemention is largely determined by the amount of pre- and/or post-processing Of the data required. The integral methods are all implemented with relative ease. With the exception of the Cornell method, each requires a simple calculation of the difference between two measurements. Differential methods, on the .Other hand, require differentiation of the analytical signal. Many effective electronic differentiation methods have been 44 apogee. odouIsoSoaos sea 2:39? we eases—30 Eases—«ZEEH .3925; no .N.« 953% :oo 65:: m m :3 mm coSoEomno cOZEOm Eanm even so oc_eeeoo..alaeoa 00_>00 .33.» 5o mcmflaooom us_eeeoo..ale..a _ LOLOOLOQ acomoom 45 developed [36,37 ,43-45], but all require careful implementation to avoid degradation of precision due to noise amplification. Numerical methods [38-40] are another option, but again careful selection Of the digital filter is required to ensure sufficient noise immunity and minimal signal distortion. Thus, the necessity for a direct rate calculation is a complicating aspect of differential methods. CONCLUSIONS Each of the reaction-rate methods examined exhibits its own particular advantages and disadvantages when applied to systems following first-order kinetics. An exception to this is the variable-time method, which is more applicable to zero-order reactions. As expected, the optimized methods are less sensitive to between-run fluctuations in the rate constant than traditional methods, but generally show poorer precision in the absence Of these fluctuations. In this study, the Cornell method exhibited the best overall performance. Although it is a multi-point method which requires measurements near t=0, it is insensitive to between-run variations in the rate constant, gives precise results in the absence of those variations, and is relatively simple to implement. The analytical limitations of this method require further investigation, however. 3. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 46 LIST OF REFERENCES H.V. Malmstadt, C.J. Delaney and E.A. Cordos, CRC Crit. Rev. Anal. Chem., 2, 559-619 (1972). H.L. Pardue, Clin. Chem. (Winston-Salem, N.C.), 23, 2189-2201 (1977). P.W. Carr and L.D. Bowers, "Immobilized Enzymes in Analytical and Clinical Chemistry"; Wiley, New York, 1980; Chapter 3. H.B. Mark, Jr. and G.A. Rechnitz, "Kinetics in Analytical Chemistry"; Wiley, New York, 1968. H.L. Pardue, in "Advances in Analytical Chemistry and Instrumentation", C.N. Reilley and F.W. McLafferty, Eds.; Wiley, New S.R. Crouch, in "Computers in Chemistry and Instrumentation", H.D. Mattson, H.B. Mark, Jr. and H.C. MacDonald, Eds.; Dekker, New W.J. Blaedel and C. Olson, Anal. Chem., 36, 343-347 (1963). G.E. James and H.L. Pardue, Anal. Chem., 40, 796-802 (1968). R.A. Parker, H.L. Pardue and 3.6. Willis, Anal. Chem., 42, 56-61 (1970). H.L. Pardue, Anal. Chem., 36, 1110-1112 (1964). J.G. Atwood and J.L. DiCesare, Clin. Chem. (Winston-Salem, N.C.), 19, 1263-1269 (1973). J.B. Landis, M. Rebec and H.L. Pardue, Anal. Chem., 49, 785-788 (1977). F.J. Holler, R.K. Calhoun and S.F. McClanahan, Anal. Chem., 54, 755-761 (1982). J.E. Davis and B. Renoe, Anal. Chem., 51, 526-528 (1979). G.E. Mieling and H.L. Pardue, Anal. Chem., 50, 1611-1618, (1978). R.G. Cornell, Biometrics, 18, 104-113 (1962). P.B. Kelter and J.D. Carr, Anal. Chem., 51, 1825-1828 (1979). 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 37. 38. 39. 40. 41. 42. 43. 47 P.B. Kelter and J.D. Carr, Anal. Chem., 51, 1828-1834 (1979). P.B. Kelter and J.D. Carr, Anal. Chem., 52, 1552 (1980). P.D. Wentzell and S.R. Crouch, Anal. Chem., 58, 2851-2855 (1986). P.D. Wentzell and S.R. Crouch, Anal. Chem., 58, 2855-2858 (1986). H. Weisz and H.Ludwig, Anal. Chim. Acta, 55, 303-313 (1971). H. Weisz and S. Pantel, Anal. Chim. Acta, 68, 311-316 (1974). H.V. Malmstadt and E.H. Piepmeier, Anal. Chem., 37, 34-44 (1965). H. Weisz and K. Rothmaier, Anal. Chim. Acta, 75, 119-126 (1975). H.L. Pardue, T.E. Hewitt and M.J. Milano, Clin. Chem. (Winston- W. Rieman and J. Beukenkamp, in "Treatise on Analytical Chemistry", LM. Kolthoff and P.J. Elving, Eds.; Wiley, New York, A.C. Javier, S.R. Crouch and H.V. Malmstadt, Anal. Chem., 41, 239-243 (1969). J.B. Pausch and D.W.Margerum, Anal. Chem., 41, 226-232 (1969). P.M. Beckwith and S.R. Crouch, Anal. Chem., 44, 221-227 (1972). S.R. Crouch, F.J. Holler, P.K. Notz and P.M. Beckwith, Appl. Spec. Rev., 13, 165-259 (1977). J.D. Ingle, Jr. and S.R. Crouch, Anal. Chem., 43, 697-701 (1971). E.M. Cordos, S.R. Crouch and H.V. Malmstadt, Anal. Chem., 40, 1812-1818 (1968). J.D. Ingle, Jr. and S.R. Crouch, Anal. Chem., 42, 1055-1060 (1972). T.A. Nevius and H.L. Pardue, Anal. Chem., 56, 2249-2251 (1984). R.W. Hamming, "Digital Filters"; Prentice-Hall, Englewood Cliffs, NJ, 1983. A. Savitzky and M.J.E. Golay, Anal. Chem., 36, 1627-1639 (1964). J. Ruzicka and E.H. Hansen, "Flow Injection Analysis"; Wiley, New York, 1981. W.B. Furman, "Continuous Flow Analysis"; Dekker, New York, 1976. H.V. Malmstadt and S.R. Crouch, J. Chem. Ed., 43, 340-353 (1966). 48 44. R.H. Calicott and P.W. Carr, Anal. Chem., 46, 1840-1842 (1974). 45. E.S. Iracki and H.V. Malmstadt, Anal. Chem., 45, 1766-1770 (1973). CHAPTER IV The Kalman Filter for Reaction-Rate Methods Of Analysis "It is the discipline of science that enables ordinary people, whether chemists or biologists, to do the ordinary things, which, when assembled, reveal the extraordinary intricacies and awesome beauties Of nature. Science permits them not only to contribute to grand enterprises, but offers them a Changing and endless frontier for exploration." -Arthur Kornberg Nobel Laureate biochemist The importance Of reaction-rate methods of analysis and their susceptibility to between-run variations in the first-order rate constant have been discussed in Chapters 11 and III. SO far, four methods for reducing the dependence of reaction-rate methods on between-run variations in the rate constant have been examined: the optimized fixed- time method [1], the optimized derivative method [2], the Cornell method [3-6], and the two-rate method [7]. In addition, the regression-kinetic method of Mailing and Pardue [8] has been mentioned but not extensively examined due to the complexities of its implementation. In this chapter, yet another approach to the problem of rate constant variations is presented in the form of the Kalman filter. The approach uses the extended Kalman filter for recursive nonlinear fitting Of experimental data and most closely resembles the regression-kinetic method. The reason for the omission of the Kalman filter method from the comparison in Chapter III is partly chronological, but also because, 49 50 like the regression-kinetic method, its implementation is radically different from traditional reaction-rate methods. In addition to minimizing the effect of rate constant variations, the Kalman filter has other applications to reaction-rate methods of analysis which are discussed in this chapter. Its application to differential reaction-rate methods is one area which has been studied with simulated data. Differential kinetic methods are those methods which take advantage of differences in the rate constants for parallel reactions to determine two or more analytes simultaneously [8-10], and have been used for many years [12-17]. Currently these methods are limited by the large errors which can arise as the ratio of the rate constants for the two reactions approaches unity. It is demonstrated in this chapter that the application of the Kalman filter can substantially reduce this constraint, permitting a wider range of rate constant ratios to be used, and provide additional information such as error estimates. The Kalman filter is a method which has been adapted from electrical engineering. It may be considered a modeling method since it requires the definition of a system model and a measurement model. This chapter begins with a theoretical discussion of the Kalman filter algorithm and is followed by a description of the application of this approach to the problems described above. Other potential applications are described and conclusions are drawn regarding the utility of the methods. 51 GENERAL CONSIDERATIONS In this section, the general algorithm for the Kalman filter is given. Particular applications are then described and the corresponding models are presented. A. The Kalman Filter In recent years the Kalman filter [18-20] has been applied to a number of (areas in analytical chemistry and is the subject of several recent review articles [21-22]. The theoretical aspects Of the discrete Kalman filter have been well treated elsewhere [19-23], but are included here for the sake Of completeness. Consider a time-variant process, the system, frOm which noisy measurements are Obtained at discrete time intervals, tIx. As with other filters, the process does not really need to be time variant, but must consist of measurements made at discrete intervals. Defined succinctly, the Kalman filter is an algorithm for recursive estimation Of state parameters from the series of noisy measurements in accordance with a linear system model. The requirement Of a linear model for our system is a prerequisite for use of the Kalman filter. The parameters of this model, which are unknown prior to filtering, are the state parameters and collectively form the state vector. For example, the state vector may consist of the consentrations of several species in a solution. It is not required to be a static entity. One way in which the Kalman filter differs from ordinary digital filters is that it is a state space method. This means that the Objective of the algorithm is to Obtain the best estimate Of the state vector and, 52 consequently, Of the true signal. To do this, we must first define a measurement model Of the form, 2k = HI: xx + v: (4.1) In this model z: is an m x 1 vector which consists of 111 measurements made. at interval 1: (i.e. for k = 1,2,...). Often m = 1 and z: is a scalar consisting of a single observation of the process. The n x 1 vector x1: is the process state vector at interval k, where n is the number Of state parameters. The 111 x 11 matrix Hy, which will be called the observation matrix, provides the connection between 2: and n under ideal (noiseless) conditions. Finally, v: is the measurement error vector (m x 1), which is assumed to be a white noise sequence. In addition to the measurement model, the Kalman filter requires that a model be defined for the propagation of the state vector between discrete measurement intervals. This model is of the form, n+1 2 Fl: x1: + wk (4.2) where F1: is the n x n state transition matrix which describes how the process state vector, xx, propagates between measurements. The n x 1 white noise vector Wk accounts for noise in the state vector (parameter estimates). Often eq 4.2 is simplified when the state vector is time invariant (F1: = I, where I is the identity matrix) and considered to be noise-free (wk 1' 0). 53 Once the above models have been defined, application of the Kalman filter is relatively straightforward. First, the Kalman gain is computed according to, KI = PI' HIT (HI: Pk’ HIT + Rk)"1 (4.3) In the notation used here, the "T" superscript indicates the transpose of a matrix and the "-" superscript indicates that this is our best estimate of the entity prior to assimilating the new measurement at interval k. In eq 4.3, K: is the n x 1 Kalman gain vector which will describe how the new measurement is to be used to update the state estimate. The n x 11 matrix Pr is the current estimate of the error covariance matrix for the state parameters, defined in the usual manner [24]. Typically, since we have no a priori knowledge of the state vector, the diagonal elements of this matrix are initially taken to be large and the off-diagonal elements are taken to be zero [25]. Since the covariance matrix is updated by the Kalman filter algorithm, only the initial value needs to be provided. The matrix RI: is the m x m covariance matrix for vu. Typically, the diagonal elements of this matrix are simply the variances of the associated measurements, 21:, and the Off-diagonal elements are zero. The next step in the Kalman filter algorithm is to assimilate the measurement at interval k to update the estimate of the state vector. The correction factor for the current state estimate is the difference between the actual measurement vector, zu, and the predicted measurement vector, Hxxr, weighted by the Kalman gain vector, KR, as shown in eq 4.4. 54 x1: = xx' + K11 (z: - 1111 x11“) (4.4) In a similar manner, the Kalman gain is used to estimate the updated covariance matrix, P11, P1: = (I - Kr H11) Pk' (I - KI: H11)? + K: RI K11" (4.5) Simpler forms of eq 4.5 are sometimes given, but are seldom used in practice because Of poor numerical stability. The final step in the algorithm is to project the estimates. Of xx and P11 ahead tO the next measurement interval using the state transition matrix, F11 . n+1 = F: xx (4.6) P191 = F11 P11 F11" + Q1: (4.7) In eq 4.7, Q1: is the covariance matrix for the white noise sequence wk introduced in eq 4.2. It is Often taken to be zero. Equations 4.3 through 4.7 constitute one form Of the Kalman filter algorithm. A number of variations of this basic algorithm are used [20,21], primarily for improved numerical stability. In the work presented in this chapter, algorithms employed the Potter-Schmidt square-root algorithm [20,26] which is known to be more stable when large dynamic ranges are encountered. The recursive procedure defined by the equations, when properly employed, will provide the Optimal estimate of the state vector, x11, according to a minimum mean-square error criterion. While Optimizing by minimum mean-square error is not, 55 strictly speaking, the same as optimizing by least squares, the two criteria provide identical results under certain conditions. Indeed, the Kalman filter is Often referred to as a "least squares filter" and has been used for linear least-squares estimation [25,27-29]. Figure 4.1 is a qualitative representation Of the factors which go into the update of the state vector after each measurement. The bottom of the figure shows the representation of a two parameter state vector in state space and its updated value. The circles surrounding the points in state space represent the errors associated with the parameter estimates (i.e. the covariance matrix) and in the general case would actually be elliptical. Note that there is a reduction in the error after the update is made. The measurement sequence is represented by the noisy signal at the top of the figure. There are essentially four factors which will determine the updated state vector. The first is the error in prediction, which is the difference between the predicted and actual measurements. A large prediction error will cause a greater adjustment to the state vector than a small one, since some Of the difference is assumed to be associated with the state parameters. The second important factor is the estimated measurement error. Very noisy measurements will result in smaller adjustments to the state vector, since this would account for a significant portion of the prediction error. A third factor is the error associated with the state vector itself, the covariance matrix. If the elements Of the covariance matrix are large, as they normally are early in the measurement sequence, the uncertainty in the state vector is great, so a large adjustment will be made. As more measurements are made, the uncertainty should become less and smaller adjustments will 56 (2) Estimate of Measurement Error Measurementu (1) Error in Prediction 0 Prediction I I I I I ,’ (4) Observation Matrix (3) Estimated Error in State Estimate (Covariance Matrix) x2 Figure 4.1. Qualitative representation of factors influencing the update Of the of the Kalman filter state vector'after a new measurement. 57 be made so that outliers do not adversely affect the estimates. In addition, the magnitude of individual elements Of the covariance matrix will determine, in part, the proportion of adjustments made to each _parameter. The final factor used in the update is the observation matrix, H, since this reflects how changes in the state parameters affect the predicted measurement. While the description in the preceding paragraph is quite qualitative, it is a useful interpretation of the principles which are an implicit part Of the mathematical representation Of the Kalman filter algorithm. B. Minimizing the Effect of Rats Constant Variations The application of the Kalman filter to the problem Of between-run variations in the rate constant assumes that the system under study follows first- or pseudo-first-order kinetics. This was the same assumption made in the development of the two-rate method in Chapter II. For the theoretical development here, it is assumed that the variable to be measured during the experiment is the absorbance of the solution, although any measurement which is directly proportional to concentration could be used. Also, it will be assumed that the measurement follows the formation of the product so that a rising exponential is observed. This is the most common case for pseudo-first- order systems. Under these conditions, the absorbance at any time t will be given by, A1 : Ar(1 - e‘"“) + B (4.8) 58 where A: is the final absorbance due to the product, it is the first- or pseudo-first-order rate constant, and B accounts for any constant background absorbance which might be present. The parameter of most interest from the analyst’s point Of view is Ar, since this will be directly related to analyte concentration and should, in many circumstances, be independent Of between-run variations in the rate constant. While eq 4.8 represents the true measurement model for the system, it is nonlinear with respect to the model parameters and cannot be used in conjunction with the linear Kalman filter. Instead, it is necessary to employ the extended Kalman filter [20,21,23,30]. This is a version Of the filter which has been modified to accommodate nonlinear systems. In its most general form, the extended Kalman filter employs a Taylor series expansion for the observation matrix and the state transition matrix, with the expansion being truncated after the first term. In this case, the measurement model is given by, 2k : hu(x11,t) + v1:(t) (4.9) and the system model is, mm = fk(x11,t) + wa(t) (4.10) The Observation matrix is now redefined as, 111. = 3.11.1.1. (4.11) 59. evaluated at x = xx. Similarly, the new state transition matrix is, __ an F1: - max (4.12) again evaluated at x = x11. With these definitions, eq 4.1 is replaced by eq 4.9, eq 4.10 replaces eq 4.2, and equations 4.4 and 4.6 are replaced by equations 4.13 and 4.14, respectively. x1: = xs" -I- K1: (2!: - h11(x11,t)) (4.13) n+1 = f11(xs, t) (4.14) Other definitions in the algorithm remain the same. For the particular case under consideration, the state parameters are assumed to be constant during a measurement sequence, so the state transition matrix, F11, is the identity matrix and does not require adaptation to the general case of nonlinearity. Only the measurement model is nonlinear. In accordance with this, the particular definitions used in the Kalman filter algorithm are as given below. At x = k (4.15) E z = A1 = A: (1 - e'“) + B = h(x,t) (4.16) 60 = [ (1-e'“) (We-h) 1 I (4.17) 1 o o F=I=010 (4.18) o o 1 The equations defining the algorithm, in order of application, are 4.3, 4.13, 4.5, 4.6, and 4.7. Although variations of the equations may have been used in actual practice (i.e. Potter-Schmidt algorithm), the sequence is essentially the same. While the extended Kalman filter is an adaptation of the ordinary Kalman filter to nonlinear systems, its characteristics are quite different from those Of the linear version. Some Of the major differences are: (1) while the linear Kalman filter is guaranteed to provide Optimal estimates, the extended Kalman filter is not; (2) the extended version can be numerically unstable and is susceptible to initial parameter estimates; and (3) the error estimates (i.e. covariance matrix) Obtained with the extended Kalman filter can only be considered reliable when convergence is achieved. Because of these characteristics, the extended Kalman filter must be applied with care and limited to systems which are fairly well-behaved (i.e.. relatively few state parameters, smoothly varying functions). It is not well-suited as a general nonlinear fitting method, but does have advantages for recursive fitting in some cases. In the work presented here, initial estimates of the state parameters were Obtained as follows. The initial estimate for the rate constant k was taken to be its nominal value in the absence of any perturbations to this system. This value was obtained from a nonlinear least-squares fit Of the first-order equation to the data. To estimate AI 61 and B, a linear fit Of the first 20 data points was employed (i.e. the initial rate method). The intercept of this fit provided B, and the slope was used calculate A: from the nominal value of k. The initial values Of the diagonal elements Of the covariance matrix, P, were taken to be large so as not to bias the algorithm, but were constrained to a certain degree (i.e. standard deviations of approximately 50% of the initial estimate) to help avoid problems of instability arising from such things as negative rate constants or absorbances. The Off-diagonal elements of P were initially taken to be zero. The magnitude of the measurement error, R, was determined from the standard error of the estimate of the nonlinear fit of a typical data set. The value of the Q matrix was zero throughout the calculations. It should be pointed out that much of the pioneering work in the application of the Kalman filter to chemical kinetics was carried out by Rutan and Brown [30,31]. Much Of the above development is an adaptation of their work. They did not, however, apply their methods to the problems considered here. C. Differential Kinetic Methods The application of the Kalman filter to differential kinetic methods is more straightforward than for the problem of rate constant variations considered in the preceding section, since a linear model is imposed. For the theoretical development, the case of two parallel reactions is considered, each of which yields an absorbing product. This model is easily extended to the case Of more analytes, but is limited here for the sake Of evaluation. 62 If first- or pseudo-first-order kinetics are assumed for the model, the absorbance of the reacting mixture at any time t is given by, A: = A1(1-e'”*) + Aa(1-e-m) + B (4.19) where A: and A2 are the final absorbances due to the products of reactions 1 and 2, respectively, in and k: are the corresponding rate constants, and B is the constant background absorbance. It should be pointed out that eq 4.19. still holds if the analytes exhibit some absorbance at the wavelength Of interest. Although the interpretation of A1 and A: is different in this case, the resulting calibration curve derived from these values is still valid. In this example, the analyst is primarily interested in the values of A1 and A2. Traditional approaches to differential kinetic methods have used a number of computational strategies [9-11]. Two Of the most popular approaches have been the graphical extrapolation method and the method of proportional equations. The former method requires the graphical determination of the initial rate Of reaction and the rate of the slower reaction after the faster reaction is complete. The latter approach is a two-point method which uses the solution of simultaneous equations to Obtain analyte concentrations (the background term in eq 4.19 is either non-existent or subtracted prior to calculation). Acceptable rate constant ratios for these two methods depend on numerous factors, including the ratio Of analyte concentrations. A complete error analysis is found elsewhere [9], but it is clear that the ultimate limitation of these two methods arises from two sources: (1) erroneous assumptions regarding the system model (eg. one reaction is complete while the other 63 is in its initial stages), and (2) the fact that only a small fraction of the data collected are used. in the computation, thereby limiting precision. Only if these limitations are removed can the full potential Of differential kinetic methods be evaluated. These weaknesses were recognized some time ago by Pausch and Margerum [15]. These authors successfully employed a method which solved the system of linear equations arising from the model using all Of the data rather than just two points. While this method effectively used all of the data, it was cumbersome to implement, especially on the digital computers of the time, since it required the manipulation of large matrices in batch mode. Margerum, Pardue, and coworkers recognized this deficiency and proposed a method for on-line regressive differential kinetic analysis [17]. This approach reduced the computational burden, but had a number of disadvantages. First, it required the use of a somewhat arbitrary weighting function to help balance the effect of bias introduced by the time and magnitude of the measurements. Second, the acquisition rate had to be Optimized for each different mixture type examined to ensure the best performance. Finally, and perhaps most importantly, the criterion for statistical Optimality with this method is uncertain. Although least squared error was used as the criterion in the theoretical development, it does not hold in the final algorithm, as evidenced by the need for a weighting function. The application of the linear Kalman filter not only alleviates the drawbacks of the methods described above, but can also provide some additional advantages. In order to establish a linear model for the Kalman filter, it is necessary to assume invariant rate constants. Use of the extended Kalman filter to compensate for rate constant variations as 64 was done in the preceding section is fraught with difficulties since the complexity Of the multicomponent model makes parameter estimation less reliable. Nevertheless, the requirement for invariant rate constants is not unreasonable, as this is normally assumed for differential kinetic methods. Given this constraint, the state vector may be represented by, x = A2 (4.20) and the measurement model in state space notation is, A1 2 = A: = [(l-e‘kl‘) (l-e"k2‘) 1] A2 + v (4.21) B where the first matrix on the right-hand side is the observation matrix, H. This definition is in accordance with the measurement model, eq 4.1. The remaining definitions for the Kalman filter algorithm are straightforward. The state transition matrix, F, is the identity matrix, since the state vector is static. Initial estimates Of the state parameters are all taken to be zero, since their values are not critical for linear. Kalman filter. The values of the rate constants for the observation matrix and the magnitude Of ‘ the measurement noise, R, can all be Obtained from a single nonlinear fit of an appropriate data set. The initial values of the diagonal elements of the covariance matrix, P, were taken to be large (several orders of magnitude larger than the expected values Of the corresponding parameters) to avoid bias, and the initial 65 Off-diagonal elements were set to zero. The Q matrix for errors in the state parameters was zero throughout. Application of the Kalman filter eliminates the two major limitations Of traditional differential kinetic methods. It uses the system model directly and is therefore free Of any erroneous assumptions related to relative reaction rates, save those arising from the model itself. It also uses all Of the data for improved precision and will provide the best estimate in accordance with a least squared error criterion. Further advantages are discussed in a later section. EXPERIMENTAL SECTION A. The Molybdophosphate Reaction In order to experimentally test the effectiveness Of the extended Kalman filter for reducing the susceptibility of reaction-rate methods to between-run variations in the first-order rate constant, the molybdophosphate reaction for the determination Of phosphate was used [32,33]. Calibration data were Obtained with five solutions with concentrations ranging from 1 ug/mL to 10 ug/mL P at three temperatures, 21.5 °C, 26.0 °C, and 31.2 °C. The temperature variations were used to bring about changes in the pseudo-first-order rate constant. Data were collected on an automated stopped-flow instrument constructed in our laboratory. Further experimental details appear in Chapter II and elsewhere [7] and will not be repeated here. 66 B. Simulated Data In order to fully evaluate the Kalman filter for differential kinetic methods, simulated data were used so that a wide variety of conditions could be tested. One hundred data points were generated in accordance with eq 4.19. The rate constant ratio, k2/k1, was varied from 1.5 to 10 in the study and data were always generated for a period Of time equal to 2/k1 (i.e. 21 for the slower reaction). To ensure rigorous testing, values of 1.0 and 0.1 were used for A1 and A2, with the order of assignment being varied to test the versatility Of the method. Differential kinetic methods are Often tested under conditions where an equal response is observed for both analytes, but this is a "best case" situation, so a ratio Of 10:1 was selected here. A value Of 0.1 was used for the background absorbance, B. After the data were generated, a Gaussian noise sequence was added (see Appendix II), normally at a level of 1% of the maximum signal. C. Data Processing All calculations were carried out on a DEC LSI-11/23 minicomputer with programs written by the author in Fortran-77. Single precision arithmetic was used for all calculations. All of the Kalman filter programs used the Potter-Schmidt square-root algorithm for improved numerical stability [20,26]. Nonlinear fitting Of the data was performed with a simplex-based fitting program also written by the author (see Appendix I). 67 RESULTS AND DISCUSSION Results are presented in this section to demonstrate that the Kalman filter is an effective tool for minimizing the effects of between- run variations in the first-order rate constant in reaction-rate methods of analysis. The Kalman filter is also shown to enhance the capabilities of differential kinetic methods. A. Minimizing the Effect of Rats Constant Variations When the extended Kalman filter is applied to data from the molybdophosphate reaction in the manner already described, estimates Of the three state parameters, A1, k, and B, are updated as each new measurement is assimilated into the algorithm. One way to characterize the Kalman filter is by Observing the evolution of parameter estimates as the measurements are assimilated. Such an evolution is shown in Figure 4.2 for a typical data set. The data set employed was that for the highest phosphate concentration (10 ug/mL P) at the middle temperature Of 26.0 °C. The plots show how each Of the parameter estimates changes during the measurement process. One can see from the figure that, for all of the parameters, the estimates evolve from the initial values to relatively stable final values, although the paths taken may be quite different. As one might expect, the background absorbance value, B, stabilizes quickly and relatively smoothly. This parameter is not closely connected to the other two and its corresponding element in the Observation matrix is as simple as possible, so it is natural for it to converge fairly early in the measurement sequence. On the other hand, there is considerable 68 Estimate of A, Estimate of k OOOOOOOOOOOOO O! J lllllll‘lll l A 0.00. : 1:: :1: 0.0 1.0 2.0 3.0 4.0 Time (sec) Figure 4.2. Typical evolution of state parameter estimates with the extended Kalman filter in accordance with the model A:=Ar(l-e"")+B. 69 correlation between the other two parameters and their evolution is closely connected. The convergence of these parameters is also not as rapid due to a problem with observability in the early stages of the reaction. Observability refers to the distinguishabilty of state parameters in the context Of the model equation and available data. For example, if we were to use k = ka + in, then It. and kb are not independently observable with the given model and data. In the same way, observability of k and A: is hindered in the initial stages of the reaction when the curve is still fairly linear. This subject is discussed in more detail later. Another feature which is apparent in the evolution Of these two parameters is the discontinuities in the curves. These small perturbations are characteristic of the instability Of the extended Kalman filter algorithm, but do not appear to be a problem since convergence was Observed in all cases. In order to Obtain more reliable parameter estimates, Rutan and Brown [30] used multiple passes of the extended Kalman filter, employing the final estimates Of one pass as initial estimates on the next. Better results can be Obtained with this iterative method, but it reduces the algorithm to a batch procedure rather than a recursive one. In this Operational mode, nonlinear least squares methods are expected to perform better still, so the justification of this approach (when a static state vector is used) is questionable. In the work presented here, a single pass was used, since the recursive approach is the expected mode of implementation. The critical test for the Kalman filter method was to extract the final estimates of A1 for all five concentrations at the three temperatures (rate constants) employed, and to used these to construct 70 calibration curves. The calibration curves generated in this manner are shown in the top half Of Figure 4.3. For comparison, the bottom half Of the figure shows calibration curves obtained with the initial rate method. Comparison Of the two sets of data clearly shows that the Kalman filter method is much less sensitive to variations in the rate constant, since the slope Of the calibration curve is less dependent on temperature. The Kalman filter results also exhibit good linearity. There is, however, a small correlation Of the slope of the calibration curves with temperature. From nonlinear fitting Of the reaction curves, part of this correlation was found to be due to a shift in the equilibrium absorbance with temperature. Other reasons must also contribute, however, particularly for the large low temperature deviation. The reasons for the deviation of the calibration curves from the ideal case can be appreciated, at least in part, by comparing Figures 4.4 and 4.5. Figure 4.4 shows the results Of fitting a first-order model to the data by nonlinear least squares. The reaction curves are those for the high phosphate solution at each temperature. There is a good fit in each case, indicating good first-order behavior. Figure 4.5 shows fits Obtained with the extended Kalman filter. The best results are Obtained for the middle temperature, which is expected since the initial estimates assumed this rate constant. A poorer fit is Observed for the high- temperature data and the worst fit is obtained for the low-temperature data. Thus, the Kalman filter is not exactly a recursive analog of nonlinear least squares and this is evident in the final result. Another useful perspective is Obtained from Figure 4.6, which shows the evolution Of A1 for the three temperatures at the highest phosphate 71 0-5 " EKF _ o T=31.2°C 0‘4 ' + T=26.0°C .. x T=21.5°C 0.3 -- Q.- < I 0.2 -.- 0.1 -- / 0.0 - A ‘ Initial Rates 0 0 0) \ a :5 .3 O n: E E / l I 1 1 . 1 4|.l.l I 0.0..:.=.=.=..:.u1u i 0. 1. 2. 3. 4. 5. 6. 7. 8. 9. 10. Concentration (ppm P) r Figure 4.3. Comparison of calibration curves generated with the extended Kalman filter (EKF) method and the initial rate method at three temperatures. 72 0.5 - T=31.2°C A . - ' Mo 3 0 4 T-26.0 c < . v , . . o 0.3 - T=21.5°C 0 c . O D 0.2 - L. 0 . U) .0 < 0.1 - 0.0 | i i ‘r : : § ; J] 0. 1000. 2000. 3000. 4000. Time (msec) Figure 4.4. Fits of molybdophosphate reaction curves at three temperatures to a first-order model by nonlinear least squares ([P]=10 nslmLio 73 O U! I ‘ ~ For. ' “fl ,9 W.f:.jd".: ‘ w‘ , A O o 4 - 3’“...- . D 17-", V at" ’ S, d y’ ' :3...\'Mw .“'o ‘ I...‘ a) O . 3 c- ff *w\"‘" '5' o I’ r C " f C ’3'. Q 0 . 2 " :' . L. ’l ‘ o d (D .0 < O . 1 '1 . ; «(v . ' f 00 l i l i l i If i l 0. 1000. 2000. 3000. 4000. Time (msec) Figure 4.5. Fits of molybdophosphate reaction curves at three temperatures to a first-order model by the extended Kalman filter ([P]=10 us/mbl- 74 0.7 r . e . . . S 0.6d T=31.2°c .. 3 * * 0.5 1 .2" T=26.0°C . ._ 0.4 1 O o 0.3 T 215°C ‘ cu = . T E 0.2 - g; 0.1 - - 0.0 4. = i : i : 0.0 1.0 2.0 3.0 4.0 Time (see) Figure 4.6. Evolution of final absorbance estimates at three temperatures. 75 concentration. Note that the curves are approaching the correct equilibrium value in all cases, but convergence has 'not yet been achieved for the high and low temperatures. The relatively large deviation of results for the low-temperature case may be investigated further. Qualitatively, it is clear that there is a high correlation between A: and 1: early in the reaction when the reaction curve is very nearly linear. As indicated earlier, this can cause problems concerning the independent observability of the two parameters. The low temperature data consist of fewer half-lives and are therefore more susceptible to this problem. Hence, this set of data exhibits the largest deviation. These qualitative observations can be made more quantitative through an entity associated with the Kalman filter known as the information matrix. For a static system (i.e. the state parameters are constant) the information matrix is defined according to eq 4.22 [21,30]. In(m) = [ HT(1')-R‘1-H(i) ] (4.22) 1 use: In this equation, 1n(m) is the information matrix after measurement m. If the determinant of the information matrix is zero, then at least two of the state parameters are not independently observable from the measurements processed to that point. The magnitude of the determinant of the information matrix is therefore directly related to the observability of the system. Figure 4.7 shows the evolution of det(In) for the three data sets shown in Figure 4.5. The most observable system is the high-temperature case since the data here represent the 76 001(1) (arbitraty units) T I 0.0 1.0 2.0 3.0 4.0 Time (see) Figure 4.7. Evolution of the determinant of the information matrix three temperatures. at 77 largest number of half-lives. The results for the low—temperature data, on the other hand, are coincident with the abscissa throughout the measurement sequence, making the observability of this system subject to question. While det(In) does have a non-zero value for the low- temperature case, it is clear that the information content of these data is less and this accounts for the large deviation of the low-temperature calibration curve. Some additional points are worth noting. First, for the linear Kalman filter, the covariance matrix, P, is a useful entity, since it reflects with good reliability the errors in the estimated parameters. For the extended Kalman filter, this was not found to be the case. During the evolution of the state vector, the actual error in the parameter estimates often exceeded the predicted error by a considerable amount. Only upon convergence of the filter did the two quantities attain similar values. Finally, the suboptimal behavior of the extended Kalman filter for modeling nonlinear systems may be attributed to two sources [19]. It may be due to the nonlinear form of the model itself, or due to poor predictions by the observation and transition matrices due to truncation of the Taylor series expansion after the first term. The latter problem may be mitigated by the use of second- or third-order extended Kalman filters which truncate the Taylor series after the second or third terms. This improves the predictive abilities of H and F. A second-order extended Kalman filter was used at one stage of this work, but there was virtually no change in the results obtained. Therefore, the nature of the model nonlinearity must be the dominant factor. 78 B. Differential Kinetic Methods A typical set of simulated data generated for this study is shown in Figure 4.8. In this particular case, the ratio of rate constants used was Ira/k1 = 3, with A1 = 1.0, A2 = 0.1, and B = 0.1. It is clear from the figure that the two reactions are not easily distinguished and traditional differential kinetic methods, particularly graphical extrapolation, become less effective under' these conditions. The Kalman filter can be effectively employed, however. i To test the utility of the Kalman filter for differential kinetic methods of analysis, data similar to those shown in Figure 4.8 were generated for Its/1n ranging from 1.5 to 10.0 in steps of 0.1. In each case, the data were analyzed with the Kalman filter and the final estimates of A1 and A2 were recorded. The results of this study are summarized in Figure 4.9. The dashed lines represent the actual values used for A1 and A: (1.0 and 0.1 respectively), while the points indicate the values computed with the Kalman filter. One can see that there is good agreement over the entire range of ratios used, although larger errors are evident for small rate constant ratios. To ensure a completely rigorous evaluation, the order of assignment of A1 and A: was reversed, i.e. A1 = 0.1 and A2 = 1.0, and the calculations were repeated. These results are shown in Figure 4.10. Again, good agreement is observed in every case. It will also be noted that the error sequences observed are the same as those in Figure 4.9, except that they are reversed. An interesting property of the linear Kalman filter is that, since the observation matrix (H), the covariance matrix (P), the Kalman gain vector (K), and the state transition matrix (F) are not functions of the 79 Response l , : .0 0.0,::::§:::: e ..I 0.0 0.5 1 15 2.0 Time (in units of 1/k1) l L P l l _l I l ' ' C Figure 4.8. Typical simulated data for differential kinetic analysis by the Kalman filter in accordance with the model: A1 = A1(l-e"“‘) + A:(l-e'”‘) + B. Values used were A1=l.O, Aa=0.1, kzlk1=3, 1% noise. 80 N I N < «(r- : 1 . O .—------3':.'O..’a..fl‘epraqfi..flv.;“rad-9HHQ‘.:.J.‘~.~;Q°O~FUQQd‘fihl 0 cu- f O . 8 ~- 3 T ,, o . 6 ~- 0 w up D s 0.4-w '3 .. U - O . 2 '- O .E 1.-----gc=5;.&4'.~;b.~."4°—:.£u£0.’1¢k~ub.~?q“0";arfiq‘,.«l‘4g..--:'..'.°.. LL ' l i l I l O s O 1r T l 'L I : —[ l J 0.0 2.0 r 430 6.0 8.0 .1013 k2/k1 Figure 4.9. Results of differential kinetic analysis by the Kalman filter as a function of rate constant ratio. A1=l.O, A2=O.l. 81 lo I l s s s _ -- - - o c :5 ..;O-..~ 9.5.2Lfiu'.u-.gggudbqeflne'.c.c 3"qu —"...-. ‘fl--....... . . . . . and A2 2: l l P 0.81- J V O O) 1 i l I L F )- O N I O . O D.... c -- ..'- '5'.:‘Cepr‘qfifiw..3rond-r‘-dq‘.-'.'.dqo~'..'.~..vfi...;ghl s . s 0,0 : i t i : A} : i : 4 0.0 2.0 4.0 6.0 8.0 10.0 k2/k1 Final Eshnates of A1 0 .p 1 I Figure 4.10. Results of differential kinetic analysis by the Kalman filter as a function of rate constant ratio. A1=0.1, A2=1.0. 82 state parameters or measurements, the evolution of the covariance matrix (i.e. the error estimates for the state parameters) can be determined without any experimental results [20,22]. Only the system models and the magnitude of the measurement error need to be defined. This property is extremely useful since it allows a proposed application to be evaluated before real experiments are conducted. In performing such an evaluation, equations 4.3, 4.5, and 4.7 are used. In the present study, this approach was used to assess the potential of the method by examining the final value of the covariance matrix for each rate constant ratio used. As a specific example, the standard deviation associated with parameter A: (when A2 = 0.1) was recorded at each ratio. The predicted error, in the form of lo' boundaries, is plotted in Figure 4.11 (solid, curved lines). The points plotted represent the magnitude of the error actually observed. It will be noted that the majority of points fall within the 16' boundaries, indicating that the predicted error is reliable. More importantly, however, the plot indicates the limitations of differential kinetic methods under the measurement conditions imposed in this example. As the rate constant ratio becomes smaller, the predicted error shows relatively little change until a ratio of about 3 is reached, at which point the errors begin to increase rapidly. This observation is corroborated by the experimental points as well. While the nature of this behavior will depend to some degree on the magnitude of the experimental noise and the measurement interval used, it does show that differential kinetic methods are not currently being pushed to their full potential. Margerum, Pardue and coworkers reported that their regressive differential kinetic analysis method was useful for rate constant ratios 83 "5 0.05 N 0u04 “ (L03 5 o 0.02 h 5’ 0u01 , , .. 3 _ 0 (L00 . ..- .e _E_-o.01 13‘0““ £4.03 _-o.04 S-o.os. :1' = : r: : i = : 15 0.0 2.0 4.0 6.0 843 10.0 l"Iv/"1 Figure 4.11. Comparison of actual (points) and predicted (curved solid lines) errors for the Kalman filter approach to differential kinetic analysis as a function of rate constant ratio. 84 of 1.5 or less [17]. This figure is misleading, as the conclusion was reached using "perfect simulated data". Since "perfect data" is an oxymoron in reality, the value of the simulation is questionable. This is further demonstrated by the experimental results they presented. These showed significant errors (several percent) even when the ratio of rate constants in a binary mixture was about 15 [15,17]. Currently, no chemical systems have been employed in this study. Experiments are now underway in our laboratory to provide data for this purpose, however. There are numerous advantages of the Kalman filter approach to differential kinetic methods. First, as already noted, it uses all of the data and is consistent with kinetic models employed. Second, it is guaranteed to provide the best estimate of the concentrations possible in accordance with the model, with the possible exception of nonlinear least squares methods. The computations carried out with the linear Kalman filter produce essentially the same results as a linear least squares analysis, which is generally considered the optimal approach. The third advantage of the Kalman filter is that it is a recursive approach rather than a batch method (i.e. traditional linear least squares) and is better suited to reaction-rate methods. Computations are simplified, it allows estimates to (be made as the experiment progresses (for relatively slow reactions), and the number of measurements employed is not constrained to a particular value. Finally, the Kalman filter provides error estimates for the concentrations, and can easily be used to evaluate the potential of new differential kinetic methods before their implementation. 85 C. Additional Applications There are numerous additional applications of the Kalman filter to reaction-rate methods of analysis. For example, Corcoran and Rutan [34] have investigated the potential of the Kalman filter for removing the effect of rate constant variations during a run (as opposed to between runs) by simultaneously making temperature measurements. Another relatively simple application is considered here which may have some utility. Although no studies of this application have been conducted, it is included as a possible direction of future research. The problem of kinetic interferences in reaction-rate methods is a general one, affecting all methods. Often such an interference takes the form of a slow reaction or reactions proceeding in the presence of the reaction of interest (for example, the interference of cepha antibiotics on the clinical determination of creatinine by the Jaffé method [35]). If the interfering reaction is sufficiently slow, as it often is since the reaction of interest is generally more favorable, then it can be represented as a linearly increasing background signal. For a first- or pseudo-first- order reaction, this is represented by eq 4.23. A: = A:(l-e"“) + ct + B (4.23) In this equation the term "ct" represents the linearly increasing background, where "c" is a constant. The contribution of the linear term may vary from sample to sample. If the rate constant is invariant, then eq 4.23 represents a linear system which is amenable to solution with the Kalman filter. This would provide a method to detect and 86 account for the linearly changing background, and may be a useful approach to systems where kinetic interferences are a problem. CONCLUSIONS It has been shown that the extended Kalman filter can reduce the effect of between-run variations in the first-order rate constant on reaction-rate methods of analysis. This was demonstrated for phosphate determinations at three different temperatures, conditions which gave rise to large rate constant differences. The method as implemented was completely stable, but did show some dependence on initial parameter estimates. This, in addition to the smaller number of half-lives used for the low temperature case, causes a small correlation between the temperature and the slope of the calibration curves, although good linearity was maintained at each temperature. For smaller, more reasonable fluctuations in the rate constant, the variation in results is expected to be insignificant. However, implementation of the extended Kalman filter may be hindered in practice by its perceived complexity. In addition, a good estimate of the experimental noise is required for the method. The application of the Kalman filter to differential kinetic methods of analysis should significantly extend the limits of traditional approaches to these methods. In addition, it does not suffer from the drawbacks of earlier multipoint methods. The required rate constant ratio with this approach may be as low as 3 or less. This would extend the number of systems to which differential kinetic methods could be 87 applied and also increase the number of components which could be simultaneously determined. The linear Kalman filter ' method is very stable, is straightforward to implement, uses all of the available data, and provides error estimates for the concentrations. It does, however, require reasonably accurate values for the rate constants and measurement noise, and is limited (so far) to systems following first- order kinetics. 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 88 LIST OF REFERENCES J.B. Davis and B. Renae, Anal. Chem., 51, 526-528 (1979). F.J. Holler, R.K. Calhoun and S.F. McClanahan, Anal. Chem., 54, 755-761 (1982). R.G. Cornell, Biometrics, 18, 104—113 (1962). P.B. Kelter and J.D. Carr, Anal. Chem., 51, 1825-1828 (1979). P.B. Kelter and J.D. Carr, Anal. Chem., 51, 1828-1834 (1979). P.B. Kelter and J.D. Carr, Anal. Chem., 52, 1552 (1980). P.D. Wentzell and S.R. Crouch, Anal. Chem., 58, 2851-2855 (1986). G.E. Mieling and H.L. Pardue, Anal. Chem., 50, 1611-1618 ((1978). H.B. Mark, Jr. and G.A. Rechnitz, "Kinetics in Analytical Chemistry"; Wiley-Interscience, New York, 1968; Chapters 6 and 7. S.R. Crouch in "Computers in Chemistry and Instrumentation", H.D. Mattson, H.B. Mark, Jr. and H.C. MacDonald, Eds.; Marcel Dekker, New York, 1973; Vol. 3 ("Spectroscopy and Kinetics"), pp 107-207. D. Perez-Bendito, Analyst, 109, 891-899 (1984). C.G. Kircher and S.R. Crouch, Anal. Chem., 55, 248-253 (1983). E. Pelizzetti, G. Giraudi and E. Mentasti, Anal. Chim. Acta, 94, 479-483 (1977). D. Benson and N. Fletcher, Talanta, 13, 1207-1209 (1966). J.B. Pausch and D.W Margerum, Anal. Chem., 41, 226-232 (1969). J.B. Pausch and D.W Margerum, Anal. Chem., 41, 233-238 (1969). 8.6. Willis, W.H. Woodruff, J.R. Frysinger, D.W. Margerum and H.L. Pardue, Anal. Chem., 42, 1350-1355 (1970). R.E. Kalman, Trans. ASME, Ser. D: J. Basic Eng., 82(March), 35-45 (1960). A. Gelb (Ed.), "Applied Optimal Estimation"; MIT Press, Cambridge, MA, 1974. R.G. Brown, "Introduction to Random Signal Analysis and Kalman Filtering"; John Wiley and Sons, New York, 1983. S.D. Brown, Anal. Chim. Acta, 181, 1-26 (1986). 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 89 W.S. Cooper, Rev. Sci. Instrum., 57, 2862-2869 (1986). 8.0. Rutan, J. Chemomet., 1, 7-18 (1987). W.H. Press, B.P. 'Flannery, S.A. Teukolsky and W.T. Vetterling, "Numerical Recipes"; Cambridge University Press, New York, 1986; Chapter 14. H.N.J. Poulisse, Anal. Chim. Acta, 112, 361-374 (1979). T.F. Brown, D.M. Caster and S.D. Brown, Anal. Chem., 56, 1214-1221 (1984). H.W. Sorenson, IEEE Spectrum, 7(July), 63-68 (1970). T.A. Brubaker, R. Tracy and C.L. Pomernacki, Anal. Chem., 50, 1017A-1024A (1978). P.D. Wentzell and S.R. Crouch, Anal. Chem. accepted for publication. S.C. Rutan and S.D. Brown, Anal. Chim. Acta, 167, 23-37 (1985). 8.0. Rutan and S.D. Brown, Anal. Chim. Acta, 175, 219-229 (1985). W. Rieman and J. Beukenkamp in "Treatise on Analytical Chemistry", I.M. Kolthoff and P.J. Elving, Eds.; Wiley, New York, 1961; Part II, Vol. 5, pp 348-351. A.C. Javier, S.R. Crouch and H.V. Malmstadt, Anal. Chem., 41, 2139-243 (1969). C.A. Corcoran and 8.0. Rutan, Anal. Chem., submitted for publication. CHAPTER V Modeling Chemical Response Surfaces with the Kalman Filter "I would like to emphasize strongly my belief that the era of computing chemists, when hundreds if not thousands of chemists will go to the computing machine instead of the laboratory, for increasingly many facets of chemical information, is already at hand. " - Robert S. Mulliken Modeling of response surfaces has long been recognized as an important aspect of research in chemistry and other disciplines. The ability to predict the magnitude of some response, such as percent yield or absorbance, as a function of experimental variables is often the primary goal of a research study. Ultimately, one hopes to obtain a mechanistic model which adequately reconciles experimental results with theoretical expectations. Where such a mechanistic model cannot be realized, an empirical model based on an appropriate function is often used. Empirical response surface modeling has been used for many years [1-3] and is still the subject of much interest, as evidenced by a recent book on the subject by Box and Draper [4]. The importance of empirical modeling arises from the need to describe multivariate surfaces which are not easily treated in any other manner. Such situations are a common occurrence and response surfaces in analytical chemistry are a case in point. Quite often the analytical chemist wishes to optimize or model the response (sensitivity, selectivity, analysis time, etc.) of an 90 91 experimental system which is a function of several variables (temperature, pH, reagent concentration, etc.). It has been established [5] that the univariate approach to optimization is inadequate where variable interaction occurs and alternate methods must be used. Simplex Optimization [5-8] has proven successful for response optimization, but fails to describe the nature of the surface in a direct manner. When the response surface is relatively smooth, empirical modeling can yield useful descriptive results. Numerous examples of empirical response surface modeling can be found in the analytical chemistry literature [9-13]. The typical approach to empirical response surface modeling involves a linear least squares fit of a multivariate polynomial. Normally, this is a second-order polynomial with appropriate cross terms for each variable [3,4]. Higher-order polynomials are generally too sensitive to minor variations in the data and may mask the major trends. Also, the number of terms (basis functions) increases rapidly with the degree of the polynomial. Cross terms in the model are essential in accounting for interactions between variables. The coefficients of the modeling equation can be determined by solution of the normal equations using Gauss-Jordan elimination or other standard methods [14-16]. Currently, for systems with a large number of coefficients 0 10) the preferred method of solution is singular value' decomposition (SVD) since it is more stable and less prone to round-off errors [16]. In this chapter, an alternative method for response surface modeling based on the Kalman filter is described [17]. The Kalman filter [18-20] was introduced in Chapter IV and has been applied to a number 92 of areas in analytical chemistry in recent years [21-23]. We demonstrate that, in modeling applications, the Kalman filter yields results which are virtually identical to those of singular value decomposition, even for relatively large systems, and offers a number of advantages over SVD for common modeling problems. The Kalman filter method is applied to two flow injection analysis (FIA) systems and an electrochemical system. The advantages of the Kalman filter approach and its relationship to simplex optimization are also discussed. THEORY The theoretical aspects of the discrete Kalman filter have been presented in Chapter IV and elsewhere [18-23], and will not be repeated here. The symbols used in reference to the algorithm are the same to ensure consistency. This section begins with a discussion of the conditions under which the Kalman filter yields the same results as linear least squares. This is followed by a description of the nature of the models used in this work. The theoretical aspects of singular value decomposition are well established [16] and are not treated here. A. Relationship to Deterministic Least Squares For polynomial response surface modeling, the elements of the Kalman filter state vector are the coefficients of the polynomial model. As already mentioned in Chapter IV, the recursive procedure defined by the linear Kalman filter algorithm will provide an optimal estimate of the 93 state vector in accordance with a minimum mean-square error criterion. Polynomial response surface modeling is normally performed by linear regression methods, which obtain the optimal estimate of the coefficients (i.e. the best fit) according to a minimum residual sum of squares, or least squares, criterion. Strictly speaking the two criteria for optimality, least squares and minimum mean-square error, are not identical, although the practical differences are often negligible. The distinction is quite subtle, and is probably best exemplified by the difference between fitting data to an equation (least squares) and performing ensemble averaging (minimum mean-square error). Under certain conditions, the results obtained with the Kalman filter are identical to those obtained by a least squares method. For this reason the Kalman filter is often referred to as a "least squares filter" and has been used for linear least-squares estimation [24—26]. This is the manner in which it is to be employed here. The conditions under which the Kalman filter yields the same results as deterministic least squares are: (1) the system state vector is assumed to be constant (i.e. F: = I), (2) the measurement sequence defines an overdetermined system of linear equations (analogous to a non-singular matrix in deterministic least squares), and (3) no prior knowledge of the state vector is assumed (i.e. the diagonal elements of Po’ are assumed to be infinite). All of these conditions are met here. B. Model Definition An example is presented here to demonstrate how least squares modeling is carried out with the Kalman filter. Suppose we have a series of measurements of some dependent variable, z, as a function of 94 an independent variable, y. We wish to fit this data to the second- order polynomial, sy'2 + by + c (5.1) N II where a, b, and c are the parameters to be obtained from the fit. In this case, the state vector is, a x = b (5.2) c and the observation matrix is, H = [ ya y 1 ] (5.3) Since the state vector is constant, the transition matrix, F, is the identity matrix. For similar reasons we can assume Q = 0. Since only one measurement is made at each interval, we also have R = 823, where 3:3 is the variance associated with each measurement (8:2 is an estimator for 612). The values of R assume the same role as the weighting factors in deterministic least squares and, in this sense, only the relative magnitudes are important. Ideally, the initial value for the error covariance matrix should be, 95 Po" = (5.4) oog ogo goo Since infinite values are not a computational reality, the diagonal elements can simply be very large relative to R [24]. Typically, we employed values of about 1010 for R = 0.1 to 1. Finally, the initial value of the state vector was taken to be zero, although this should not be important if the elements of the initial covariance matrix is large enough. The above example contains all of the definitions necessary to carry out the least squares fit with the Kalman filter. Beginning with eq 4.3, the fitting parameters can be recursively estimated for each new value of z (the modeled response) and the final values will be optimal in the least squares sense. This is essentially the procedure used in this work, except that larger, multivariate polynomial expressions were employed rather than the univariate quadratic of eq 5.1. For example, the model used for the two-variable case is, z : ayg + by: + 0m2 + dyz2 + ey1y2 + f (5.5) The state parameters with this model are a, b, c, d, e, and f, and the observation matrix consists of the terms y1, ya, 3'12, ya”, y1y2, and 1. This can be directly extended to more variables, with the number of cross terms increasing rapidly. 96 EXPERIMENTAL SECTION In this section the chemical systems used for the modeling studies are briefly described and the computational aspects of the modeling procedures are discussed. Since details of the three chemical systems have appeared in the literature, along with the pertinent data, only the essential information is given here. A. Analytical Systems The first system modeled was an FIA system for the determination of glucose by an immobilized enzyme/Trinder Reaction procedure reported elsewhere [27]. Optimization of this system was achieved with the composite modified simplex method of Betteridge et al. [7,8]. The data resulting from this optimization were modeled with the Kalman filter and with singular value decomposition. The optimization procedure involved seven experimental variables, including the pH of the carrier and reagent streams, the flow rate, the temperature of the immobilized enzyme single bead string reactor, and three reagent concentrations. The modeled response, the time-corrected absorbance ratio (TCAR), was a composite function of time and analytical peak height to ensure an adequate compromise between sensitivity and analysis time, and is given by) TCAR 1'- [80 8/(60 S + (41)] Aexp/Abase (5.6) In eq 5.6, Aexp is the peak absorbance for a given set of conditions, Abase is the peak absorbance for the reference set of conditions, and tp 97 is the time from injection to the peak maximum. A total of 69 points were used in the fit. The modeling equation was a seven-variable quadratic function with cross terms for each pair of variables for a total of 36 fitting parameters. The second implementation of the Kalman filter algorithm utilized results from a previously reported optimization for the determination of isOprenaline by FIA [28,29]. In this paper five experimental variables were considered: the reactor coil length, flow rate, pH, reagent concentration, and sample loop volume. The data set consisted of 168 iterative univariate experiments and 37 simplex experiments. Peak height values in mm were used as the response function. For this system, a five-variable quadratic with cross terms was used as the modeling function. The final application of the Kalman filter modeling approach employed simplex optimization data from the polarographic determination of uranyl ion [7,30]. The results of 26 experiments, in which the concentrations of nitric and hydrochloric acids were varied, were used. Since the merit of each set of experimental conditions depended on both the height and shape of the polarographic wave observed, a composite response function incorporating both parameters was used in the optimization and in the modeling studies presented here. The form of this function is, R = h case (5.7) where R is the response function, h is the height of the catalytic wave and 8 is the angle between the horizontal and the tangent to the top of 98 the wave. The model function was a simple two-variable quadratic with a cross term. B. Computational Aspects Two algorithms were employed for response surface modeling, the Kalman filter and singular value decomposition (SVD). The latter was used for comparison studies. The Kalman filter algorithm was implemented in accordance with the equations present in Chapter IV. Calculations for the two FIA systems were carried out on a DEC LSI 11/23 minicomputer with programs written by the author in Fortran-77. For the electrochemical system, programs written in Microsoft QuickBASIC were used on an IBM PC compatible microcomputer. Double precision arithmetic was used in both cases. In addition, all independent experimental variables were normalized to fill the range 0 to 1. This ensured that the model coefficients generated reflected the significance of the corresponding basis function rather than the scale of the corresponding variable. The response values (dependent variable) were not adjusted. The SVD algorithm used in calculations for the seven variable FIA system was adapted from Fortran routines provided by Press and coauthors [16]. Minor modifications were made to accommodate aspects of the particular modeling problem. The calculations were carried out on a DEC LSI 11/23 and also incorporated double precision arithmetic and normalized experimental variables. 99 RESULTS AND DISCUSSION In this section, each of the three systems modeled is used to demonstrate a particular point. The seven-variable FIA system is used to compare the results of the Kalman filter with those of SVD for a complex model; the five-variable FIA system employs a large number of points and is used to demonstrate some of the aspects of this approach to modeling; and the polarographic system is used for a direct comparison of three-dimensional surfaces where some of the points must be ignored in the model. The advantages of the Kalman filter method over SVD are also discussed, and its relationship to simplex optimization is examined. A. Seven-Variable FIA System. While the theoretical relationship between Kalman filtering and deterministic least squares is well established, no examples of the farmer to the modeling of complex multivariate systems has been encountered. Therefore, we first wished to demonstrate that the results obtained with the Kalman filter are the same as those obtained through traditional least squares methods such as SVD. The seven-variable FIA system for the determination of glucose was well suited as a test since it required a model with 36 basis functions. Success with this would ensure that the method would be reliable for systems of lower dimensionality. The data from previously conducted optimization experiments were processed by both the Kalman filter and SVD algorithms. To compare the results," the 36 coefficients obtained from the SVD algorithm were plotted against those obtained with the Kalman filter, as shown in 100 Figure 5.1. The coefficients ranged from -196 to 249 and a least squares fit (by traditional means) gave a slope of 1.0009 and an intercept of -0.0009. The correlation coefficient of the fit was 1.00000. Since even small differences in the coefficients can cause significant variations in the response surface, another useful test of the equivalence of the methods is to compare the residual sum of squares for the modeled surface, which was 327.627 in both cases. Thus it is clear that the results obtained with the Kalman filter are virtually identical to those obtained with the more traditional SVD approach, even for this large system. Unfortunately the quadratic model obtained for the seven-variable FIA system was inadequate for predictive and diagnostic purposes. This was evident from several points of view: ( 1) poor convergence of the model coefficients as the final measurements were assimilated; (2) the absence of a predicted optimum (there was a saddle point instead); and (3) generally poor predictive and morphological qualities for the model. The failure was due to an insufficient number of experimental points (less than twice the number of coefficients estimated), and was not because of any deficiency in the method of modeling the response surface. For such a complex system, many more points (and/or an alternate model) would be required for satisfactory modeling. Reducing the dimensionality of the problem was considered, but not enough variables could be reasonably rejected to make this possible. Nevertheless, the results demonstrated that SVD and the Kalman filter yield equivalent models, so only the latter was applied to the two systems that follow. 101 “a .2 .2 a: Q) 0 K) K .2 it Q U E D k -200 = t = i = t = : = i -200 —100- 0 100 200 300 SVD Coefficient Figure 5.1. Comparison of coefficients determined by the Kalman filter and singular value decomposition (SVD) for the seven-variable FIA system. 102 B. Five-Variable FIA System This system was chosen as a further demonstration of modeling by the Kalman filter because it still required a relatively complex model but consisted of a larger number of data points than the seven-variable FIA system discussed above ([number of data points]/[number of coefficients] as 10). In addition, 168 of the 205 points were obtained from univariate optimization data and therefore were more evenly distributed over the surface. One difficulty associated with this type of modeling is that the chemical system under study usually only approximates a quadratic surface at best. Gaussian surfaces generally provide a better approximation, but unfortunately constitute a nonlinear system and are therefore incompatible with the Kalman filter. The extended Kalman filter can be used for nonlinear systems [20,21], but is not as well- behaved as the linear version, particularly for large state vectors, and requires good initial parameter estimates. If the actual response surface is radically different from the form of the model, the predictive capabilities are diminished. Therefore, it is necessary to test the stability and predictive abilities of the final model equation. To do this, the Kalman filter was first used to generate the model coefficients using only the 168 points from the univariate experiments. The remaining 37 points, obtained from the simplex optimization, were used to test the model. The simplex data were not entirely representative of the data space, favoring the region around the optimum, but were a convenient test set since points fell between the regions used to generate the model. Also, the predictive abilities of the model are the most critical near the optimum. 103 The Kalman filter fit of the 168 univariate data points yielded a standard error of estimation (SEE) of 9.3 mm over a response range of 0 to 192 mm. The SEE is a common figure of merit for regression methods [31,32] and is defined according to, SEE = 2 (3’1 "J’cglc)2 (5-8) Ncal‘ Npar where y. is the observed response, ycaic is the calculated response, Ncsl is the number of points in the calibration set (Neal 2 168), and Npar is the number of fitted parameters (Npsr = 21). The summation is over all of the points in the calibration set. The SEE is a measure of the deviation of the calculated responses from the observed responses and ideally should be approximately equal to the standard deviation of the measurements (about i 1 mm in this case) for an accurate model. For empirical response surface modeling this seldom happens due to inadequacies in the model and normally the model is considered useful if it reliably reflects the trends in the data. Thus, the value of the SEE obtained in this case is quite reasonable. The standard error of prediction, SEP, is used to evaluate the predictive abilities of a model [32] and is defined by, SE? = 3 (YJ-YCalc)2 (5-9) Ntest where y: is the measured response of the test point, ycanc is the calculated response, and Ntest is the number of points in the test set (Num 2 37). The summation in this case is over the test set of data. 104 Ideally the SEP should be equal to the SEE. A value of SEP = 19.9 mm was found for the test data here. Although this is about twice the value of the SEE, it is still reasonable since many of the test points were in regions for which there were no points in the calibration set. Since the mean value of test set responses was 136 mm, the SEP also indicates that predictions can be made with some reliability. When all 205 points were used in the fit, the resulting SEE was 11.1 mm, indicating some perturbation in the model was brought about by the inclusion of the simplex data set. This is anticipated when there are errors in the form of the model, but predictive abilities should be improved. One advantage of using the Kalman filter for this type of modeling is that its recursive nature allows evaluation of convergence without the use of a test set. By observing the evolution of the parameter estimates as each new measurement is assimilated and examining the reliability of the predicted measurements, the stability of the model can be dynamically assessed. Figure 5.2 is an example of the results obtained with the model for this system. It is included to illustrate the capabilities and limitations of this type of modeling. Of course, direct comparison of experimental and modeled surfaces in six dimensions is impossible. Even reduction to a three dimensional representation does not make a direct comparison possible, since the variables which are held constant may have many levels. Nevertheless, trends may be compared. The top two plots in the figure (a and b) show how the experimental univariate response changes as a function of coil length and ferricyanide (reagent) concentration. The different sets of curves (symbols) correspond to 105 i. N 8 8 PEflK HEIGHT (mm) s a PEflK HEIGHT (mm) $0 450 0 1 2 J 4 s o FERRICYHNIDE (an use 230 ‘330 COIL LENGTH (cm) C. 71) 61)- 3 . \E— 5.0 ‘\ 1?0 It. . :9 4.0 d\ t d\ Q g. 3.0 A j K .. E 2.0 - N 4 1.0- (L0 . a I fir 200 300 400 Coil Length (cm.) I 100 Figure 5.2. Response curves for the determination of isoprenaline by FIA. (a) Experimental univariate responses as a function of coil length (other variables at four levels, as indicated by symbols). (b) Experimental univariate responses as a function of ferricyanide concentration. (a) Contour map of the two-variable response surface in accordance with the derived second-order model. Contours represent 10 mm differences in in response from 100-180 mm. Other variables are held at their optimum values. 106 different values for the other experimental parameters (pH, flow rate, sample loop volume). The contour plot in Figure 5.2c shows the two-variable surface generated when the other experimental variables are held constant at the optimum values determined by the simplex procedure [28]. For the purpose of comparing trends, the univariate responses may be viewed as "slices" of the modeled surface parallel to one of the axes. Note that the model agrees quite well with the behavior exhibited when the coil length is changed. For low responses, the response increases in a fairly linear fashion with coil length, but at higher responses a maximum is observed which shifts to shorter and shorter coil lengths. The agreement with the effect of changes in the reagent concentration is also good, but shows some discrepancies. The univariate plots show an exponential rise in the response with the reagent concentration which levels off at high concentrations. The contour plot also exhibits such a rise and leveling off, but since we are trying to fit an exponential with a quadratic model, a decrease in response is observed after the maximum is reached. This is a common deficiency of quadratic models. It should be pointed out that the region where the decrease in response is observed is one in which few or no data points were available. As long as these limitations are kept in mind, this method of modeling response surfaces can be useful. Contour plots such as that in Figure 5.2c are useful for predictive and diagnostic purposes. For example, the figure clearly shows the negative correlation between coil length and reagent concentration, i.e. the optimum coil length becomes shorter as the ferricyanide concentration is increased. This is not completely evident from the iterative univariate plots, since several variables are changed for each 107 curve. Such conclusions can also be reached by examining the coefficients of the cross terms in the model (relatively large and negative in this case). 0. Two-Variable Polarographic System As a final example of modeling with the Kalman filter, a two-variable system is presented so that a direct comparison of the actual and modeled response surfaces may be made. The polarographic determination of uranium (VI) in the presence of nitric and hydrochloric acids was used as the model system [7,30]. The response function of this system was chosen to reflect both the sensitivity and definition of the polarographic wave, as already discussed. Such a response was found to be affected by the concentration of both acids in a .complex manner. The actual response surface obtained from the optimization experiments is shown in Figure 5.3. Note that although the units for the response are uA, the response shown is actually the composite function mentioned above rather than the diffusion-limited current. The surface shows a sharp optimum which is somewhat elongated in one direction. There also appears to be a fairly strong interaction between the two acid concentrations (i.e. total [Ht] may be a fundamental factor). The surface exhibits tendencies which are more Gaussian than quadratic. Because of this, attempting to model the surface with a second-order polynomial will flatten the peak maximum in an attempt to accommodate all of the points. To minimize this effect, all of the points with responses less than 100 (points 1,2,3,7,8, and 9) were not used in the fit. This modification was trivial with the Kalman filter method, since it 108 ’°°" 4.3%: ° "/ .. "1%.? 2° 24 0 10 20 30 Figure 5.3. Three-dimensional display of points obtained from the simplex optimization of the polarographic determination of uranium (VI). The labels on the data refer to the order in which the points were obtained in the simplex procedure. 109 only required making the assimilation of a new measurement conditional on its measured response. A three-dimensional representation of the response surface generated from the model is shown in Figure 5.4. It will be noted that the optimum response of the modeled surface is still short of the observed maximum of 610 pA, even with the rejection of several points, and that points far from the optimum are poorly modeled. Again, these problems arise from the fact that the surface is not really quadratic. Indeed, the optimum is so sharp that the original workers feel it would have been missed totally by the univariate search procedure they had originally envisaged [30,33]. Better modeling of the region near the optimum could be obtained if more points were rejected, since the higher responses exhibit more quadratic behavior. Nevertheless, the general shape of the modeled surface is the same as the experimental one and the interaction between the two acid concentrations is more clearly indicated in this representation. The equation for the modeled surface (unnormalized variables) is, z : 1402.9 [HCl] + 887.1 [HNOa] - 36.9 [H0112 (5.10) - 13.9 [HNOa]3 - 38.1 [HCl] [HNOa] - 14502.9 where z is the response in ”A and the concentrations are in mM. The interaction of the two variables is clearly indicated in the fifth term of this equation. The equation predicts an optimum response of 474 ”A when [HCl] 2 8.7 mM and [HNOa] = 20.0 mM. The largest response observed from the simplex data is 610 “A at [1101] : 8.5 mM and [HN03] = 23.2 mM. While the magnitudes of the optimum responses are 111 not in complete agreement due to the fact that a quadratic model was used, they appear at approximately the same location on the surface. D. Advantages of the Kalman Filter Approach While the results of the Kalman filter method and SVD are essentially the same, the Kalman filter exhibits certain practical advantages over SVD. First and foremost among these is the fact that the Kalman filter is a recursive algorithm. This means that as each new measurement is made, it can be readily assimilated into the existing model and a new model equation can be obtained. This dynamic nature allows the evolution of the surface to be observed as experiments are conducted and permits the modification of planned experiments if necessary. Convergence of the model can be more easily evaluated by examining the changes in the model parameters, covariance matrix, and predicted measurements as more and more measurements are processed. The omission of measurements from the fit is a trivial matter, simply requiring that they not be assimilated into the existing model. The Kalman filter approach also offers numerous advantages from a programming perspective. First, the algorithm described earlier is relatively simple and involves straightforward matrix manipulations. Singular value decomposition, on the other hand, involves a more complex algorithm which is often implemented blindly, and mistakes in the transcription of code are easily made. Problems with the Kalman filter may be more easily diagnosed for this reason. Program size requirements are also less with the Kalman filter method, both in terms of the amount of code and the size of arrays required. The reduction in the amount of code is due to the simplicity of the algorithm. The 112 Kalman filter programs typically required about one-third of the amount of code required by SVD. Array storage requirements are less because the Kalman filter is recursive and does not need to store all of the original data. The size of the arrays required by the Kalman filter increases only with the complexity of the model and not with the number of measurements. With SVD on the other hand, an array containing all of the basis functions for each data point needs to be maintained. Finally, there is a speed advantage with the Kalman filter. While the SVD algorithm is about five times faster than the Kalman filter algorithm from start to finish, the Kalman filter requires far less time to compute the new model when a new measurement is made (for example, about 0.75 s for the two-variable case on an IBM PC-XT with double-precision arithmetic, an 8 MHz clock, and no math co-processor). The SVD algorithm needs to be started from the beginning whenever a new point is to be added. Thus, the Kalman filter shows greater potential for interactive environments. The linear Kalman filter algorithm described here is very stable. It is not sensitive to the initial guesses for the state vector or covariance matrix, as long as the latter is sufficiently large. Successive iterations of the filter are not required, since it provides optimal estimates on the first pass. System observability [21] is not a problem due to the nature of the model. Perhaps the only disadvantage of the Kalman filter may be when an extremely large number of measurements are to be made. In this case, it is conceivable that round-off errors make the reliable assimilation of new measurements difficult. In such cases, however, any method is expected to have difficulties. Also, variations of the Kalman 113 filter using square-root algorithms [20,21] minimize this problem. Response surfaces in analytical chemistry should not experience this difficulty since the number of measurements is usually limited by practical considerations. As demonstrated here, more than 200 measurements can be easily accommodated with the Kalman filter method. . E. Relationship to Simplex Optimization Since its introduction [34,35] simplex optimization has been widely used in analytical chemistry and many variations of the basic method have evolved (see [8] and references therein). Generally speaking, simplex optimization is still the fastest and most reliable way to locate the response optimum on an unknown surface. However, simplex provides little direct information regarding the nature of the response surface and the algorithm must be restarted if the form of the response function is changed. Modeling, on the other hand, describes the shape of the response surface and can still be used if the form of the response function changes. While the optimum response can be determined through modeling approaches, as shown here, simplex is faster and more amenable to complex surfaces. Thus, simplex optimization and surface modeling are complementary methods. The advantages of their combination have been noted in the literature [8]. Smith et a1 [36] have recently described an algorithm for Optimization, called optiplex, that is based on polynomial fitting. While this method should have some utility for smooth surfaces with no variable interactions (no cross terms are used in the model), it is not as robust as the simplex method. Also, the use of traditional modeling algorithms could make the optiplex method cumbersome to use, especially 114 for complex models where interaction terms are included. The Kalman filter algorithm described here, however, is much more compatible with simplex methods because it is naturally recursive. Furthermore, it could be used in conjunction with simplex to generate complementary modeling information as an experiment progresses. Efforts in this area are anticipated. Such a combination would allow the continuous updates of the predicted optimum to be displayed (when possible) and surfaces to be analyzed while the optimization is proceeding. CONCLUSIONS For least squares modeling of response surfaces, the Kalman filter has been shown to produce results identical to those obtained with the more widely accepted singular value decomposition method. This was demonstrated with a second-order polynomial fit to a seven-variable system in which 36 coefficients were estimated. Further aspects of this type of modeling were illustrated with two additional chemical systems. The Kalman filter method was found to have a number of advantages over traditional approaches to response surface modeling, including its recursive implementation, speed, and simplicity. The recursive nature of the filter makes it ideal for real-time applications in areas such as simplex optimization. 100 ll. 12. l3. 14. 15. 16. 17. 18. 115 LIST OF REFERENCES G.E.P. Box, Biometrics, 10, 16-60 (1954). R.H. Myers, "Response Surface Methodology"; Allyn and Bacon, Boston, MA, 1971. R. Mead and DJ. Pike, Biometrics, 31, 803-851 (1975). G.E.P Box and N.R. Draper, "Empirical Model Building and Response Surfaces"; John Wiley and Sons, New York, 1987. S.L. Morgan and S.N. Deming, Anal. Chem., 46, 1170-1181 (1974). S.N. Deming and LR. Parker, Jr., CRC Crit. Rev. Anal. Chem., 7, 187-202 (1978). D. Betteridge, A.P. Wade and A.C. Howard, Talanta, 32, 709-722 (1985). D. Betteridge, A.P. Wade and A.G. Howard, Talents, 32, 723-734 (1985). K.M. Cellier and H.C.T. Stace, Appl. Spectrosc., 20, 26-33 (1966). S.L. Morgan and C.A. Jacques, J. Chrom. Sci., 16, 500-505 (1978). S.N. Deming and S.L. Morgan, Clin. Chem., 25, 840-855 (1979). R.J. Matthews, S.R. Goods and S.L. Morgan, Anal. Chim. Acta, 133, 169-182 (1981). S.R. Goode and K.W. Baughman, Spectrochim. Acta, 38B, 75-80 (1983). K.J. Johnson, "Numerical Methods in Chemistry"; Marcel Dekker New York, 1980; Chapter 4. A.C. Norris, "Computational Chemistry"; John Wiley and Sons, New York, 1981; Chapter 5. W.H. Press, B.P. Flannery, S.A. Teukolsky and W.T. Vetterling, "Numerical Recipes"; Cambridge University Press, New York, 1986; Chapter 14. P.D. Wentzell, A.P. Wade and S.R. Crouch, Anal. Chem., accepted for publication. R.E. Kalman, Trans. ASME, Ser. D: J. Basic Eng., 82(March), 35-45 (1960). 19. 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 35. 36. 116 A. Gelb (Ed.), "Applied Optimal Estimation"; MIT Press, Cambridge, MA, 1974. R.G. Brown, "Introduction to Random Signal Analysis and Kalman Filtering"; John Wiley and Sons, New York, 1983. S.D. Brown, Anal. Chim. Acta, 181, 1-26 (1986). W.S. Cooper, Rev. Sci. Instrum., 57, 2862-2869 (1986). S.C. Rutan, J. Chemomet., 1, 7-18 (1987). H.N.J. Poulisse, Anal. Chim. Acta, 112, 361-374 (1979). H.W. Sorenson, IEEE Spectrum, 7(July), 63-68 (1970). T.A. Brubaker, R. Tracy and C.L. Pomernacki, Anal. Chem., 50, 1017A-1024A (1978). C.L.M. Stults, A.P. Wade and S.R. Crouch, Anal. Chim. Acta, 192, 155-163 (1987). D. Betteridge, T.J. Sly, A.P. Wade and J.E.W. Tillman, Anal. Chem., 55, 1292-1299 (1983). A.P. Wade, Ph.D. Dissertation, University College Swansea, Swansea, West Glamorgan, SA2 8PP, United Kingdom, 1985. D. Betteridge, A.P. Wade, E.A. Neves and I. Gutz, An. Simp. Bras. Electroquim. Electrosnal., 3rd, 2, 411 (1982). N.R. Draper and H. Smith, "Applied Regression Analysis"; Wiley, New York, 1966; p. 118. D.E. Honigs, J.M. Freelin, G.M. Hieftje and TB. Hirschfeld, Appl. A.P. Wade, Anal. Proc., 20, 523-527 (1983). W. Spendley, G.R. Hext and F.R. Himsworth, Technomet., 4, 441-461 (1962). J.A. Nelder and R. Mead, Computer J., 7, 308-313 (1965). 8.8. Smith, Jr., M.A. Sainz and R.G. Schleicher, Spectrochim. Acta, 42B, 323-332 (1987). CHAPTER VI Recycle Heartcut Chromatography: Theory and Evaluation "When you can measure what you are speaking about, and express it in numbers, you know something about it. But when you cannot - your knowledge is of a meagre and unsatisfactory kind. " -Lord Kelvin In this chapter, a model is developed for the evaluation of recycle heartcut chromatography (RHC). While this technique has been used in our laboratory for some time, claims regarding its potential have been not been quantitatively established. With the aid of the model developed here, the true capabilities of RHC are critically assessed. Generally speaking, it is the goal of the chromatographer to obtain the maximum resolution in the minimum amount of time. Equation 6.1 is a well known relationship [1] which describes how chromatographic resolution is related to a number of important parameters. _ 11.1—V.- - .431. - 117...... Rs — 4 a [P1 (6.1) In this equation, R3 is the resolution, N is the number of theoretical plates, a is the separation selectivity, and k is the capacity factor. It is clear that resolution in chromatographic systems may be improved by increasing the selectivity, the capacity factor, or the number of 117 118 theoretical plates. While changes in the first two variables are generally the preferred means to effect a better separation, limitations are often reached to the adjustments which can be made. Separation selectivity can only be controlled to a certain degree by the chromatographer, who is limited by column characteristics, and changes in the capacity factor become less effective when It becomes large (>10). When enhancements to the selectivity and capacity factor have reached their practical limits, resolution can only be improved by increasing "the number of theoretical plates. One of the ways to increase the effective value of N is by recycling the partially resolved mixture through the same column. The recycling method for improving resolution was originally suggested for gas chromatography [2] and later extended to liquid chromatography [3]. Because of differences in flow dynamics and other factors, the two areas of application are quite distinct. This chapter will be devoted entirely to recycling in high-performance liquid chromatography (HPLC). The advantages of recycling methods over an extension in column length for HPLC have been discussed in the literature [4]. While recycling methods in HPLC have been known for some time, their application has been somewhat limited, possibly due to instrumental requirements. Recycle chromatography has traditionally been carried out by one of two general approaches, the closed loop (CL) method [4-10] and the alternate pumping (AP) method [4,11-13]. The instrumental configurations for these two methods are shown in Figure 6.1. The closed loop configuration consists of the pump, injector, column, and detector (D) connected in a continuous loop. This is 119 "UV-‘10" Column Closed Loop Method Pump Waste fl Solvent Reservoir _ Column 2 _0— Alternate 2 *1 Pumping Method v A @— Circulation Valve Method Waste ‘. Solvent Reservolr Figure 6.1. Instrumental configurations for recycling methods in HPLC. 120 conceptually the simplest of the two arrangements, but the major drawback of this approach is the band broadening which occurs when the partially separated sample passes through the pump chamber. This requires the internal volume of the pump to be small (<5 pl). Martin et al. have published a thorough theoretical and experimental study which focuses primarily on this approach [4]. The alternate pumping configuration utilizes two columns, two detectors, and a six-port valve. The sample is alternately switched between two columns until an adequate separation is achieved. The AP approach minimizes extra-column band broadening and does not require a special pump, but does require two identical columns and two high— pressure detectors. It has also been shown that the resolution attainable with the AP method is ultimately limited by the change in capacity factors which occurs when column pressure is increased [13]. A variation of the AP principle is used for boxcar chromatography [14], which attempts to simultaneously improve resolution and increase throughput. A variation of the closed loop method, which will be designated as the circulation valve (CV) method here, was introduced by Minarik et al. [15]. This is the third configuration shown in Figure 6.1. The configuration employed is essentially the same as that for the CL approach, but avoids the necessity of passing the sample through the pump by using a circulation valve driven by a stepper motor. The valve consists of a number of short channels arranged circumferentially on a disk which can rotate to connect adjacent pairs of channels to connecting ports. In this way, small segments of eluent may be continuously switched to the head of the column. This method avoids 121 band broadening introduced by the pump, but requires the use of the specialized valve. An alternate approach to recycle chromatography, recycle heartcut chromatography, was developed in our laboratory [16,17] and is considered in this chapter. The term "heartcut" arises from the fact that a only selected portion of the eluting components is recycled. The terminology is perhaps more appropriate when only a small section of the eluting peak profile is recycled, but its usage has persisted for the case where the cut contains more than one peak. The RHC method utilizes a heartcut 100p onto which part of the eluting sample is directed. After all of the components have passed through the detector, those in the heartcut loop are directed back to the head of the column for recycling. This process is continued until adequate resolution is achieved or the limit recycle number is reached. The principal advantage of this method over other approaches is the simplified instrumental requirements, which make it easily adaptable to existing chromatographic systems. Although RHC has been the subject of other studies [16,17], this is the first study of the theoretical and practical limitations of the method. A model is developed to investigate the resolution enhancements possible with this method and to simplify the complex interactions of variables involved. The model allows the potential and limitations of RHC to be evaluated. Recycle heartcut chromatography is compared with the CL and AP methods to assess its usefulness. 122 EXPERIMENTAL S ECTION A. HPLC Apparatus A block diagram of the apparatus employed for recycle heartcut chromatography is shown in Figure 6.2. The HPLC system is assembled from a Spectraphysics SP8700 ternary solvent delivery system (Santa Clara, CA), a Rheodyne 7120 injection valve (Cotsti, CA) with a 20 111 sample loop, a Dupont 4.6 mm x 25 cm Zorbax ODS column (Wilmington, DE), two Rheodyne valves (7010 and 7040) with pneumatic actuators, and a Chromatronix 220 UV detector. The heartcut loop consisted of 0.02 in. i.d. stainless steel tubing of various lengths and coiling diameters. The system was operated under microcomputer control [17-20]. The valve configurations for the system in the recycle heartcut mode are shown in Figure 6.3. Although experiments were conducted with a system constructed from two six-port valves, four-port valves may be used instead and this configuration has been hsown in Figure 6.3 for clarity. The start-up configuration is the normal arrangement for HPLC, with the sample passing through the column and detector to waste. In the heartcut configuration, the eluting components are directed onto the heartcut loop after passing through the detector. The components isolated in this way remain in the heartcut loop in the column flush configuration while the remaining components elute from the column. Finally, the components are directed back to the head of the column in the recycle heartcut configuration. The heartcut/flush/recycle is then repeated as necessary. 123 .hgqaaanaa—oaao espouse: £0.32 cacao—me: 3 pom: satay—amazon. o>~a> use 8393 0.3: .N.@ 0535 _:o>Lonom ocozom I go Loooo_c_ — \Gc A v hY1 _ aEam _ Looooooo _ - :Eazo 124 C E H Start-u " 011 Pump Waste C E Heartcut - Gil-j 4. ‘ t5 Pump Waste C E Column Flush , 011 Pump Waste C E Rec cle Heartcut " ’ M 01¢ Pump Waste Figure 6.3. Valve switching sequence used in recycle heartcut chromatography. 125 Note that this arrangement places no internal volume constraints on the pump and does not require an additional column or high pressure cells. Only the two four-port valves are required over the basic HPLC system. B. Separations The performance of RHC was evaluated with solutions of napthalene and biphenyl in methanol. Separations reported here employed isocratic elution with 92.5% methanol in water at a flow rate of 1 ml/min. Dispersion in the heartcut loop was evaluated from several recycles of a pure napthalene solution (400 pM). A computer program (CPACK.FTN) written by the author on a DEC LSI-11/23 minicomputer was used for the variance calculations. Separation enhancements with RHC were examined with a mixture of napthalene (250 11M) and biphenyl (200 11M). Spectral grade HPLC solvents (Burdick and Jackson, Muskegon, MI) were passed through a 0.45 pm filter and degassed with helium prior to use. Solutions were prepared from these solvents and also filtered prior to use. RESULTS AND DISCUSSION In this section, equations for resolution enhancement with recycling methods are developed and applied to recycle heartcut chromatography. An equation for the limit recycle number is derived 126 and shown to depend only on the length of the heartcut loop. Limitations to the length of the heartcut loop are considered. The maximum resolution enhancement attainable is given as a function of length. Finally, the advantages of RHC are examined from the perspective of improved selectivity. A. Resolution Improvement For an ideal system, the increase in the effective number of theoretical plates, N, is directly proportional to the number of passes through the column, 11. It is easily shown that, in this case, the ratio of the resolution after n passes to the initial resolution is given by, mm m gm 3 = Jr? (6.2) This relationship is approximately valid for the AP method where extra- column band broadening is small, but for other recycling systems the band broadening introduced by the recycling process needs to be considered. If the volume standard deviation of the peak after one pass through the column is 61 and the standard deviation of an infinitely narrow plug passing through the recycling process is 6:, then we can define, 8 : (fr/0'1 (6.3) For comparison purposes, it should be noted that Martin et al. [4] use the notation 1:82. Since the standard deviation of a peak is more 127 easily visualized than its variance, ll will be employed here rather than 1. Using this definition, a modified form of eq 6.2 can be derived for non-ideal systems [4]: ELLE)- - ’7 (6.4) 113(1) I. 11+ 11- This equation reduces to eq 6.2 as 8+0. In addition it shows that, regardless of the value of 8, an improvement in resolution can always be achieved if the number of passes through the column is large enough. In order for this resolution improvement to be realized, we must have n>83. For the CL method, the value of 6’: is fixed by the characteristics of the pump and the number of passes through the column is ultimately limited by the column void volume. For RHC, a more complex relationship exists since both the dispersion and number of recycles depend on the geometry of the heartcut loop. This relationship is now investigated. B. Limit Recycle Number For recycle heartcut chromatography, the limit recycle number will ultimately depend on the volume of the heartcut loop. As the number of passes through the column increases, the fraction of the heartcut loop occupied by the components being separated increases due to further broadening and separation. Eventually it is no longer possible to obtain a heartcut without losing significant portions of the partially separated peaks. While there is no fundamental reason why partial peaks cannot 128 be recycled, it can be shown through convolution models and other means [21] that resolution enhancements are degraded when this happens since there is an undesirable shift in the peak maxima. In the limit, for a very small cut, the resolution will be the same as it was on the first pass through the column. Because of this, it will be assumed that the last recycle occurs when the condition V11 2 m6 is met, where V: is the volume of the heartcut loop, m is a constant, and 0' is the volume standard deviation of the peak in the heartcut loop (we assume in the interest of simplicity that all of the overlapped peaks have the same standard deviation). If n is taken to be the limit recycle number, then prior to the nth pass through the column, the sample will have passed through the column and the heartcut loop n-l times. The ~standard deviation of a single peak prior to the nth pass through the column will be given by, a.-. : J(n-l)+(n—l)83 a. (6.5) Substituting this relationship into the constraining equation V11 2 01611-1 and rearranging, it is easily shown that the limit recycle number is, VE/m2 + 0‘12 + Jrz m in = (6.6) 0'12 + drz In order to proceed further it is necessary to examine the standard deviation introduced by the heartcut loop, 6:. Since laminar flow should prevail in the heartcut loop, we can describe the dispersion 129 through a slightly modified form of the well known relationship of Taylor for dispersion in open tubes [22-24]: 1‘13. 2 ,,. = g ”“2471” (6.7) In this equation, d'r is the volume variance in mlz, l is the length of the heartcut loop in cm, r is the internal radius of the loop in cm, F is the flow rate in ml-S'l, and D is diffusion coefficient of the solute in the mobile phase in cmzs'l. The variable 3, which will between 0 and 1, has been introduced to the original equation to accommodate reductions in dispersion brought about by other factors, such as coiling. For straight tubes 3 = 1, while for coiled tubes fractional values of R will be obtained. The effect of coiling is discussed in more detail in a later section. The results of eq 6.7 for a given set of parameters can be substituted into eq 6.6 to yield the limit recycle number for particular values of 6'1 and m. From this the resolution improvement may be computed from eq 6.4. For small heartcut volumes, it is not possible to obtain more than one pass through the column. As the volume becomes larger multiple passes are possible and resolution improvements can be realized. These improvements are mitigated by two opposing effects, however. For small values of 0'1 the value of mm is large, but so is the value of 8, and no resolution enhancement is observed. For large values of 61, min-*1 and resolution is again limited to that obtained on one pass through the column. The greatest improvement in resolution is observed for intermediate values of 6'1. 130 An important observation regarding eq 6.6 is that all of the terms except 61 are proportional to r‘. This means that any changes in 61 can, in theory, be compensated by changes in the radius of the heartcut loop so that n11- can be held constant. For example, if 61 doubles, r can be increased by a factor of J2 so that the result remains the same. In addition, 8 remains constant under these conditions, since 61 and 6: will change by the same factor. This implies that the theoretical maximum improvement in resolution which can be obtained for a given flow rate and diffusion coefficient depends only on the length of the heartcut loop. Furthermore, this maximum enhancement in resolution is observed to increase with length (Va increases with l, 61- increases with J1). The length of the heartcut loop cannot increase indefinitely, however, since it is limited by practical considerations. These are considered next. C. Length Limitations Beyond the logistical considerations arising from the need to use long sections of tubing for the heartcut loop, the length is limited by three other factors: dispersion, pressure drops, and the column dead volume. It has already been stated that increased dispersion resulting from a longer heartcut is counteracted by improved resolution, so this is not an important factor. The other two factors do warrant consideration, however. Pressure drops are important mainly from the point of view of the detector cell and, in some cases, the column being used. If RHC is to be used with a conventional chromatographic system, a high pressure detector cell cannot be assumed. In the heartcut configuration, the 131 heartcut loop is connected between the detector and waste, and the absolute pressure on the cell at this point may cause leaks to occur. The pressure drop across the heartcut loop is given by the Poiseuille equation as shown in eq 6.8. - filial - sr‘ (6.8) In this equation, AP is the pressure drop in dynes'cm'z, F is the flow rate in ml-seC'l, h is the viscosity in poise, and r is the radius in cm. For the purposes of practical calculations, we considered that the tubing for the heartcut loop had an internal diameter of 0.01 to 0.02 inch. The detector cell employed in this work is an older model and the maximum recommended pressure drop was 50-100 psi. Currently available detectors (non-high-pressure) have higher ratings, typically in the range 500-1000 psi [25]. The length of the column is limited to the point where the pressure drop across the heartcut loop is equal to the maximum which can be tolerated by the detector cell. This limitation is found to be most important for small diameter loops. The other factor limiting the length of the heartcut loop is the dead volume of column. After multiple passes through the column and the heartcut loop, the peaks will eventually broaden so that they occupy the dead volume of the column. When this happens, it is only possible to take one more cut, since further recycles would require starting one out before the last one has been flushed onto the column, and this is not a valid configuration. It would still be possible to recycle partial peaks, but it has already been stated that this is not a desirable option. Therefore, there is no advantage to increasing the volume of the 132 heartcut loop beyond the dead volume of the column as eq 6.6 becomes invalid at this point and there is no increase in the limit recycle number. Further increases in volume will increase the dispersion in the heartcut 100p and therefore degrade resolution enhancements which might otherwise have been achieved. D. Evaluation of Resolution Enhancements In order to evaluate the potential of RHC for improving resolution, calculations were carried out to determine the maximum resolution enhancement (Rs(l'l)/Rs(1)) possible as a function of the length of the heartcut loop. To do this, fixed values of F, D, m, and r were first assumed for a given length of heartcut loop. By employing equations 6.4, 6.6, and 6.7, a numerical procedure was then used to find the value of 61 which generated the maximum resolution enhancement. A numerical method of solution was employed because the order and complexity of the resulting equations precluded the development of an analytical expression. The enhancement values obtained in this manner are plotted in Figure 6.4. The calculations used F = 1 ml min“, D : 1x10'5 cm2 s'l, and m = 8. These values were considered reasonable to make a qualitative assessment of the method. Note that the value of r is not of critical importance to these calculations since, as mentioned earlier, it will affect only the value of 61 for which the optimum occurs, and not the optimum value itself. The maximum length shown in Figure 6.4 is determined by the point at which the limits set by the maximum pressure drop and the dead volume of the column meet one another so that no further increase in length is possible. Values used to determine this point were APmax = 500 psi and V0 = 3 ml, with 133 p=o.2 (1’) lllllLlLlllllllllll. Maximum Resolution Enhancement '0 fi = O .5 (CoHed) l3=1 1. / (Straight) 0. : 3 : 3 : 3 : 3 3 3 0. 10. 20. 30. 40. 50. Length of Heartcut Loop (m) Figure 6.4. Theoretical maximum resolution improvement attainable with RHC as a function of the length and geometry of the heartcut loop. 134 methanol as the solvent. Discontinuities in the figure are real and arise from the fact that only integral values of the limit recycle number can be used. The bottom curve in Figure 6.4 (3 = 1) gives the expected maximum resolution enhancement for a straight heartcut loop in which laminar flow prevails. In this case, it is clear that resolution improvements are marginal even if very long lengths of tubing are used. Furthermore, the actual resolution improvement observed may be less than the optimum, since the tubing radius is not continuously variable and a compromise would have to be reached. The bottom curve does not represent the most common situation, however, since coiling of the heartcut loop is essential and is known to reduce dispersion. The reduction in dispersion brought about by coiling is caused by a disruption of the normal laminar flow pattern in straight tubes and has been investigated by a number of workers [26-29]. We have found that by moderate coiling (0.02 inch i.d. stainless steel tubing on a ca. 0.5 inch diameter), the value of F in eq 6.7 can be reduced to 0.51 (napthalene in a 92.5% methanol/water mobile phase at 1 ml/min). This calculation assumed that all of the broadening was due to the flow and did not isolate contributions by other factors, such as tubing connections. The result is in almost perfect agreement with the value of 3 = 0.53 obtained by Selavka et al. [29] under similar conditions (0.02 inch i.d. Teflon tubing on a 0.8 inch diameter, with sodium benzoate in 50% methanol/water at 1 ml/min), and in reasonably good agreement with the value of 0.41 found by Scott and Simpson [28] under somewhat different conditions (0.01 inch i.d. stainless steel tubing on a 0.5 inch diameter, with sodium chloride in water at 1 ml/min). Employing the 135 value 3 = 0.5 in eq 6.7 and repeating the calculations gives the middle curve in Figure 6.4. In this case, greater resolution enhancements are observed for shorter lengths of tubing, suggesting that RHC is a more viable option. Resolution improvements are still 'not unbounded, however. In order to fully evaluate the potential of RHC for resolution improvements, calculations were also carried out with R = 0.2, and the results of these are shown as the top curve in Figure 6.4. The predicted resolution enhancements observed in this case are quite significant and show that RHC can become a very useful technique if dispersion in the heartcut 100p can be reduced still further. While a reduction of this degree has not yet been achieved experimentally and reductions by simple coiling are probably near their limit, other possibilities exist. For example, the use of an analog to the single bead string reactor (SBSR) employed in flow injection analysis [30] may reduce broadening without significantly diminishing volume or increasing pressure drops. Other studies have found that "knitted" [29,31,321 or geometrically deformed tubes [33] can also significantly reduce dispersion. Values of B = 0.22 to 0.24 have been experimentally achieved by other workers with knitted geometries [29] and, while these may be more difficult to implement with stainless steel, it demonstrates the technical possibilities. Such possibilities have not yet been investigated, but at least the model indicates the level of performance which must be met. The validity of the model has been verified only to the extent that no significant resolution improvement has been observed with relatively short (3.7 m) heartcut loops. Improvements in resolution will only be 136 realized with longer heartcut loops, which have not been available. Experiments to verify resolution improvements are now planned and await the arrival of longer lengths of tubing. E. Selectivity Enhancements While improvements in resolution with RHC are currently limited by the dispersion that occurs in the heartcut loop, the potential of the method can be evaluated from another perspective. Quite often in a chromatographic system, the ability to quantitate a given component is not limited by the resolution, which may actually be sufficient, but rather by the fact that an overlapping interferent may be present in large excess. In such a case, RHC may be used with a relatively short heartcut loop to improve the selectivity of the method for the component of interest. To do this, a cut is made which favors the component of interest and excludes the majority of the interfering component. The heartcut taken in this manner can then be recycled. Since the cut in this instance is small and dispersion will be minimal, the resolution on the second pass should approach that on the first pass, but the dominance of the interfering component will be removed. This principle is demonstrated in Figure 6.5, which shows a suboptimal separation of napthalene and biphenyl. Although the concentrations are approximately equal, the higher sensitivity of the detector for biphenyl causes this signal to predominate on the first pass. This makes it difficult to quantitate the napthalene even though the two peaks are substantially resolved. By performing a heartcut which selects the first portion of the two peaks, the contribution of the biphenyl peak is significantly reduced after just one recycle. To 137 .1 ,JL r r I r I I I I I— I 200. 300. 400. Time (sec) Figure 6.5. Selectivity enhancement in the separation of napthalene and biphenyl with RHC after two passes through the column. 138 ensure that the second peak was not simply an artifact of valve switching, other chromatograms were also obtained which favored a larger fraction of the biphenyl peak. It is apparent that integration of the napthalene peak can be performed more accurately after recycling. The process described is essentially the same as fraction collection followed by re-injection except that it is completely automated and retains some of the separation which was achieved on the first pass. The improvement in selectivity by RHC is analogous to that described by Cooks and Busch [34] for tandem mass spectrometry, since the heartcut acts as a "gate" to help remove chemical noise. This is done at the expense of random measurement noise, which will increase proportionately because we are removing some of the component of interest. Therefore detection limits will increase when recycling is used. As shown in Figure 6.5, however, the reduction in the original signal can be small. The viability of RHC to improve selectivity will depend on a number of factors, such as the extent of the interference, the magnitude of the S/N, and the separation of the components. In any case, a flexible and reproducible method of valve switching is needed to ensure good precision. Work in our laboratory has demonstrated that reliable results can be attained when the instrument is placed under computer control [16]. CONCLUSIONS The model developed in this chapter indicates that recycle heartcut chromatography, which has fewer instrumentation requirements 139 than other recycle methods, can be used to enhance chromatographic resolution. The method only becomes useful, however, if long heartcut loops are used and/or dispersion in the heartcut loop can be reduced by coiling or other means. The method may be more valuable for providing selectivity enhancements which are not strictly dependent on resolution, such as the removal of major interferents in the determination of minor components. For such applications, the need for a flexible system of automated control to provide good precision is evident. An additional advantage of the heartcut method is that it allows compounds with a greater retention to be removed from the column before the recycle phase is started. This is not true with other recycle methods where the potential exists for co-elution of highly retained peaks with the recycled peaks. 10. 110 12. 13. 14. 15. 16. 17. 18. 19. 140 LIST OF REFERENCES E.L. Johnson and R. Stevenson, "Basic Liquid Chromatography"; Varian Associates, Palo Alto, CA, 1978; p. 39. A.J.P. Martin, in "Gas Chromatography", V.J. Coates, H.J. Noebels and LS. Fagerson, Eds.; Academic Press, New York, 1958; p. 237. J. Porath and H. Bennich, Arch. Biochem. Biophys., Suppl. 1, 152-156 (1962) M. Martin, F. Verillon, C. Eon and G. Guiochon, J. Chromatogr., 125, 17-41 (1976). K.J. Bombaugh, W.A. Dark and R.F. Levangie, J. Chromatogr. Sci., 7, 42-47 (1969). S. Nakamura, S. Ishiguro, T. Yamada and S. Moriizumi, J. Chromatogr., 83, 279-288 (1973). B. Coq, G. Cretier, J.L. Rocca and J. Vialle, J. Liq. Chromatogr., 4, 237-249 (1981). H. Kalasz, Chromatographia, 20, 125-128 (1985). A. Hirose and D. Ishii, J. Chromatogr., 363, 391-393 (1986). D. Alba and R. Meira, J. Liq. Chromatogr., 9, 1141-1161 (1986). R.A. Henry, S.H. Byrne and D.R. Hudson, J. Chromatogr. Sci., 12, 197-199 (1974). P. Kucera and G. Manius, J. Chromatogr., 219, 1-12 (1981). Sj. van der Wal, Chromatographia, 22, 81-87 (1986). L.R. Snyder, J.W. Dolan and Sj. van der Wal, J. Chromatogr., 203, 3-17 (1981). M. Minarik, M. Popl and J. Mostecky, J. Chromatogr. Sci., 19, 250-252 (1981). S.K. Ratanathanawongs, Ph.D. Thesis, Michigan State University, East Lansing, MI, 1986. P.M. Wiegand, Ph.D. Thesis, Michigan State University, East Lansing, MI, 1985. M.D. Wiegand, P.M. Wiegand and S.R. Crouch, Anal. Chim. Acta, 187, 241-253 (1986). B. Newcome and C.G. Enke, Rev. Sci. Instrum., 55, 2017—2022 (1984). 20. 21. 22. 23. 24. 25. 26. 27. 28. 29. 30. 31. 32. 33. 34. 141 S.V. Mediln, Research in progress, Michigan State University, East Lansing, MI. P.S. Williams, Dept. of Chemistry, University of Utah, Unpublished communication. G. Taylor, Proc. R. Soc. Ser. A, 219, 186 (1953). J. Ruzicka and E.H. Hansen, "Flow Injection Analysis"; John Wiley and Sons, New York, 1981. G. Horvai, CRC Crit. Rev. Anal. Chem., 17, 231-264 (1987). Spectraphysics (Santa Clara, CA) technical information, unpublished. K. Hofmann and I. Halasz, J. Chromatogr., 173, 211-228. R. Tijssen, Anal. Chim. Acta, 114, 71-89 (1980). R.P.W. Scott and C.F. Simpson, J. Chromatogr. Sci., 20, 62-66 (1982). G.M. Selavka, K.S. Jiao and LS. Krull, Anal. Chem., 59, 2221-2224 (1987). J.M. Reijn, W.E. van der Linden and H. Poppe, Anal. Chim. Acta, 123, 229-237 (1981). M. Uihlein and E. Schwab, Chromatographia, 15, 140-146 (1982). H. Engelhardt and U.D. Neue, Chromatographia, 15, 403-408 (1982). K. Hofmann and I. Halasz, J. Chromatogr., 199, 3-22 (1980). R.G. Cooks and K.L. Busch, J. Chem. Ed., 59, 926-933 (1982). CHAPTER VII Conclusions and Future Perspectives "The Butcher would gladly have talked till next day, But he felt the lesson must end, And he wept with delight in attempting to say He considered the Beaver his friend. While the Beaver confessed, with affectionate looks More eloquent even than tears, It had learned more in ten minutes than all the books Would have taught it in seventy years. " - Lewis Carroll "The Hunting of the Snark" The preceding chapters have attempted to illustrate the importance of modeling methods in analytical chemistry through a variety of applications. In this chapter, the major conclusions of Chapters II through VI are summarized and ideas for future work are presented. In Chapter II, a new reaction-rate method of analysis, the two-rate method, was introduced. It was shown that this method can help solve the problem of between-run variations in the rate constant for systems following first-order kinetics. A model for the propagation of random errors in the rate measurements with the method was developed and shown to be valid with experimental data. This model indicates that the precision of the two-rate method is comparable to that of the optimized derivative method under selected measurement conditions. Further applications of the two-rate method to other first-order systems, including enzyme reactions, are anticipated. 142 143 Application of the method in commercial laboratories is also being investigated. Se‘veral reaction-rate methods of analysis for first-order systems were compared in Chapter III. This comparison is unique for two reasons. First, it considered "optimized" methods as well as traditional methods. Second, it considered the problem of between-run variations in the rate constant in the comparison. Models were developed to predict the significance of this problem for each method. As expected, the optimized methods show less susceptibility to between-run variations in the rate constant, but were generally less precise in the absence of those variations. Although the choice of a reaction-rate method depends on many factors, the best overall performance was exhibited by the Cornell method. Remarkably, this method has been virtually ignored in the literature. Further work needs to be done to examine its potential and practical limitations (i.e. its dependence on noise, number of data points, etc.). The application of the Kalman filter to reaction-rate methods of analysis was investigated in Chapter IV. The extended Kalman filter was found to be an effective means to minimize the effect of between-run variations in the first-order rate constant. Further applications of this method should be considered in the future. The linear Kalman filter, through its application to simulated data, was also shown to have great potential in the area of differential kinetic methods of analysis. The results of this method are equivalent to those obtained with linear least squares, which is the optimal approach. The advantages of the Kalman filter are that it is recursive and easy to implement. The method still 144 needs to be applied to experimental data to demonstrate its effectiveness, and experiments are planned for this purpose. Chapter V considered another application of the Kalman filter to analytical chemistry in the area of response surface modeling. The Kalman filter was shown to produce results identical to those of more traditional linear least squares modeling methods, but has a number of advantages arising from its recursive nature. These include speed, simplicity, adaptability, and the ability to observe the evolution of the model as new experimental points are obtained. The method was applied to three chemical systems to demonstrate its utility. Future plans are to utilize the recursive nature of the Kalman filter in conjunction with simplex optimization methods. Finally, in Chapter VI a model was developed to examine resolution enhancements attainable with recycle heartcut chromatography (RHC). This study is significant because no previous investigation into the potential of RHC had been performed. The model shows that RHC is only useful for enhancing resolution if long heartcut loops are used and dispersion in the heartcut loop can be significantly reduced. It was also demonstrated, however, that RHC can provide significant enhancements to selectivity where adequate resolution is not the real problem. Experimental verification of the model developed is to be carried out in the future. Also, techniques for the reduction of dispersion in the heartcut loop need to be investigated to ascertain whether or not RHC can be a useful tool for resolution enhancement. APPENDIX I APPENDIX I XYFIT - A Program for Nonlinear Curve Fitting This appendix describes XYFIT, a nonlinear curve fitting program based on a simplex algorithm. This program was written by the author in FORTRAN-77 and designed to run on DEC PDP-ll (or LSI—11) computers under the RSX-llM operating system, although it is not restricted to this environment. It has been used for the nonlinear fitting applications described in Chapters 11, III, and IV. The intent behind the development of XYFIT was to create a simple but useful nonlinear curve fitting program. It was not intended to replace more fully evolved software packages such as KINFIT [1] or MINPACK [2,3], but rather to act as alternative with PDP-ll operating environments. Historically, XYFIT is the second simplex-based fitting program written in our laboratory. The first, PFIT06, was written in FORTRAN-IV by Wai Tak Law for use on a PDP—8 computer [4]. This program became outdated when the PUP-8 was replaced with an LSI-11/23 as the principal laboratory computer. At the time XYFIT was developed, KINFIT was being adapted from a CDC Cyber to 3 DEC VAX environment. At that time, there existed no general nonlinear curve fitting package for use under the RSX operating system. This posed a problem, since it required the time consuming task of transporting data to other systems. This motivated the development of XYFIT. While some of the ideas for 145 146 XYFIT arose from PFIT06 [4] and other sources [5,6], the program is a unique entity which has evolved from specific design goals. Shortly after XYFIT was implemented, the MINPACK package [2,3], based on the Levenburg-Marquardt algorithm [3,7,8], was adapted to the RSX operating system by Mark Phillips of this department [9]. This still posed problems, however, since the routines were not easily modified for specialized applications and required the definition of the Jacobian matrix [10] of derivatives. The latter task was tedious when multiple models were to be tested and prone to errors with complex functions. For these reasons, the MINPACK package has not been extensively used in our research group. The primary objective in the creation of XYFIT was simplicity. This meant that the code should be easy to follow and modify, and the program should be easy to use. This goal was achieved in a number of ways. First, the modified simplex algorithm of Nelder and Mead [11] was used to obtain the best fit by minimizing the sum of the squares of the residuals (SSR). This algorithm is relatively straightforward and does not require the definition of the Jacobian matrix. Second, the original version of the program assumed only one independent and one dependent variable (the word "original" is used since the has been modified by several users). Third, uniform weighting was assumed in the original version, requiring no variance estimates to be included with the data. While this limited the rigor of the algorithm, it was found to be a useful simplification for the majority of applications, where equal variances would be assumed in any case. Finally, since a reliable algorithm for the determination of convergence with the simplex method is difficult to implement, the burden of determining convergence was left 147 to the user. To aid in this determination, the program requests the number of iterations from the user and then outputs the SSR at ten. equally spaced intervals during the fitting process. To further test convergence, the user may choose to restart the fitting procedure with different initial estimates and/or more iterations. The program XYFIT is also written to be compatible with the general purpose plotting package, MULPLT, written in this department by Tom Atkinson. It reads data files which are MULPLT compatible and produces a MULPLT compatible output file which contains calculated values of the dependent variable and the residuals. This allows easy visual inspection of the quality of the fit. To use the program, the user first modifies the assignment of YCALC, the calculated value of the dependent variable, in the four line subroutine YFUNC. The value of YCALC is defined in terms of the independent variable, X, and several parameters, PAR(1)...PAR(n), that are being fitted. After the modification is made, the program is then compiled and executed. The user must supply the name of the MULPLT compatible file containing the data (a maximum of 500 points are allowed in the original version) and the number of parameters to be fit. For each parameter, the user is asked for an initial estimate and an initial step size. These values are used by the subroutine VERTEX in establishing the initial simplex. This subroutine creates an "equilateral" simplex centered around the initial estimate in parameter space. "Equilateral" is in quotation marks because the dimensions of the simplex are actually scaled according to the step size given, so it is not a true equilateral simplex. The algorithm used to determine the initial simplex is adapted from the recursive method used by Long [6], 148 although he does not explicitly describe the mathematical details. These may be found in the program listing which follows. (The routine used in PFIT06 [4] to establish the initial simplex was not used because it did not generate an equilateral simplex and could conceivably result in premature contraction of the simplex, especially for higher dimensions. After the initial simplex is established, the user provides the desired number of iterations and the simplex optimization proceeds. Normally 10 to 100 iterations are used, but this depends on the number of data points and other factors. The modified simplex algorithm [11], which is limited to simple reflection, expansion, and contraction, is used in XYFIT. Law [4] had employed the super modified simplex algorithm [12], but some of the initial tests indicated premature contraction with this method. Other modifications of the simplex method [13] were not used as it was felt that they detracted from the simplicity of the program. The program reports the SSR as the minimization proceeds and when finished reports the best parameter estimates obtained and the corresponding SSR. The user may then repeat the fitting procedure with new initial estimates, or generate a MULPLT compatible data file with the calculated points and corresponding residuals. Because of the simplicity of XYFIT, it has been easy to implement and easy to modify by users. The first conclusion is supported by the number of individuals who have employed it for a variety of purposes. In addition to the applications by the author, it has been used to fit reaction curves [14,15], rate laws for chemical kinetics [15], fluorescenae/phosphorescence lifetime data [16-19], models for absorbance-corrected fluorescence [20], laser-enhanced ionization signal- 149 to-noise ratio expressions [21], and models for dual laser ionization [22]. The ease of modification is demonstrated by numerous changes which have been made by individual users to suit specific needs. Modifications have included the addition of "windowing" to select the range of data to be used [16,20], extension to multiple independent variables [15], addition of reciprocal variance weighting [20], and adaptation to an IBM PC environment [19]. Despite its success, the program XYFIT has certain disadvantages. First, since simplex is an exploratory method rather than a "steepest descent" approach, convergence can take more time. This is normally insignificant, however, compared to the time required to transport data to another machine. The major drawback of the simplex algorithm is that it does not provide estimates of the error in the parameters, as do KINFIT and MINPACK. Such estimates are difficult to generate in the absence of a Jacobian matrix. Numerical methods can be used to obtain the covariance matrix [10] and this was attempted. The attempt was unsuccessful, however, and complicated the program, defeating the original objectives. To summarize, the program XYFIT has proven to be a successful nonlinear curve fitting program for relatively simple applications. A listing of the program and its associated subroutines follows. 150 PROGRAM XYFIT SIMPLE! fitting pecksge Written by : Peter Wentzell Dept. of Chesistry Michigan State University 5 June, 1985 Copyright Reptilisn Lifeforss Inc. References: l) I. Spendley et al., Technosetrics, 4, 441 (1962). 2) J.A. Nelder and R. Mead, Conputer J., 7, 308 (1965). 3) S.N. Dening end S.L. Morgen, Anal. Chem., 45, 278A (1973). 4) M.H. Routh et al., Anal. Chen., 49, 1422 (1977). 5) S.R. Dening end L.R. PArker, CRC Crit. Rev. Anal. Ches.,7,187 (1978). 6) P.P.A. van der Niel, Anal. Chin. Acta, 122, 421 (1980). 7) 0.1. Creik et al., J. Chen. Inf. Conp. Sci., 23, 31 (1983). This progren is designed to fit sinple functions by the Sinplex nethod. Date to be fit must consist of only one independent variable (is x,y pairs) and the nexisun nunber of points which can be fit is 500. Up to five persseters say be fit with this progras. To use this program the user nust provide the following. First, the subroutine YFUNC sust be provided or sodified to calculate a value for y (YCALC) for the function to be fit. The progren can then be conpiled and run. The user will be asked for the nane of the data file containing the date to be fit. This data needs to be in MULPLT forest. After reading in the date, the progrss asks the user to provide the nusber of psrsseters, initial estinetes for each paraseter, and a step size for each parsneter. The step size is used to detersine the dinensions of the initial sisplex. It is, to ease extent, a seasure of confidence in the initial estisete for that psreseter. Finally, the user is asked to provide the nusber of iterations to be carried out. This places the burden of detersining convergence on the user. This was done because no reliable and consistent scans for testing for convergence could be found which was also relatively sisple. To aid in the detersinstion of convergence, the sun of the squares of the residuals is output at intervels during the fit. After the fit has been carried out, the user has the option of outputting the calculated curve to a MULPLT file. This file is in the forsat RD x ycalc resid , where ’ycslc' is the calculated y value and ’resid’ is the residual at that point. The date are written with e 3615.6 fornet and the third colunn is included in case e plot of the residuals is desired. finally, the fitting process say be carried out again in order to test for true siniss. The subroutines called by this progrss are : VERTEX - establishes initial sinplex (its vertices) based on initial estisetes of the pereneters and step sizes given by the user. OGO0000000GOO00OO0000O0000000000000OOOOOOOOOOOOOOOOOO O 151 C YEUNC - coaputes a calculated y (YCALC) based on the values of x and the C paraaeters passed to it. The user aust provide the function in C this subroutine. C RPUNC - conputes the sun of squares of the residuals (response function) C for a given set of paraaeters. C SIMPLE - conputes the new point in the sinplex based on the Modified C Siaplex Method, which allows expansion and contraction. C c33888838388=8=88888==83883338888338:888338888888888.88888833388.88888832883238 C C Set up house C DIMENSION X(500),Y(500),PAR(5,4),SSIZE(4),SSR(5),PTMP(4) CRARACTERt32 PNAME BYTE LINE(80),IASE C C Now get soae info froa the dope stooped over the keyboard C TYPE 3, ’SIMPLEX FITTING PROGRAM’ TYPE 3, ' ' TYPE 100 100 FORMAT (’ NAME OF INPUT PILE?’,$) ACCEPT 110,PNAME llO FORMAT (A32) OPEN (UNIT81,PILE=ENAME,STATUS=’OLU’) TYPE t,’READING DATA.’ C Read in file with MULPLT for-at. Ignore conaents. Count as we go. 181 120 CONTINUE READ (1,130,END=200,ERR=120) LEN,LINE 130 FORMAT (0,80Al) IE (LINE(1).LT.’A’) GOTO 120 DECODE (LEN,140,LINE,ERR=120) X(I),Y(I) 140 FORMAT (2X,2615.6) I=I+l GOTO 120 200 CONTINUE NPTS=I-l TYPE t,’FINISNED READING FILE. NUMBER OF POINTS READ =’,NPTS CLOSE (UNIT=1) C C Got the data. Now get soae other info. Subroutine VERTEX will take care C of getting the initial paraneter estiaates. C 205 TYPE 210 210 FORMAT (’ NUMBER OF PARAMETERS (l-4)?’,$) ACCEPT S.NPAR IP ((NPAR.LT.1).OR.(NPAR.GT.4)) GOTO 205 213 CONTINUE CALL VERTEX(PAR,SSR,NPAR,X,Y,NPTS) 215 TYPE 220 220 FORMAT (’ONUMEER OE ITERATIONS?’,$) ACCEPT ¥,NMAX 152 IF (NMAX.LT.1) GOTO 215 NINT!INT((FLOAT(NMAX)+O.5)/10.0) C C All set. Now do siaplex optiaization. Print results at each tenth of the C way there. C C Shit on string - fan enable C TYPE t,’SIMPLEX OPTIMIZATION PROCEEDING.’ NOUT=NINT ICNT=O 300 CONTINUE CALL SIMPLX(PAR,SSR,NPAR,X,Y,NPTS,ICNT,ID) ICNTSICNT+1 IF (ICNT.EO.NMAX) GOTO 400 IF (ICNT.EO.NOUT) TEEN TYPE t,’AFTER ’,ICNT,’ ITERATIONS SSR = ’,SSR(ID) NOUT=NOUT+NINT END IF GOTO 300 C C All done with fit - print results C 400 CONTINUE TYPE t,’ ’ TYPE t,’FINISRED.’ TYPE t,’ ' TYPE t,’SUM OF SQUARES OF RESIDUALS = ’,SSR(IE) TYPE t,'DEST FIT PARAMETERS :' DO 410 I=l,NPAR TYPE 420, I,PAR(IB,I) 420 FORMAT(' PARAMETER ’,12,’ = ’,Gl$.6) 410 CONTINUE C C Now check to see if a MULPLT file is to be generated C TYPE 430 430 FORMAT (’ODO YOU RANT A MULPLT FILE OF THE CALCULATED FIT (Y/N)?’,$) ACCEPT 440, IASR 440 FORMAT (Al) IF (IASE.NE.’Y') GOTO 500 C C Generate MULPLT file. Calculated fit has RD tags, residuals are in third C coluan of file. C TYPE 445 445 FORMAT (’ MULPLT FILE NAME?’,$) ACCEPT 110, FNAME OPEN (UNIT=1,FILE=FNAME,STATUS='NEW’,CARRIAGECONTROL=’LIST’) TYPE t,’NORRING.' DO 450 I=l,NPAR PTMP(I)=PAR(IE,I) 450 CONTINUE DO 460 I=l,NPTS CALL YFUNC(X(I),PTMP,YCALC) 153 RESIDSYCALC- -Y(I) NRITE (l, 470) X(I),YCALC, RESID 470 FORMAT (’RD’, 3G15. 6) 460 CONTINUE CLOSE (UNITII) TYPE t,’MULPLT FILE NR1TTEN.’ C C Play it again, San? C 500 CONTINUE TYPE 510 510 FORMAT (’ODO YOU WANT TO TRY AGAIN (Y/N)?',$) ACCEPT 44D, IASE IF (IASE.EO.'Y') GOTO 213 C C All done - go hone C TYPE 4,’IT”S EEEN REAL.’ STOP END 154 SUIROUTINE SIMPLE(Y,R,NPAR,X,Y,NPTS,ICNT,IE) Subroutine to perfora the siaplex fitting of data. The variables in the paraaeter list are : V - is the two-diaensional array containing the vertices (ie the paraaeters) of the current siaplex. R - is the one-diaensional array containing the response functions (sun of squares of residuals) for each vertex of the siaplex (each set of paraaeters). NPAR - is the nuabar of paraaeters being fit. The nuaber of vertices is equal to the nunber of paraaeters plus one. I - is the one-diaensional array containing values of the independent variable. Y - is the one-diaensional array containing the observed values of the dependent variable. NPTS - is the nusber of points being fit. ICNT - is the value of the iteration count IR - is the index of the vertex (set of paraaeters) giving the best fit (lowest sun of squares of residuals) on return froa the subroutine. 000000000000000000000 00000 C This subroutine uses the Modified Siaplex Method. Nhen passed a set of C siaplex vertices it will generate a new vertex to replace the worst one C according to the aodified siaplex algoritha. After this new point is C generated, the best point is deterained and its index is returned in ID. C The sun of the squares of the residuals is used as the criteria for C the best response. This is evaluated with the subroutine RFUNC. In the beginning.... DIMENSION V(5,4),R(5),X(500),Y(500) DIMENSION PC(4),PR(4),PE(4),PCR(4),PCN(4),DELTA(4) C C Check to see if first tiae through (ICNT=O) and if so reset LASTN to C zero to indicate no new points found yet. Find best, worst, and next-to- C worst vertex. If worst point is the saae as last worst point, replace C worst point with next-to-worst point (huh?). C IF (ICNT.E0.D) LASTN=D IN=1 IE=1 INN=1 DO 100 I=l,(NPAR+l) IF (R(I).GT.R(IM)) IN=I IF (R(I).LT.R(IE)) 13:1 100 CONTINUE 155 IF (INN.EO.IN) INN=2 DO 110 I=l,(NPAR+l) IF ((R(I).GT.R(INN)).AND.(I.NE.IN)) INN=I llO CONTINUE IF (LASTN.EO.IN) IN=INN LASTNSIN C C Find centroid and reflected point, then find response at reflected point C DO 120 1:1,NPAR PC(I)80.O DO 130 J=l,(NPAR+l) IF (J.NE.IN) PC(I)=PC(I)+V(J,I) l30 CONTINUE PC(I)=PC(I)/FLOAT(NPAR) DELTA(I)=PC(I)-Y(IN,I) PR(I)=PC(I)+DELTA(I) 120 CONTINUE CALL RFUNC(PR,X,Y,NPTS,RPR) C C If reflected point is between worst and next-to-worst, replace worst C point with reflected point. C 200 CONTINUE IF ((RPR.GE.R(IE)).AND.(RPR.LE.R(INN))) TEEN DO 205 I=l,NPAR V(IN,I)=PR(I) 205 CONTINUE R(IN)=RPR C C If reflected point is better than best, try expansion and if expanded C point is better than best, retain it as new vertex. Otherwise, retain C reflected point. ELSE IF (RPR.LE.R(IE)) TNEN DO 210 I=l,NPAR PE(I)=PR(I)+DELTA(I) 210 CONTINUE CALL RFUNC(PE,X,Y,NPTS,RPE) IF (RPE.LT.R(IE)) TEEN DO 220 I=l,NPAR V(IN,I)=PE(I) 220 CONTINUE R(IN)=RPE ELSE DO 225 I=l,NPAR V(IN,I)=PR(I) 225 CONTINUE R(IN)=RPR END IF Point is worse than next-to-worst. If it is worse than worst point, contract in direction of worst point. Otherwise, contract in direction of reflected point. Replace worst vertex with appropriate contraction. 00000 156 ELSE IF (RPR.GT.R(IN)) TEEN DO 230 I=l,NPAR PCN(I)8PC(I)-O.53DELTA(I) Y(IN,I)=PCN(I) 230 CONTINUE ,CALL RFUNC(PCN,X,Y,NPTS,RPCN) 'R(IN)=RPCN ELSE DO 240 I=l,NPAR PCR(I)8PC(I)+O.53DELTA(I) Y(IN,I)8PCR(I) 240 CONTINUE CALL RFUNC(PCR,X,Y,NPTS,RPCR) R(IN)8RPCR END IF END IF C C New point has been found. Now look through to find best point. C DO 300 I=1,(NPAR+l) IF (R(I).LT.R(IE)) 13:1 300 CONTINUE C C That’s all folks 2!! C RETURN 0000000000000000 00000 110 00000000000000000 00000 157 SUDROUTINE RFUNC(PAR,X,Y,NPTS,SSR) This subroutine evaluates the sun of squares of residuals between the user supplied function and the observed data. The variables in the paraaeter list are: PAR - is the one-dinensional array containing the paraaeters being fit. I - is the one-diaensional array containing values of the independent variable. Y - is the one-diaensional array containing observed values for the dependent variable. NPTS - is the nunber of points being fit SSR - is the sun of the squares of the residuals returned by the prograa This subroutine requires the user supplied subroutine YFUNC to deteraine calculated y values. DIMENSION PAR(4),X(500),Y(500) SSR=0.0 DO 110 I=l,NPTS CALL YFUNC(X(I),PAR,YCALC) SSR=SSR+(Y(I)-YCALC)I(Y(I)-YCALC) CONTINUE RETURN END SUEROUTINE YFUNC(X,PAR,YCALC) This subroutine coaputes a value for y according to the function being fit and the paraaeters being fit. This subroutine must be supplied by or aodified by the user. The variables in the paraaeter list are: X - the value of the independent variable (a single value, not an array) for which YCALC is to be conputed. PAR - the one diaensional array containing the values of the paraaeters to be used in the calculation. YCALC - is the calculated value of y returned froa this subroutine. To aodify this subroutine for a certain function, the user siaply equates YCALC to the function in terns of X and the paraaeters in PAR. For exaaple, a sinple linear fit with the intercept as paraaeter l and the slope as paraaeter 2 would be represented as : YCALC = PAR(l) + PAR(2)8X DIMENSION PAR(4) YCALC=PAR(1)3(l.O-EXP(-X#PAR(2)))+PAR(3) RETURN END 0000000000000000000000000000000000000000000000000 00000 158 SUDROUTINE YERTEX(Y,R,NPAR,X,Y,NPTS) This subroutine sets up the initial vertices of the siaplex used by the siaplex optiaization subroutine SIMPLX in conjunction with the siaplex based curve fitting prograa XYFIT. The user is required to input the nuaber of paraaeters for the fit (i.e. nunber of diaensions) and initial values and step sizes for each paraaeter. The paraaeters returned froa this subroutine are described below. Y - is the array containing up to 5 vertices of 4 diaensions (paraaeters) each. R - is the one-dinensional array containing the response function at each vertex. These are calculated by this subroutine each tine a new siaplex is for-ed. NPAR - is the nuaber of paraaeters being fit (i.e. diaensionality of the siaplex space). X - is the one diaensional array containing values for the independent variable Y - is the one dilensional array containing observed values for the dependent variable being fit. NPTS - is the nuaber of points being fit ICNT - is a counter which keeps track of how aany tines the sun of squares of the residuals is evaluated. For the calculation of the vertices, this subroutine uses a recursive aethod eaployed by Long (see Anal. Chin. Acta, 46, 193 (1969)) although not explicitly described by his. This aethod was chosen over that eaployed by Wai Tak Law (see thesis) since the latter does not result in an equilateral siaplex. With the recursive aethod, the Nth vertex requires N-l coefficients to be aultiplied by the first N-l step values in order to obtain the increaents to be added to the first N-l parameters in the initial vertex. The coefficients for the initial (starting) vertex are, of course, all zeros. For the second vertex, the first coefficient is unity and the following ones are all zeros. For subsequent vertices, the coefficients are cal- culated recursively. Each new vertex requires calculating the N-l and N-2 coefficients (for the Nth vertex). If we let 0(N) represent the Nth coefficient, then this aeans that for the Nth vertex we need to calculate 0(N-1) and 0(N-2), with all of the preceding coefficients reaaining the sane and all of the following ones set to zero. Furthernore, if we let 2 equal the previous value of 0(N-1) (now O(N-2)) then it is easily shown that the new coefficients are, 0(N-2) s z - l/(ZtZ) and, 0(Nel) = SORT(Ztt2 - 0(N-2)8¥2) so that 00000000000000 159 Use of these coefficients will result in an equilateral siaplex with the distance between all of its points equal to unity if all of the step sizes are unity. After this siaplex has been set up using the aethod of Long, it is noved the first point (initial estiaates of paraaeters) is at the center of the siaplex. This reaoves any probleas which night be caused by having the siaplex always to one side of the initial estiaates. The response function (sun of squares of residuals) at each vertex is then calculated using the subroutine RFUNC. Initialization section DIMENSION V(5,4),R(5),SSIZE(4),X(500),Y(500),O(4),VTMP(4) O(l)81.O DO 0000 DO 120 130 110 C 90 182.4 O(I)=0.0 CONTINUE Input paraaeters for vertex calculation froa user 110 I=l,NPAR TYPE 3,’ ' TYPE l20,I FORMAT (’ ENTER INITIAL VALUE FOR PARAMETER',12) ACCEPT t, V(1,I) TYPE 130,1 FORMAT (' ENTER STEP SIZE FOR PARAMETER’,12) ACCEPT t,SSIZE(I) CONTINUE C Now calculate initial vertices. C TYPE 8,’COMPUTING INITIAL SIMPLEX.’ V(2,l)=V(l,1)+SSIZE(l) DO 140 D0 160 150 C C Vertices calculated - now offset so that center of siaplex is at starting 140 I=2,NPAR V(2,I)=V(l,I) CONTINUE 150 I=3,(NPAR+1) 2:0(1-2) 0(I-2)=Z-O.5/Z 0(I-l)=SORT(Z#Z-O(I-2)#O(l-2)) D0 160 J=1,NPAR V(I,J)=V(l,J)+O(J)¥SSIZE(J) CONTINUE CONTINUE C point input by user. C DO I70 I=l,NPAR CENTER=0.0 DO 180 J=l,(NPAR+l) CENTER=CENTER+V(J,I) 160 180 CONTINUE CENTERSCENTER/FLOAT(NPARtl) OFFSET8V(1,I) DO 190 J81,(NPAR+1) Y(J,I)8Y(J,I)-CENTER+OFFSET 190 CONTINUE 170 CONTINUE C C Siaplex calculated and centered - now deternine response function at each C vertex. C DO 200 I=l,(NPAR+l) DO 210 J81,NPAR VTMP(J)=V(I,J) 210 ' CONTINUE CALL RFUNC(VTMP,E,Y,NPTS,RTMP) R(I)=RTMP 200 CONTINUE RETURN END 10. 11. 12. 13. 14. 15. 16. 17. 18. 19. 20. 21. 161 LIST or REFERENCES J.L. Dye and V.A. Nicely, J. Chem. Ed., 48, 443-448 (1971). J.J. More, "NTIS Publication ANL 80-74"; Argonne National Laboratories, Argonne, IL, 1980. J.J. More, in "Lecture Notes in Mathematics 630: Numerical Analysis", G.A. Watson, Ed.; Springer-Verlag, New York, 1977; p 105. W.T. Law, Ph.D. Thesis, Michigan State University, East Lansing, MI, 1978. D.J. Craik, A. Kumar and 6.0. Levy, J. Chem. Inf. Comput. Sci., 23, 30-38 (1983). D.E. Long, Anal. Chim. Acta, 46, 193—206 (1969). C. Levenburg, Quart. App]. Math., 2, 164 (1944). D.W. Marquardt, SIAM J. Appl. Math., 11, 431 (1963). G.M. Phillips, Ph.D. Thesis, Michigan State University, East Lansing, MI, 1985. W.E. Wentworth, J. Chem. Ed., 42, 96, 162 (1965). J.A. Nelder and R. Mead, Comput. J., 7, 308 (1965). M.W. Routh, P.A. Swartz and M.B. Danton, Anal. Chem., 49, 1422-1428 (1977). D. Betteridge, A.P. Wade and A.G. Howard, Talents, 32, 723-734 (1985). R.G. Harfman, Ph.D. Thesis in writing, Michigan State University. M. Lopez-Nieves, Research in progress, Michigan State University. F.R. Kraus, Ph.D. Thesis in writing, Michigan State University. N. Herron, Ph.D. Thesis, Michigan State University, East Lansing, MI, 1987. H.D. Kim, Research in progress, Michigan State University. M. Newsham, Research in progress, Michigan State University. T.P. Doherty, Ph.D. Thesis in writing, Michigan State University. M.F. Hineman, Ph.D. Thesis in writing, Michigan State University. 162 22. Y.Y. Wu, Ph.D. Thesis in writing, Michigan State University. APPENDIX II APPENDIX II Generating Normally Distributed Random Numbers for Simulation Computer simulation of instrumental responses has long been an important part of chemistry. It is particularly valuable in research where new methods of signal processing or data analysis are to be evaluated, as with the simulations performed in Chapter IV. It is remarkable, however, that a large number of the simulations used in research and teaching inadequately model experimental noise, despite its importance in overall performance evaluation. In some cases (ref. 1, for example), experimental noise is ignored altogether, leading to unrealistic expectations for the data analysis technique under consideration. In other cases (ref. 2, for example), random noise with a uniform distribution is assumed. While this is preferable to no noise, true experimental noise typically has a normal or Gaussian distribution about the mean. The use of uniformly distributed noise can greatly alter the statistics resulting from analysis of the data and lead to the incorrect assessment of a method’s potential. The appearance of uniformly distributed noise is also distinctive and not typical of actual experimental noise. In this appendix, several methods for generating random numbers with a normal distribution are briefly described [3]. Details of the specific method used in the research presented here are given, along with a FORTRAN-77 subroutine listing. 163 164 Most computer languages support pseudo-random number generators which provide uniformly distributed values, typically between 0 and 1. Some languages support the generation of random numbers with a normal distribution, but these are not yet commonplace. The process of obtaining new random numbers from old is known as transformation and the subject has been thoroughly treated by Morgan [4]. Only transformations from a uniform to a normal distribution are considered here. Also, the discussion here assumes that the noise to be generated is white noise (i.e. a flat frequency spectrum). Correlated noise sequences may be generated from white noise by applying a digital filter with the appropriate frequency response. In describing these transformations we will adopt the notation that U(a,b) represents a uniformly distributed random variable between a and b, and N01 ,0'2) represents a normally distributed random variable with mean u and variance of 0". One of the simplest ways to generate Gaussian random variables is known as the Box-Muller method [5]. If U1 and U2 are independent uniformly distributed random variables between 0 and 1 (U(O,l)), then two normally distributed random variables can be obtained from the equations below. 3‘ I " (-2¢f2 11‘1U1)1/2 008(2'IU2) (A2.1) N2 (—262 an1)1/3 sin(21U2) (A2.2) In this equation, N1 and N2 will be normally distributed about zero with a variance of 0’2 (N(0,d’2)). The basis of this method is the transformation of a point in polar coordinates to two-dimensional 165 Cartesian coordinates. Details of, the proof are given elsewhere [4] and will not be reproduced here. The transformation is theoretically exact, however, and not an approximation. It is also very simple to implement. A disadvantage of the Box-Muller method if large numbers of random numbers are to be generated is the use of time consuming trigonometric functions. A modification known as the Polar Marsaglia method [6] avoids this problem. As before, we first obtain two U(0,1) random variables, U1 and U2, from a random number generator and then convert them to V1 and V: which are U(-1,1). This is easily accomplished with the relation V = ZU—l. It should be clear that if V1 and V2 represent Cartesian coordinates, then all of the points so defined will fall inside a box with sides of length 2 around the origin. The Polar Marsaglia method requires rejection of all points so generated which fall outside of a circle of unit radius inscribed inside the box, i.e. any case where (V13+V23)>l. If this is done, then the pairs of random variables which will remain can be easily transformed to N(0,o‘2) with equations A2.3 and A2.4. N1 V2[(-2¢f2 an)/W]1/2 (A2.3) z N u Vl[(-2¢'3 lam/mm (AZA) where W : V12+V22. The principle of this method is essentially the same as for the Box-Muller method, but the approach is somewhat different. In this case, sines and cosines of the polar angle are calculated from their definitions rather than their functional form (note, for example, that sinB : Vz/Wl/z). Incidentally, the proportion of values rejected with this method is 1 - 11/4. 166 Perhaps the fastest way to generate normally distributed random variables is to use a look-up table stored in an array. To do this, the inverse normal cumulative probability distribution is employed. A uniformly distributed variable between 0 and 1 is taken to represent a value of the cumulative normal distribution function, i.e. the area under the Gaussian distribution. By tabulating this function, the corresponding deviation from zero can be determined. For example, a value of 0.5 would correspond to a deviation of zero, 0.84 to a deviation of 16’, 0.16 to a deviation of -ld', and so on. Tabulations of cumulative distribution functions are available in most statistics texts and math handbooks. This method suffers somewhat in terms of resolution in the results obtained, since the look-up table is by its nature discrete. This can be improved through interpolation between table entries, but data for the table still needs to be entered by hand, a somewhat tedious process. A much more attractive alternative to the tabulation of the inverse normal cumulative distribution function is to use a mathematical approximation to this function, a method which has been described by Kaufman [7]. Although the inverse normal cumulative distribution function is not analytically solvable, some very good approximations to it have been given [8]. Using one of these, we can obtain a normally distributed random variable from the equations below. p : 0.5—U (A25) t : [-2 1n(ABS(p))]1/2 (A2.6) x : L-[(Co+C1L+CzL2)/(1+D1L+D2L2+D3L3)l (A2.7) N : SIGN(p)-d-x (A2.8) 167 where: U is U(0,l) Co=2.515517 D1=1.432788 01:0.802853 D2=0.l89269 C2=O.OIO328 Dar-0.001308 The result of this computation, N, is N(0,d'3), a normally distributed variable with the desired variance. Equation A2.5 transforms U to p which is U(—0.5,0.5). This is a necessary step because the approximation to the inverse cumulative normal distribution function is single-sided; i.e. it only treats the positive side of the probability distribution function. Equations A2.6 and A2.7 perform the required approximation. Equation A2.8 scales by the desired standard deviation and gives the result the appropriate sign. It is necessary to adjust the sign at this point since otherwise the equation would produce only positive values (as mentioned already, the approximation is single-sided). If the SIGN function is not available in the computer language used, it is easily implemented with a test for p>0. Since positive and negative values of p are equally likely, the random numbers will be evenly distributed on both sides of the mean. Although not indicated in the above equations, it is also necessary to check for the condition p=0 to avoid an error in eq A2.6. If this condition occurs, p can be replaced with a very small value. The simple algebraic form of this method makes it easy to implement and compatible with spreadsheet programs which are becoming more commonly used in data analysis. Although the above method is an approximation, it is an extremely good one, quite adequate for virtually all simulations. The method has been successfully employed in our laboratory for some time. Figure A2.1 shows a histogram of the actual distribution obtained from 10,000 168 0-05 firrrTTI‘l'ryir 0.04 -- 0.03 -- Probability 0.02 -- 0.01 ‘- 0-00' 'lrl'I'T'l' " -4. -3. -2. -1. O. 1. 2. 3. 4. Deviation from Mean (0) Figure A2.1. Comparison of actual (histogram) and theoretical (smooth curve) distributions obtained with the inverse cumulative normal distribution approximation. Results are based on 10,000 samples. 169 samples generated by this method. For comparison, the smooth curve shows the theoretical normal distribution. The agreement is quite good. In order to facilitate the easy addition of normally distributed noise in simulations, the FORTRAN-77 subroutine ADDNSZ was written. This subroutine is based on the last algorithm given above and will add simulated noise with a given standard deviation to a specified data array. The random number seed is specified by the user to allow different noise sequences to be generated. A listing of this subroutine follows. 00000000000000000000000000000000 0 100 0000 110 C 170 SUDROUTINE ADDNSZ(DATA,IDIM,SIGMA,ECNGRD,ISEED,IMODE) This subroutine adds randoa noise with a Gaussian distribution to the data in an array. The paraaeters involved are : DATA - is the array containing the data to which the error is to be added or the array which will contain the errors generated, depending on the setting of IMODE. IDIM - is the size of the data array. SIGMA - is the desired standard deviation of the errors to be added (value is absolute). ECEGRD - is a constant background level to be added to all of the errors generated. ISEED - is the randoa nuaber seed to be used in generating the randoa errors (INTEGERt4). IMODE - is the node of error generation : IMODE=1 indicates errors generated are to be added to the data already in the array. IMODE=2 indicates errors generated are to be put directly into the array (not added). To generate the nor-ally distributed randoa errors, this subroutine uses the randoa nuaber generator, which gives a rectangular distribution, in conjunction with the cuaulative Gaussian distribution. The randoa nuaber generated is used as the arguaent in an approxination to the inverse of the standardized cuaulative nornal distribution function (see M. Abraaowitz and I.A. Stegun, ’Eandbook of Matheaatical Functions’ (1968), Equation 26.2.23). This aethod is briefly described by L. Eaufaan in ’Trends in Anal. Chen.’, vol. 2, p. 244. Set up house INTEGER34 ISEED DIMENSION DATA(IDIM) C0=2.515517 Cl=0.802853 C2=0.010328 D1=l.432788 D2=0.189269 D3=0.001308 Zero data array if node 2 is being used IF (IMODE.EO.2) TREN DO 100 I=1,IDIM DATA(I)=0.0 CONTINUE END IF Generate errors and add to array Generate randoa nuaber. Needs to be in range 0 < P <8 0.5. If nuaber is in range 0.5-l, subtract 0.5 and aake sign of error negative. Ensure nuaber <> 0. Coapute error, scale, and add background. DO 110 I=1,IDIM ' P=RAN(ISEED) IF (P.GT.O.5) TREN SIGN=-l.0 P=P-O.5 ELSE SIGN=1.0 IF (P.E0.0.0) P80.00001 END IF T=SORT(-2.0tLOG(P)) XP=T-(C0+CltT+CZtTtT)/(l+DltT+DZtTtT+D33TtTtT) DATA(I)=DATA(I)+SIGN!SIGMA8XP+ECEGRD CONTINUE All finished, now go hoae RETURN END 171 LIST OF REFERENCES B.G. Willis, W.H. Woodruff, J.R. Frysinger, D.W. Margerum and H.L. R.Q. Thompson, J. Chem. Ed., 62, 866 (1985). P.D. Wentzell and S.R. Crouch, J. Chem. Ed., submitted for publication. B.J.T. Morgan, "Elements of Simulation"; Chapman and Hall, New York, 1984; pp. 77-136. G.P. Box and M.E. Muller, Ann. Math. Stat., 29, 610 (1958). G. Marsaglia and T.A. Bray, SIAM Rev., 6, 260 (1964). L. Kaufman, TRAC, 2, 244 (1983). M. Abramowitz and LA. Stegun, "Handbook of Mathematical Functions: National Bureau of Standards Applied Mathematics Series no. 55"; U.S.- Government Printing Office, Washington, D.C., 1964; p. 933.