APPLICATION OF THE EXTENDED KALMAN FILTER TO ENZYME REACTIONS By Henry Kwok-Hei Lai B Sc. (Physics) Imperial College, University of London, London, England A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF MASTER OF SCIENCE in THE FACULTY OF GRADUATE STUDIES DEPARTMENT OF MATHEMATICS, INSTITUTE OF APPLIED MATHEMATICS We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA August 1993 © Henry Kwok-Hei Lai, 1993 In presenting this thesis in partial fulfilment of the requirements for an advanced degree at the University of British Columbia, I agree that the Library shall make it freely available for reference and study. I further agree that permission for extensive copying of this thesis for scholarly purposes may be granted by the head of my department or by his or her representatives. It is understood that copying or publication of this thesis for financial gain shall not be allowed without my written permission. (Signature) Department of m zafiginta;. The University of British Columbia Vancouver, Canada Date DE-6 (2/88) ^2 g Abstract In monitoring biological processes, measurement of key variables is often impeded by the lack of reliable on-line measurements of biomass, substrate, and product concentrations, and the difficulty to properly model biological activity in processes that are nonlinear and time-varying. One approach to solving this problem involves the development of state estimation techniques that use available measurements to reconstruct the evolution of state variables or to estimate the bioprocess parameters. The use of filtering theory for state estimation provides a means of incorporating a deterministic model into a method of forecasting future states which includes the probabilistic uncertainties of both the system and the measuring devices. The state estimation technique we use is the discrete extended Kalman filter. This method allows the use of an approximate model and partial measurements of process variables. To study the application of the discrete extended Kalman filter to a biological problem, we used two different enzyme reactions as model systems. The first model problem is a simulation of cellulose hydrolysis with enzyme inactivation. The second model problem is a simulation of the separation of a chiral substance into its respective enantiomers. In the first problem, the hydrolysis of cellulose, a slowly varying parameter is tracked after correctly choosing the model error covariance matrix. In the second problem, we reconstruct the complete composition of the system from a partial set of measurements. Although these model problems are simple cases involving only single enzymes, the extended Kalman filter can be applied to more complex systems of enzymes. Abstract ^ Table of Contents ^ List of Figures ^ vii List of Tables ^ Acknowledgement ^ xi 1. Introduction ^ 1 2. 1.1 Enzymes ^ 1 1.2. Biological Processes ^ 2 Derivation of the Rate Equation ^ 5 2.1. The Law of Mass Action ^ 5 2.1.1 Single Enzyme-Substrate Reaction ^ 5 2.1.2. Two Enzyme System ^ 7 2.2. Quasi-Steady-State Approximation ^ 9 2.2.1. Single Enzyme-Substrate Reaction ^ 9 2.2.2. Two-Enzyme System ^ 10 2.3. Singular Perturbation Treatment ^ 11 2.4. Validity of the Quasi-Steady-State Approximation ^ 14 2.5. Enzyme-Substrate-Inhibitor System ^ 3. Filtering Theory ^ 16 20 3.1. Probability Preliminaries ^ 21 3.2. Kalman Filter Theory ^ 25 3.2.1. System Description ^ 26 3.2.2. Noise Description ^ 26 3.2.3. Derivation of Kalman Filter Algorithm ^ 27 iv 3.2.3.1. Lemma 1 ^ 27 3.2.3.2. Lemma 2 ^ 28 3.2.3.3. Lemma 3 ^ 29 3.2.3.4. Lemma 4 ^ 30 3.2.3.5. Lemma 5 ^ 33 3.2.3.6. Lemma 6 ^ 34 3.2.4. Kalman Filter Algorithm ^ 36 3.2.5. Alternative Method of Derivation of Kalman Filter ^ 37 3.2.5.1. Innovation Process ^ 37 3.2.5.2. Lemma 7 ^ 37 3.2.5.3. Theorem 5 ^ 39 3.3. Extended Kalman Filter Theory ^ 40 3.4. More About Error Covariance Matrices ^ 41 3.5. Exact Filter ^ 42 4. Enzyme Deactivation ^ 4.1. Reactions and Kinetic Equations ^ 45 46 4.2. On-Line Estimation of Enzyme Activity Using the Extended Kalman Filter ^ 47 4.2.1. System model ^ 48 4.2.2. Noiseless Measurements ^ 48 4.2.2.1.Generation of Measured Data ^ 48 4.2.2.2. Numerical Results ^ 49 4.2.3. Small Noise Measurements ^ 50 4.2.3.1. Generation of Measured Data ^ 50 4.2.3.2. Numerical Result ^ 4.2.4. Large Noise Measurements ^ 51 52 4.2.4.1. Generation of Measured Data ^ 53 4.2.4.2. Numerical Result ^ 53 4.2.5. The Choice of Model Error Covariance Matrix ^ 54 4.2.5.1. Noiseless Measurement ^ 56 4.2.5.2. Small noise measurement ^ 57 4.2.5.3. Large Noise Measurement ^ 58 4.2.5.4. Conclusion ^ 4.3. Iteration of Extended Kalman Filter ^ 4.3.1. Numerical example ^ 5. Separation of Enantiomers ^ 60 60 62 65 5.1. Enzyme-Substrate reaction rate equation ^ 65 5.2. Enantiomeric Ratio ^ 67 5.3. Conversion ^ 68 5.4. Perturbation Analysis ^ 68 5.5. Construction of Difference Equation From 1st Order Perturbative Solution ^ 71 5.5.1. Forward Propagating Equation ^ 72 5.5.2. backward propagating equation ^ 73 5.6. On-Line Estimation of Substrate Concentration by Extended Kalman Filter ^ 76 5.6.1. Defining the problem ^ 76 5.6.2. Generation of Measured data ^ 77 5.6.3. Numerical Result ^ 77 5.6.4. Filter Divergence ^ 79 5.6.5. Correction of Filter Divergence ^ 81 5.7. An Example Of Filter Divergence ^ 5.7.1. Example ^ 5.8. Estimation by Exact Filter Algorithms ^ 84 84 86 vi 5.8 A . Dynamical^and Measurement Model in One Dimension ^ 86 5.8.2. Numerical Results ^ 87 6. Discussion and Future Research 89 6.1. Discussion ^ 89 6.2. Future Research ^ 91 Appendix 1 ^ 95 Appendix 2 ^ 101 Bibliography ^ 105 vii List of Figures Fig. 1 Noiseless measurements of concentration of substrate A ^ 49 Fig. 2 Error in estimated values of A from exact measurement of A ^ 50 Fig.3 Error in estimated values of rm ^ 50 Fig. 4 Small Noise measurements of concentration substrate A ^ 51 Fig. 5 Error in estimated A using small noise measurements ^ 52 Fig. 6 Estimated rm and exact rm in small noise measurement ^ 52 Fig. 7 Large noise measurements of concentration substrate A ^ 53 Fig. 8 Estimated A and exact A in large noise measurement ^ 54 Fig. 9 Estimated rm and exact rm in large noise measurement ^ 54 Fig. 10 Absolute percentage error of estimated rm as function of covariance ^ 57 Fig. 11 Relative percentage error of estimated rm as function of covariance ^ 57 Fig. 12 Absolute percentage error of estimated rm as function of covariance ^ 58 Fig. 13 Relative percentage error of estimated rm as function of covariance ^ 58 Fig. 14 Absolute percentage error of estimated rm as function of covariance ^ 59 Fig. 15 Relative percentage error of estimated rm as function of covariance ^ 59 Fig. 16 The estimated values of rm for the original system model (4.5) and the modified system model (4.6) ^ Fig. 17 Numerical solution of rate equation of the enantiomer ^ 67 18 Exact and perturbation solutions of A ^ FIR. 19 Exact and perturbation solutions of B ^ 71 71 viii Fig. 20 Numerical solutions of substrate A from the forward propagating equations and the 4th-order Runge-Kutta scheme ^ 73 Fig. 21 Numerical solutions of substrate B from the forward propagating equations and the 4th-order Runge-Kutta scheme ^ 73 Fig. 22 Numerical solutions of the forward and backward propagating equations for substrate A ^ 75 Fig. 23 Numerical solutions of the forward and backward propagating equations for substrate B ^ 76 Fig. 24^Enantiomeric ratio calculated from the estimated concentration of substrates A and B ^ 78 Fig. 25 Exact and estimated concentrations of substrate A ^ 79 Fig. 26 Exact and estimated concentrations of substrate B ^ 79 Fig. 27 Residue of A and predicted standard deviation from the EKF ^ 80 Fig. 28 Residue of B and predicted standard deviation from the EKF ^ 80 Fig. 29 Exact and estimated concentrations of substrate A ^ 82 Fig. 30 Exact and estimated concentrations of substrate B ^ 82 Fig 31 Residue of substrate A and predicted standard deviation from the EKF ^ 83 Fig. 32 Residue of substrate B and predicted standard deviation from the EKF ^ 83 Fig. 33 Estimated enantiomeric ratio vs. conversion ^ 84 Fig. 34 The estimated concentractions of substrate B obtained by the exact filter and the EKF ^ 88 Fig. Al Numerical solutions of the forward and backward equations using (A.2.11) and (A.2.12) ^ 103 ix Fig. A.2 Numerical solutions of the forward and backward equations using (A.2.11) and (A.2.13) ^ 104 List of Tables Table 1 Absolute and relative percentage errors of estimated A and rm ^ 64 xi Acknowledge First, I would like to extend my thanks to my supervisor Prof. Robert M. Miura. Thoughout my research, he guided me at every step and spent hours on reading my thesis. My thanks also go to Dr. Tom Stokes, he not only introduced enzyme kinetics to me, but helped me with the technical aspects of enzyme kinetics thoughout my research. I am grateful to Professor Brian Seymour for his valuable time in reading this thesis. I am also grateful for financial support from the B.C. Science Council and Pacific Al in the form of the G.R.E.A.T. award. 1 1. Introduction 1.1. Enzymes Enzymes are highly specialized proteins that catalyze the chemical reactions in all living cells. They operate as organic catalysts by lowering the energy barrier for chemical reactions. Indeed, in the absence of enzyme catalysts, many biochemical processes would proceed far too slowly to maintain life. Not all enzymes are within cells, e.g. enzymes aid in the digestion of food in the stomach and intestines. For more details the reader is referred to [1]. Since enzymes are proteins, they are sensitive to temperature. Increasing the temperature generally increases the reaction rate. However, above a certain temperature, the enzyme becomes damaged or denatured, and loses its catalytic ability. All chemical reactions are reversible. If an enzyme can break a substrate (a substrate consists of the molecules that the enzymes act upon) A into its products B and C, then the same enzyme is equally capable of catalyzing the reverse reaction. Enzymes are true catalysts and do not alter the chemical equilibrium. The enzyme catalytic action merely reduces the time needed to reach this state. There are several theories to explain why enzymes are specific and how they operate. All of these theories focus on the three-dimensional structure of the enzyme. The simplest theory assumes the so-called "lock-and-key" hypothesis. It postulates that a substrate molecule (the "key") attaches itself to an active site (the "lock") of an enzyme molecule, forming a temporary complex. The active site has a particular shape, so only a substrate with the complementary shape can attach itself to this site. Similarly, a lock only accepts a key with the correct shape. Hence, molecules with different shapes cannot attach to the active site. Excessive heat causes the structure of enzyme 2 to change. This distorts the active site cause the structure to change. This distorts the active site causing a loss of enzyme activity. (Note: a denatured enzyme may still bind to the substrate.) That is why the catalytic ability of an enzyme is decreased if the temp erture is raised too much. However, there may be other compounds which are close enough in shape to the substrate to fit into the active site of the enzyme. These alien substrates react slowly and compete with the true substrate for active sites. This inhibits the enzyme's activity - the ability of the enzyme to combine with the substrate to form an enzyme-substrate complex (enzyme and substrate are held together by physical forces). There are two types of inhibition: reversible inhibition and irreversible inhibition. In reversible inhibition, an alien molecule only forms a temporary bond with an enzyme. In irreversible inhibition, an alien inhibitor molecule either permanently blocks an enzyme's active site or affects this site by binding to a site elsewhere on the enzyme. Many inhibitors, particularly of the noncompetitive type, act as poisons. Cyanide, which binds with an enzyme necessary for cellular respiration, is a good example. 1.2. Biological Processes Monitoring and controlling biological processes is often impeded by the lack of reliable on-line measurements of biomass, substrate, and product concentrations, and by the difficulties in properly modeling biological activity in processes that are nonlinear and time-varying. One approach to solving these problems is to develop state estimation techniques that use available measurements to reconstruct the evolution of state variables or to estimate the bioprocess parameters. The use of filtering theory for state estimation provides a means of embedding a deterministic model into a method of 3 forecasting future states which incorporates the probabilistic uncertainties of both the system and the measuring devices. The state estimation technique we use is the discrete extended Kalman filter. This method allows the use of an approximate model and partial measurements of process variables. To study the application of the discrete extended Kalman filter to a biological problem, we use two different enzyme reactions as model systems. The first model problem is a simulation of cellulose hydrolysis with enzyme inactivation. The second model problem is a simulation of the separation of a chiral substance into its respective enantiomers. In the problem of hydrolysis of cellulose, the value of a slowly varying parameter is estimated from a set of measurements that assumes the model equations of the dynamical system contain errors. For this problem, we have studied the effects of different amounts of noise in the measurements and iteration of the extended Kalman filter. The accuracy of the estimated value of the slowly-varying parameter decreases as the noise in the measurements increases and, cannot be improved by a better choice of model error covariance matrix (the variance of the error in the dynamical equations). Secondly, we find iteration of the extended Kalman filter on its output does not improve the accuracy of the estimated value of the slowly-varying parameter. However, we show in Chapter 4 that iteration of the model equations from the output of the extended Kalman filter could improve significantly the accuracy of the estimated value of the slowly-varying parameter. In the second problem, we reconstruct the composition of a chemical system from a partial set of measurements using model equations of the dynamical system which contain no error. Both the extended Kalman filter and 4 exact filter are applied to this problem. In applying the extended Kalman filter to monitor this process, we find that if the model error covariance is chosen to be too small, then it will cause the filter to diverge. This phenomenon is known as filter divergence which is discussed in Chapter 5. In applying the exact filter to this process, we find a faster rate of convergence to the exact value than the extended Kalman filter when there are errors in the initial concentration of the enantiomers. However, it is feasible to compute the exact mean and variance of these estimated concentrations only in scalar problems. Although these model problems are simple cases involving only single enzymes, the extended Kalman filter can be applied to more complex systems wity several enzymes. 5 2. Derivation of the Rate Equation Enzyme kinetics is a tool for studying enzyme mechanisms and the rate equations, which represent the behavior of the enzymatic reactions, are generally expressed as a set of nonlinear differential equations. To obtain the analytical solution of these rate equations is generally impossible. This makes it necessary to introduce approximations in order to analyze the experimental data in enzyme kinetics. One of the most common approximation schemes is the quasi-steady-state approximation. In this chapter we will look at the quasisteady-state approximation from the singular perturbation point of view [2] and examine the validity of this approximation scheme [3]. 2.1. The Law of Mass Action Chemical kinetics deals with the rate of chemical change of moluclues in a reaction system. A chemical reaction in this system leads to an instantaneous difference in concentration of reactant molecules. A chemical reaction occurs when two or more molecules collide under appropriate conditions. According to [4], trimolecular reactions are unlikely. The law of mass action states that the rate of molecular collisions of two chemical species in a dilute gas is proportional to the product of the concentrations of the two chemical species. As a result a rate equation of a chemical reaction can be derived based on the law of mass action. Throughout this chapter, we adopt the following notation. Roman capital letters represent chemical species and a square bracket around a Roman capital letters represents the concentration of that chemical species. 2.1.1. Single Enzyme Substrate Reaction - The study of the single enzyme-substrate reaction dates back to 1902 when Brown (see Fersht [1]) published his work on the enzyme-catalyzed hydrolysis 6 of sucrose by invertase. In 1988, Segel [3] derived a validity condition for the quasi-steady-state approximation of this reaction. According to Michaelis and Menten in 1913 (see [1]), the basic enzyme-substrate reaction is described by the reactions k, E+ S ^ C k2 >E+P where E is the enzyme, S is the substrate, C is the intermediate enzymesubstrate complex, P is the product, and kb k_1, k2 are the rate constants. The rate equation obtained by using the law of mass action are d[E] dt d[S] dt d[C] dt d[P] dt = k1[E][S]+k_1[C]+k2[C], (2.1a) = k [E][S] + k„[C], (2 .1b) = k1[E][S] — k 1 [C] — k,[C], (2.1c) = k2[C], (2.1d) with initial conditions [E](0) =E0, [S](0) =S0, [C1(0) = 0, [P1(0) = 0. We note that the system (2.1) is nonlinear and cannot be solved analytically. However, there are two conservation laws given by [NW +[C](O = E0,^ (2.2a) [S](t)+[C](t)+[P](t)= So.^ (2.2b) Employment of (2.2) to eliminate [E] in (2.1) yields the following simplified system of differential equations, 7 d[S] dt (E0 — [CP[S] + k_, [C], (2.3) d[C] ki (Eo — [C101— (1c-1 k2)1C1, dt IS1(0) = S0, [C](0) = 0. Note that (2.1d) can be solved fpr [P] once [C] is known, or [P] can be obtained algebraically from (2.2). The system (2.3) is still nonlinear and cannot be solved analytically. However, in the next section a quasi-steady-state approximation can simplify this system to a manageable form. 2.1.2. Two Enzyme System An example of a two enzyme-substrate reaction is cellulose hydrolysis which produces glucose from cellulose and in which there are at least two different enzymes involved in the reaction. To simplify the discussion, we consider the following schematic reaction. A+E k., B+E k3 B+E^BE 1_13_4 ki B +E' X2 k2 > C E' CE' < k_,3^ where A, B, and C are the substrates; X1, X7, BE, CE', are the enzymesubstrate complexes; kl, k..1, k7, k3, k_3, k'1, k'.1, k'2, k'i, k'_3 are the rate constants. The rate equations obtained by using the law of mass action are d[A] = —ki[A][E]+ k_i[Xi], dt (2.4a) 8 d[E] dt = k,[E][A]+k i[Xi] + k21X1]— k3[E][B1+ k_3[BE],^(2.4b) d[X1] —ki[A][El— (k_, + k2)[X,],^ dt d[B] k,[Xi] + k 3[BE] —k3[BE]— k;[B][El+ k'1[X2], ^(2.4d) dt d[E'] dt (2.4c) = k[0][B]+k'1[X,]+k;[X,]— le3[C][E']+ k3[CE1,^(2.4e) d[X,1 k, [Elm _ k'^ — k; [X2],^ dt d[C] dt k;[X2]+1e_3[CE1— k[C][E'],^ d[BE] = k [B][E] — k 3 [BE],^ dt 3 d[CE'] dt le3[C][E']— k' 3 [CET^ (2.40 (2.4g) (2.4h) (2.4i) with the initial conditions: [A](0) = Ao, [B](0) = B0, [C](0) = 0, [E](0) = Eo, [E1(0) = [BE](0) [CE1(0) = [X1](0) = [X2](0) =0. Although the system (2.4) is a set of complicated nonlinear equations, there are two conservation laws which can be verified by a direct application of the differential equations (2.4). Adding (2.4)(b) ,(c), (h) and integrating, we get the first conservation law [E] + [X,]+ [BE] = Eo.^ (2.5) Similarly, adding (2.4)(0, (g), (i) and integrating, we obtain the second conservation law [El F[X2]+[CET= Po.^ - (2.6) Since (2.4) have many parameters and is very complicated, we will apply the quasi-steady-state approximation to simplify it. 9 2.2. Quasi-Steady-State Approximation The quasi-steady-state approximation was proposed by Briggs and Haldane in 1925 (see [1]). The steady-state treatment is based on the hypothesis that the concentration of certain reaction intermediates, such as the enzyme-substrate complex does not change rapidly during the course of the reaction. Mathematically speaking, if [ES] is the concentration of the intermediate enzyme-substrate complex then d[S]idt))dIESVdt 0. In the next section, a singular perturbation technique will be applied to justify this statement for a single enzyme-substrate reaction. Below we will apply the quasi-steady-state approximation to the rate equations for single and double enzyme-substrate reactions developed in the last section. 2.2.1. Single Enzyme-Substrate Reaction Since the rate of the change of the enzyme-substrate complex is close to zero, as a first-order approximation, we can set d[C]/dt=0 in (2.1c) and use the conservation law (2.2a) to express [E] in terms of [C] and E0. Then [C] can be expressed as [Cl ^k, E0 [S] (k + k2) + k, [S] (2.7) Substituting (2.7) into the rate equation (2.1d), we get dill _ V [S] ^dt^K. + [S] (2.8) where K. = (k2 + kilki and V =k, Eo. Finally combine d[C]idt=0 and the conservation law (2.2b) to obtain ^d[S] ^d[P]^V[S] ^dt^dt^K. +[S] • (2.9) This equation can be solved exactly using the method of separation of variables and yields the implicit solution 10 1 (S, [S])+ (se —LS-1)}. t = --{K lnV V^m^ (2.10) 2.2.2. Two-Enzyme System To apply the quasi-steady-state approximation to (2.4), we have to set the rates of change of the concentrations of all intermediate enzyme-substrate complexes to zero, i.e., d[X1] = d[X2] d[BE]d[CE'l 0 dt^dt^dt^dt (2.10 The trick is to express the concentrations of all the intermediate enzyme complexes in terms of [A], [B], [C], [E], and [El. From 2.4(c), (f), (h), (i), and (2.11), we get 1^[A][E], [X1] = ^ k_, +k2 (2.12a) k; [X21= ^ [B][E], ' +k', (2.12b) [BE[ = k3 [B][E], k_3 (2.12c) k' [CE']= 3 [C][E'[. (2.12d) If we combine the conservation law (2.5) with 2.12(a), (c), then [E] can be expressed in terms of [A], [B], and the rate constants to give [E] = Eo 1+ [A] + [B] K, (2.13) where K. (k_, +k2)/k, and Ki = k_3/1c3 .^Similarly, ^combining^the conservation law (2.6) with 2.12(b), (d) gives 11 [E] = [B1 [C1 1+ ^ + [Cm K: ^(2.14) where K'm (k'^k;)/k; and K: = k'3/k; . To complete OUT calculation, substitute (2.13), 2.12(a) and (2.14), 2.12(b), into 2.4(a) and 2.4(d), respectively. The result is d[A]^I-. [A] ^(2.15a) dt [A]+ K (1+ [13] Ki) d[C]^[13] dt [13] + K' 11 + ^ K' (2.15b) ' where rm. = k2 Eo, and rr:, = k; E. Finally, applying the two conservation laws (2.5) and (2.6) as well as (2.12) to 2.4(d) gives another conservation law for the substrates, namely [A] + [B]+ [C]=constant.^ (2.16) This completes the derivation of the simplified rate equations. However, the conservation law (2.16) is distinct from the other two conservation laws in that it is a consequence of the quasi-steady-state approximation. 2.3. Singular Perturbation Treatment The implication of using the quasi-steady-state approximation to derive expressions for the enzyme mechanisms can be examined by considering the validity of the approximation when it is applied to the simplest enzymecatalyzed reaction E + S^C k2 >E + P (2.17) 12 Even if the quasi-steady-state approximation works for (2.17), there is no guarantee that the approximation will be valid for a more complex mechanism. However, if the quasi-steady-state assumption is a poor approximation for (2.17), then it wll probably be a poor approximation for more complex enzyme mechanisms. Mathematically speaking, the quasi-steady-state approximation changes the character of the mathematical problem from a system of first-order ordinary differential equations to a system of first-order ordinary differential equations coupled to an algebraic equation. It seems prudent to apply a careful analysis to the rate equations in order to determine when this approximation can be applied. In 1967, Heineken, Tsuchiya, and Aris [2] justified this assumption by using singular perturbation theory. In the following, we borrow their idea but we use a different non-dimensionalization scheme. Applying the conservation law (2.2a) to (2.1b) and (2.1c) one obtains d[S] dt k (E0 — [C])[S]+ k [C],^ d[C]^(k^ k2)1C1+ ki(E, - [C])[S], dt (2.18a) (2.18b) with the initial conditions [S](0) = So, and [C](0) = 0. We nondimensionize of (2.18) by introducing the following dimensionless variables [S][C]^E0 S=^, C=^6=^ K^Eo^K = kiEot. (The dimensionless substrate concentration s, according to Heineken et al. [2], is s=[S]/So. Our transformation, results in only one parameter, a, in the dimensionless equation in contrast to two parameters obtained in their 13 dimensionless equation. However, their initial condition on s is normalised to 1, but our initial condition depends on the parameter so. In conclusion, there is still a total of two parameters in both dimensionless formulations, but for the purposes of deriving the quasi-steady-state approximation, the dimensionless transformation employed in this chapter simpifies the algebra.) In terms of the dimensionless variables, (2.18) becomes ds — = a c – s (1– c), dr 6 dc = s – c (1+ s), dr where^a =^K. , s(0) so, and c(0) = 0. In the language of singular perturbation theory, we assume (2.19) governs the long time behaviour and define a short time scale, T, by T=6 T. In terms of T, (2.19) becomes ds ^ = c(a c – s (1– c)),^ dT (2.20a) dc = s – c (1+ s).^ dT (2.20b) A complete perturbation solution can be obtained by the method of matched asymptotic expansions. To derive the quasi-steady-state approximation, we need only the leading-order perturbation solution of (2.19) which is obtained by setting c=0 in (2.19). Solving the resulting algebraic equation, we obtain c= (2.21) 1+s Substituting c into (2.19a), we get ds dr k,2 k, +k s l+s If the definitions of s, c ,and (2.22) T are used in (2.22) and (2.21), then they become 14 [C] = E [S] 0 , + [S1. d[S]^V [S] dt^K+ IS] • (2.23) This system is a good approximation of (2.18) provided e—E0/K,, is small. (The small parameter obtained by Heineken et al. [2] is 6=E0/S0. Consequently, the criterion for the quasi-steady-state appeoximation is E0/S0<<1. However, in the next section, we show that both criteria can be interpreted as a first-order perturbation expansion of the criterion + So) «1 if either So or IC, is small.) This perturbation argument is difficult to apply to more complicated enzyme systems. In the next section, another approach based on the special character of singularly perturbed systems is employed to determine a condition for the validity of the quasi-steady-state approximation. 2.4. Validity of the Quasi-Steady-State Approximation According to the singular perturbation theory [5] used in the last section, the two different time scales in (2.20) and (2.22) correspond to the inner and outer regions, respectively. In this problem the inner region is associated with a fast initial transient period in which the substrate concentration does not undergo any significant change. On the other hand, the outer region is associated with the steady-state period in which the enzyme-substrate complex concentration does not undergo any significant change. The basic hypothesis is that for an enzyme-catalyzed system, the substrate concentration does not change during the fast initial transient period and the enzyme complex does not change during the steady-state period [18]. We can apply this hypothesis to a more complicated enzyme-catalyzed reaction, for which the perturbation argument is difficult to apply, to 15 determine the condition for the validity of the quasi-steady-state approximation. An example is the enantiomer reaction described in Chapter 5. For now, we apply this hypothesis to derive a condition for the validity of the quasi-steady-state approximation for the simple enzyme-catalyzed reaction (2.17). From the system (2.18), since the substrate concentration, [S}, does not change during the fast initial transient period, 0<t<tc, we can estimate the order of of magnitude of tcs by putting [S]=So in (2.18b). This makes (2.18b) a linear equation with one unknown, [C]. The solution is given by [C] = ^ S0 (1 exp(—X t)), K. +So where X k1(S0 +K.). Consequently, the fast initial transient period, tc, is given by tc=X•1, i.e., 1 tc k1(S0 +K.)^ (2.24) Having obtained te, we can proceed to estimate the decrease in substrate concentration, AS, during the fast transient period. The criterion from the hypothesis can be formulated mathematically as, AS «1^ So (2.25) where So is the initial substrate concentration. However, if the inequality !IASI. /Sol «1 is satisied, then (2.25) is true as well. ButIASI. = tc • idS/dtl...(1AS I is the order of tc •IdSidtlin.). For the class of solutions that d[S]idt 0, (2.20) imply IdSldtl. Eo So. (For the existence of this class of solutions, see Appendix 1.) Combining these equations with (2.25), we obtain the criterion for the quasi-steady-state approximation, namely 16 E0 ^« 1. (2.26) K. +So The criterion (2.26) indicate that the quasi-steady-state approximation is valid when EQ/K. «1 is not true provided that So is large enough or when Eo/S,, «1 is not true provided that Km is large enough. Furthermore, E E K. +So °^ «1^ ^«1. Thus the criterion E0/(K. + S0) «1 takes into account both the criterion derived by Heineken et al [2] and the one derived in the last section. 2.5. Enzyme-Substrate-Inhibitor System In this final section we derive the approximate equations describing the decrease in enzyme activity due to the presence of an inhibitor and show that this mechanism of enzyme deactivation is not equivalent to the unknown mechanism involved in the simulated enzyme deactivation in Chapter 4. Roughly speaking, inhibition is the inactivation of a enzyme in a major metabolic pathway by another substance. Inhibition of an enzyme can be either reversible or irreversible. When inhibition is reversible, the activity of the enzyme can be regained by removal of the inhibitor by physical or chemical means. Irreversible inhibition, however, is characterized by complete loss of enzyme activity over a period of time even in the presence of low concentrations of the inhibitor and the enzyme activity cannot be regained by any physical treatment. Consider a competitive inhibition reaction, E + S^C, k3 E+1^C„ k2 E+P 17 The rate equation obtained using the law of mass action are d[S] = —ki[E][S]+ k_1[C1], dt (2.27a) d[E] dt (2.27b) k, [E][S] — k_3[E][I] + k_,^+ k_3[C,] + k2[C1], d[11 — k,[C,], dt (2.27c) d[C1] = [Nisi _ ki ^ — k2[C, ], dt (2.27d) d[C,i = k3[E][il— k_3[C2], dt (2.27e) d[I]^k3[E]ri±k3ic2i. dt (2.270 with the initial conditions [S[(0) = So, 111(0) = 0, [I](0) = Io, [C,[(0) = 0, [C2](0) = 0, [E](0) = E0. Combining 2.27(b), (d), (e) and (e), (0, respectively, give the two conservation laws [E]+ [C1] + [C2] = constant, ^ [II+ [C2] = constant. ^ (2.28a) (2.28b) We apply the quasi-steady-state approximation, i.e., set d[C,]/dt = 0 and d[C21/dt = 0, to obtain E0[S] [C,J= K 11+ ^)+[S] K [E] = K [C,] m [S] (2.29a) (2.29b) where K. = (k + k2)/k1^K, = k3/k3. Combining (2.29) and 2.27(c) gives, 18 d[P]^V[S] dt K1+ ^ 1+ [S] K (2.30) where V = k, E0. Substituting 2.29(b) into 2.27(a) and using the definition of Km, we have another conserved quantity [P] + IS]=constant. Using d[C21/dt=0 in the second conserved quantity in (2.28b), we have [I] = Ie so that d[S]^V[S] dt K. 1+ ° j+[S] K, (2.31) Comparing (2.31) to (2.9), we find that the only difference between these two equations is the additional multiplication factor in the denominator. Since Kin(1 +Io/K,) > Km, the rate of decrease in the substrate concentration is slower when inhibitor is present. The simulated enzyme deactivation phenomena that is discussed in Chapter 4 can be understood in terms of an inhibition mechanism. However, the mathematical equations in these two cases is not consistent and deny the interpretation that the mechanism of the enzyme deactivation discussed in Chapter 4 is competative inhibition. In indeed, consider I0/K1 as a small parameter, 6, then (2.31) can be rewritten as d[S] _ ^V[S] dt^K.(1+)+ [S] [S] K. + [S1' where (2.32) 19 ( 1 K 1+ 6 ^ m K. +IS] , The differential equation for V is given by dV ç7 2 Km1S1 dt^+1SlY Compare this with the differential equation drin/dt = (2.33) -6 rin the inhibition in the cellulose system discussed in Chapter 4. It shows that the first-order system comprising (2.32) and (2.3) is not equivalent to (4.4), cf. Chapter 4, so we conclude the underlying mechanisms of the enzyme deactivation of these two cases are different. 20 3. Filtering Theory In general terms, "filtering" means to remove something that we do not want to pass through a barrier, e.g., water purification. When the entities involved are signals, such as an electrical voltage, the barrier becomes a filter in the sense of signal processing. To determine whether a system is performing properly, and ultimately to control the system's performance, we must know what the system is "doing" at any instant of time. In other words, we must know the state of the system; for example, in navigation, the state consists of position and velocity of the craft in question. In a bioreactor, the state may be substrate concentrations and enzyme activity. When we build measuring devices and take measurements of a system, the measurements are generally contaminated with noise caused by the electronic and mechanical components of the measuring device. The problem of determining the state of a system from noisy measurements is called estimation, or filtering. The statistical approaches to filtering postulate that certain statistical properties are possessed by the useful part of the signal and by the unwanted noisy part of the signal. Measurements are the sum of the useful part and the noisy part of the signal, and the task is to eliminate as much of the noise as possible through processing of the measurements by a filter. The earliest statistical ideas of Wiener [6] and Kolmogorov [7] concern processes with statistical properties which do not change with time, i.e., stationary processes. This assumption that the underlying useful signal and noise processes are stationary is crucial to the Wiener and Kolmogorov theory. However, a theory that does not require this assumption was developed by Kalman in his celebrated paper in 1960 [8]. The new theory is called the Kalman filter and is an optimal state estimation process which is applied to a dynamical system that 21 involves random perturbations. More precisely, the Kalman filter consists of a linear, unbiased, and minimum error variance recursive algorithm to optimally estimate the state of a dynamical system from noisy data taken at discrete realtime intervals. It has been widely used in many areas of industrial and governmental applications such as laser tracking systems, satellite navigation, ballistic missile trajectory estimation, radar, and fire control. In this chapter, we are concerned only with Kalman filter theory. This theory will be applied to enzyme systems in the next two chapters. The unifying theme in this chapter is the probabilistic or Bayesian approach. To derive the Kalman filter of a linear dynamical system a Hilbert space approach is possible [9], but it is difficult to generalize this approach to a nonlinear dynamical system. In the next section, we will summarize some results from probability theory that are required to understand and interpret the Kalman filter. Further information on probability theory can be found in [10], [11]. 3.1. Probability Preliminaries One of the fundamental concepts in the Kalman filter is the gaussian process, so we will start with it. 1. Gaussian Random Vector Let X be a random n-vector with expectation denoted by X = E[X] and covariance matrix E = ENX — XY[X — X1)]. If E has a nonsingular covariance matrix, i.e., det(E) 0, then we say that X is gaussian or normal if and only if its probability density function is of the form Px (x) = 1 ^F 1^ (27cr idet Er2 exPL^E-1(x 22 For a gaussian random vector, X, with mean X and covariance matrix E, we write X—N(X, E). If we denote (¢x to be the characteristic function (the Fourier transform of the probability density function) of X, then (I)„ is evaluated to be (1)x = exp(j sT X — —1 sT E s), 2 where j is -sri. 2. Theorem 1 Let X and Y be random n-vectors with Y=f(X). Suppose f-1 exists and both f and f I are continuously differentiable. Then the probability density function for X and Y are related by (a-1(0) Py (Y) Px (f-1(Y)) det ay ) (3.1) where det is the determinant. (See [12] for the proof of Theorem 1.) This result is important in Bayesian probability theory to compute the conditional probability density function conditional on a set of measurements. It is used in the proof of Theorem 4 which is used extensively in the derivation of the Kalman filter algorithm and is used in deriving the exact filter. 3. Jointly Gaussian Distribution To say that random n-vector X and random m-vector Y are jointly gaussian random variables is the same as saying that the random vector ZT-IX, YIT is a gaussian random vector. The probability density function and the characteristic function are given by 1^ p(z) =^ ,,, e^— x (27t)(m+n)l2 det(Rz r (3.2) 23 (1)z(u) = exp(juTz — .1uTRzu), ^ (3.3) where ZE91m+1, ZT E{ZT = [E{X11,E{YT^[XT^and Rz = ERz — 2)(z — 2)1 [E{(x — R)(x — TOT] Eit(x — g)(y — E[(y — V)(x — ERy — V)(y — Rx Rxy 4RYX Ri. 4. Theorem 2 If the jointly gaussian distributed random vectors, X and Y, are uncorrelated, they are independent. (See [12] for a proof of Theorem 2.) Theorem 2 and Theorem 3 below, are used in the proof of theorem 4. 5. Theorem 3 Consider the p-vector^Rz), and let w=CZ+a, where C is a qxp constant matrix, a is a constant q-vector, and w is a random q-vector. then . w—N(C- Z+a,C-Rz.CT).^ (3.4) (See [12] for a proof of Theorem 3.) This theorem essentially means that linear operations on gaussian random vectors produce gaussian random vectors. 6. Theorem 4 Let X and Y be jointly gaussian distributed random vectors of dimension n and m, respectively . Then the conditional density of x given Y=y is gaussian with mean X + Rxy R-y1(y — Y)^ (3.5) and covariance matrix Rx — Rxy 12-y1 Ryx^ (3.6) 24 This theorem is essential in the derivation of the Kalman filter and so we present its proof below. In the next section, we will show how to obtain the Kalman filter algorithms by repeated application of this theorem. Proof This proof combines elements found in different sources [12], [13]. Consider the transformation = X — RxyR-y1Y, (3.7) W, = Y. From Theorem 3, WT [WiT,W21 is gaussian distributed. The covariance matrix of W1 is — E[Wi])(W, — E[', DT I = Rx — Rxy Ry'R, . An easy computation shows that W1 and W2 are uncorrelated. Hence, they are independent by Theorem 2. As a result, their joint probability density function is the product of their individual densities N(X — Rxy Ry^Rx — Rxy Ry'Ryx W, NI(Y,Ry). The joint probability density function is given by Pw1,w2(w1,w2) [(2 det(Rx — R xy R D )11 ^ •exp{-4-(w1 — X — Rxy Ry'^(Rx — Rxy Ry'Ryx )(w, — X — • [(27t)m12 det(Ry )1/2 Rxy Ry-1Y)] (3.8) exp[— (w, —^Ry-1(w2 — The jacobian of the transformation (3.7) is 1, and so according to Theorem 1, we obtain the joint probability density function px,y(x,y) by eliminating W1 and W2 in (3.8) in favor of X and Y via (3.7). As a result we get 25 px,y (x,Y) = {(271)212 de Rx — Rxy Ry-lRyx )1121 exp{— (x — X — Rxy Ry-1(y — Y))T (Rx — Rxy Ry'Ryx ) 1(x — X — Rxy Ry-1(y — Y))1 • [(27c)m12 det(Ry )1/2]^„(y — Y)T Ry-1(y — 17)1 ^ (3.9) The second probability density function in (3.9) is the probability density function of y. The conditional probability density function of x on y is defined by PxtY (xIY) = P" (X' Y) (3.10) PY (Y) Comparing (3.9) and (3.10) the conditional probability density of x on y is given by Px (x ly) = [(27cf2 det(Rx — Rxy Ry -exp - (X - X Rxy YX 1/21 1 Ry-1(y — 17))T(Rx —12„yRy1Ryx) 1(x — X — R„yRy-1(y— y))1.(3. 1 1) 3.2. Kalman Filter Theory In this section, the linear Kalman filter algorithm will be derived. Before the derivation, we define the notation used in the filter equations: X denote the mean of a random vector X. Zk-ldenotes the set of measurement vectors^,} where Zi for i=1,2, ...,k-1, is the q-vector representing the measurements at time ti, denotes the mean of the conditional probability density of Xk conditioned on Z', e.g., Xk,,_, predicts the state at time tk from measurements up to time t1_1, whereas )(km updates the predicted state at time tk using all measurements up to time tk, 26 Rx, denotes the covariance matrix of random vectors X, Y, with 12„., = E[(X — 'R)(Y P, denotes the covariance matrix of a random vector X, conditioned on the measurement set Z'. 3.2.1. System Description We shall restrict our attention to discrete-time systems, or, equivalently, systems where the underlying system equations are difference equations rather than differential equations. In some cases, the difference equations correspond to discretizations of differential equations. The reason for studying a discretetime system is that observations are made and control strategies are implemented frequently at discrete times in a real system. The dynamics of the state vector X is described by the following system of difference equations: F„ X, +^k=0,1,2,3,...,^ (3.12) where Xk is a random p-vector, Fk is a pxp matrix independent of Xk, and is a random p-vector representing the model system error. The statistical properties of k are described in the next section. The measurements of the system at time tk are described by the following system of algebraic equations Z, = H, X, +^ (3.13) where Hk is a qxp matrix independent of Xk, Zk is the vector of quantities measured at time tk, and flk is a random q-vector representing the model system error. The statistical properties of 1k are described in the next section. 3.2.2. Noise Description To derive the Kalman filter (KF), we have to make assumptions about the noise processes k and Th. 27 Assumptions on the noise processes The noise processes {k} and {Ilk} are white noise, i.e., E.k and ilk are gaussian variables with zero mean and EN,A1= ö ^= ö ^E[rik0 =^0^(3.14) where Sk and Qk are qxq and pxp matrices, respectively, and Su is the Kronecker delta. In most of the literature, the covariance matrix, Qk, error covariance matrix and the covariance matrix, is called the model Sk, is called the measurement error covariance matrix, respectively. 3.2.3. Derivation of Kalman Filter Algorithm In this section, we derive the KF algorithm using Theorem 4 which is based on Bayesian probability theory. This derivation is quite long and we will break it up to into several lemmas. Because of the recursive nature of the resulting formulas, i.e., the index appears only in the subscripts, it suffices to prove the lemmas for k=0 and 1. Each lemma summarizes an important result. In breaking up the derivation of the KF, our objective is to use the known measurements of a system to both update its present state and predict its future state. In Lemma 1, we develop formulas to update the mean and covariance of the system state at time to based on measurements up to time to. The distinction between the formulas in Lemma 1 and the updating formulas in Lemma 5 is that the formulas in Lemma 1 require us to supply the initial value covariance matrix. Lemmas 2, 3„4, and 6 give results predicting the state at time tk when the measurements up to time tk_1 is given. 3.23.1. Lemma 1 The mean X010 and the covariance matrix P010 of the conditional probability density function of X0 conditioned on Zo are, respectively, ^ ^ 28 X010 = X0 +Rx0x0 HoT(H0 Rx0x0 HoT +s) '(z0 —Ho XJ,^ (3.15) \-1 H O T(H R^HT+S)^R P010 H 0^ XoX0^X0X0 ^0^X0X cuo^^—R^ 0^0^ X0X0 whereXo is the mean of the probability density of Xo. Proof Applying Theorem 4 from Section 3.1, the conditional mean X0/0 and the conditional covariance matrix are given symbolically as X010 — X0 + RX0Z0 R;0 (z0 20 )9 P = R — R^R 0/0^XoX0^X9Z0^ZOZO^Z0X0 (3.17) (3.18) We need to evaluate R^ and R from the dynamical and measurement zoz,^xozo model equations (3.12)and (3.13). Firstly, let us evaluate Rx0z0 from its definition Rx0z0 = E[(X0 — X0)(z0 — 4)1 Applying the measurement model equation we get Rx0z0 = E[(X0 —^ )(H0(X0 — X0) +ri )1' . Since X0 and 10 are independent, Rxozo can be simplified to Rx0z0 Rx0x0 H°T (3.19) Similarly, if we apply the definition of R zoz 0 , measurement model equation and E[rioriol = So, we obtain Rz0z0 = Ho Rx0x0 HoT +So. ^ (3.20) Finally, 12,0x0 is the transpose of Rx0,0 12,0,0 = Ho Rx0x0.^ (3.21) Substitution of (3.19), (3.20), (3.21) into (3.17), and (3.18) will complete the proof. 3.2.3.2. Lemma 2 29 The mean of the probability density function of X1 conditioned on Zo is given by X„0 = F0 X 010 ^ . (3.22) Proof Applying Theorem 4, the mean of the probability density function of Xi conditioned on Zo is given by X,10 =X + Rxizi R;lozo (Zo — Z0). ^ (3.23) Firstly, X, is the mean of X, and can be shown to be equal to Fo X° from the dynamical model equation. Secondly, applying definition of Rxizo to the dynamical and measurement model equations. Calculations similar to those in the proof of Lemma 1 show that R^= F0 R XoZo HT0 '^ XiZo^ (3.24) Substituting (3.24) and X, = Fo Xo into (3.23) and using Lemma 1 gives the desired result. 3.2.3.3. Lemma 3 The covariance matrix of the probability density function of Xi conditioned on Zo is given by Pito Fo P010 FOT + Q0' (3.25) Proof Applying Theorem 4 the covariance matrix of the probability density function of X1 conditioned on Zo is given by Puo = R — R^R XiZo^ZOZO^ZOX1 •^ (3.26) To get an expression for Puo, we have to calculate Rx1x1 since the other terms in (3.26) have been calculated in Lemma 2. From the definition of Roc, x and the dynamical model equation, Rxix, = E{(Fo X0 + —F0 X0)(F0 X0 + —F0 Xo)Ti. 30 Rearrangement of terms and the independence of X1 and 40 implies Rx1x1 Fo Rx0x0 FT +Q0, where Q0 = E[001. Substitution of Rxi„,. R10, R;lozo into (3.26) and application of Lemma 1 give the desire result. 3.2.3.4. Lemma 4 Let X1, Z1, Zo be jointly gaussian, then the mean of the probability density (X 1/0 ) function of the random vector (Xi) conditioned on Zo has mean ^and ‘, Z 1/0 covariance matrix T P1/0 H, Pito HiT + (3.27) where S, = Proof From Theorem 4, the mean of the probability density function of Z j conditioned on Zo is given by =^+ Rzizo 1Z-zioz0 Rz0z1. (X, Also, consider I " 1 as a random vector, then from application of Theorem 4, the mean of the probability density function ofrX1) conditioned on Zo is Z, given by rzii E zo (3.28) 31 However, from our notation R X,^X, is equal to E[{( zi —^11{Zo —11. zo^ J ) (::) Hence, (X — X )(Z — Z^(Rwo R^= E[^° (P°^(Zi^Rzizo (3.29) Substitution of (3.29) into (3.28), we obtain Erz: zo^(Xz,i )1_ (RR:,z0)1Z-z10,0(Z0 —20) izo (Xi) rxizo R-,10zo(Z0 — Z0) Rzizoll-z'ozo(Zo — Comparing the above expression to X110 and Z110, we obtain z0= (3.30) "4110 By assumption for the KF, rz: , Zo are gaussian. Hence, the probability (X, density function of^conditioned on Zo is also gaussian. The mean of this Z, conditional probability density function is given by (3.30). Now we calculate the covariance matrix of the probability density of zi conditioned on Zo. Applying Theorem 4 again, and denote. the covariance matrix of the probability density function of^ (X 1 conditioned on Zo by a, Z, 32 fRx zo Similar to the argument for deriving Rri) =^, z, zo R (XzitXzi) = (R^Rx,z1) xix, Rzixi Rzlzi) Substitution of the above expression into a and after simplification, a is shown to be a— —R R-1 zoz, R ) R.-1^R^ [R^ x,x,^Rx 120^Z020^ZoX 1 R xiz,^x,zo^zozo R — R^R z1z0 zozo 7ox1 R 71z1 —R z1z0 R-;go R zoz, From^Lemma^3,^the^first^diagonal^element^of^a,^is xx Rx1z0 Rzo Rzox, R11 From Theorem 4 the covariance matrix of the probability density function of Z1 conditioned on Zo, 0„ is given by 0„0 — Rz1z1 Rz1z0 Rziozo Rz0z1. Consequently a can be rewritten as R —Rx,z0^zozo^zoz, R-1 R R zoxi R^12-1 z1z0^zozo^ • 01/0 Application of the measurement model equation and from the definition of the covariance matrix, 01,0 can be shown to equal to 111 13;f0 HIT ±SI• Calculation similar to the proof in Lemma 3, the off-diagonal element Rxizi — Rxizo R-ziozo Rzozi is equal to P110 H,T. Finally the other off-diagonal element Rzlxi — Rz170 Rzo Rz11x1 is just the transpose of Rx1z1 — Rx1z0 R-zlozo Rzoz, and hence equal to Rzix, — Rzizo 12;104 R zox = 111 11/0 . Combining all these results a can be rewritten to be 33 PI/0 PI/0 HIT a^[^ 111 1) ^H, P11 0 H I I. T S1 j And this is the desired result. 3.2.3.5. Lemma 5 Let X1, Z1, Zo be jointly gaussian, then the probability density function of Xi conditioned on (Z 1, Zo) is gaussian. The mean and the covariance matrix of this conditional probability density function are, respectively, = X„o +Puo H1 (H1 P 111 T^) 1( Z1 H, X10),^ PI/1^PI/0^131/0 H1T (HI P110 111T + SI ) 1 (3.31) (3.32) H1 IN. P Proof The probability density function of X1 conditioned on (Z1, Zo) according to Bayes' Law is given by f(X,IZ„Zo)= f(X^ , Z1Z° , \ fkZilZ0) The probability density function f(X1,Z1IZ0) is gaussian and the mean and covariance matrix are given by Lemma 4. Also from probability theory, we can associate to this probability density function f(X„Zi Z0) a random vector. Notationally, we write the random vector as (X1lZ0,Z1lZ0). Then in terms of (X Zo,ZilZo) the probability density function f(X„ZdZo) can be rewritten as f(X1lZ0,Z1lZ0). Thus, f(XilZ0,Z,14) f(X„Z114) ^ 0,14) Hence, we can apply Theorem 4 to get the mean and covariance of f(X,1Z1,Z0) from the mean and covariance matrix of f(X, \ Zo,Z, \ Z0). 34 Firstly, from Theorem 4 the mean of the conditional probability density function of f(X11Z1,Zo) is given as X111^X110 +^H1 (H1 P„o HIT +S1) 1(z1 \ zo —Z110), where z1 \ zo is the measured value of the random vector Z1, so we will write z, \ zo as 2,. On the other hand, Z110 is the mean of the probability density function of Z1 conditioned on Zo. By Theorem 4, Z110^Z, +Rz1zo 12-zolz0(Z0 If we apply dynamic and measurement model to evaluate R10, Z110 can be rewritten as Z110 HilXi +(Fa Rxoxo HOT) R Z1020 (ZO ZO )1 • From Lemmas 2 and 3, the expression inside the curly braces is just X110. As a result, Z110^H, X„0. Consequently, this result implies X1n X„0 + P 111T (H, PI/0 HIT + SI ) 1(21 HI x110 Secondly, the conditional probability density function f(X, \ Zo,Z, \ Zo) has covariance matrix given by Lemma 4 as PI/0 [H, P PI/0 HIT P ^HiT ( R xilzo,xii, zo Rzpo,x1Izo R^, xivo,zuzo R ziizo,zilzo]) On the other hand the covariance matrix P density function f(X,14,Z1) = f((X,IZ0) (3,33) ^the conditional probability (ZilZo)) is given by Theorem 4 as Rzlizoxilzo• ^ Hence, substitution of (3.33) into (3.34) gives the desire result. 3.2.3.6. Lemma 6 The probability density function of X, conditioned on (Z1, Zo) has mean (3.34) 35 ^ (3.35) X2/1^F1 )(1, and covariance matrix ^ (3.36) P2/1^Fl P111 FIT ± Q1 Proof Firstly, let us consider (X2,Z1) conditioned on Z0. From Theorem 4, the probability density function of (X„Zi) conditioned on Z0 has mean (X2,0,Z110) and covariance matrix a, where CT = R1,1„ —R^R;1zRr oo zo^zo2). From arguments similar to those in Lemma 4, a can be rewritten as = R„ —R^ x2z0 Rz'ozo R7ox2 Rz1z0 Rzo R2ox1 — Rx2z0 Rzozo Rzaz, Rz,z, Rz1z0zoz0 R zoz,^021 (3.37) From the argument in the proof of Lemma 5, we obtain immediately a simplified expression of 022, 022 =}11 P ^+S1- Applying Theorem 4 again with (3.37), we obtain the mean, X211. and covariance matrix, P2/1, of X? on (Zi, Zo). The mean and the covariance matrix can be written symbolically as X„, = X210 + 012 0(Z, — H, E[Z, \ Z01), P ^0„ —0,2 0221 021, where E[Z, \ Zo I =H1 X110. So the proof is reduced to calculating 000(= 0T), which in turn is further reduced to calculating Rx221, Rz0z1, Rx2x2. Application of the dynamical and measurement model equations together with Lemmas 1 and 3, 0„,02„0,2(= 02,T) can be shown to be 021^Fl P1/0 H1T 1 36 011 F, P„o F,T +Q, where Q Substitution of 011,O21,012(— (3211 into P21, = 011 —012 0;1 021 with Lemma 5 results in P2„ = 0„ —0,2 02-; 021 Finally, substitution of 0„,02„0,2(= 021T) into X2/I^X210 ± 012 (322 (Z1 HI E[Z1 \ Z0 1) and applying Lemma 5 results in X,„ = F, X,„. Since X210 can be shown to be (from Lemma 2) X„0 F, This complete the proof. We can see that (3.35) in Lemma 6 are the same as (3.22) in Lemma 2 except each index is increased by 1. However, this difference is the same as observed for Lemmas 1 and 5. 3.2.4. Kalman Filter Algorithm Generalizing the indices on the formulas in Lemmas 1-6, the following KF is proposed for the dynamical and measurement model equations given by, Dynamical model Fk • Xk^7 where ERk^— °kV Qk^ (3 .3 8) Measurement model Zk Hk .Xk + Ti„ where E[rik riTk] = °kV Sk^ (3.39) The KF algorithm consists of several different pieces. To predict the state and covariance matrix of the system at time tk using measurements up to time t", we have the results (1) state prediction X k/k-1 = Fk . X k-1/k -1, (3.40) 37 (2) propagation of covariance matrix ^ 'k/k-1 = Fk Pk-1/k-1 FT Qk • (3.41) To update the state and covariance matrix of the system at time tk using measurements up to time tk , we compute the gain matrix (3) Gain Matrix Kk^Pk/k-1 ITI k[Hk 'k/k-1 HkT Sk 1, ^ (3.42) and insert it intio the state update and covariance matrix update given by (4) state update = Xkik^Kk (Zk ' k Xk/k^ (3.43) (5) estimated error covariance matrix update ^ Pk/k^Pk/k-1 Kk Hk Pk/k-1• (3.44) The KF algorithm is applied in Chapters 4 and 5 to simulated enzyme deactivation and separation of enantiomers. 3.2.5. Alternative Method of Derivation of Kalman Filte414] The Bayesian approach used to derive the KF algorithm in Section 3.2.4 has the advantage that it can be easily generalized to nonlinear filtering problems. However, the KF has the property that the covariance matrix of the error is minimized. This property is not apparent using the Bayesian approach. Therefore, we give an alternative proof in which this property is more obvious. The proof of the KF algorithm given by (3.38)-(3.44) requires the concept of innovation. We start with the definition of an innovation process. 3.2.5.1. Innovation Process The innovation un is defined by un^ = Z Hn n^ where Zn,Fin A./. are defined in Section 3.2.4. 38 3.2.5.2. Lemma 7 The process un is a gaussian process such that = 0, H 4-S.) Erunuml^ where S.,H.,P./. are defined in Section 3.2.4. Proof From the measurement model equation o= H.(X. — X.,.)+ Therefore, (3.45) E[unIZn1 = E[H. (X. — X111._1) Zn-1 + E[rinIT-1] = 0. Also from Theorem 3, un is gaussian since it is a linear combination of gaussian variables X. From the conditional expectation in probability theory, we have (3.46) E[A] = E[E[AIB]].^ Combining (3.45) and (3.46) results in qu = 0. Moreover, E[On Un = \^T E[Ifin(xn —x„,„„)+TinlIfi jxn —^f Since i is independent of Xn, we deduce immediately E[lin^] Hp. Pnin- 1 HnT Sn • Furthermore, for k=1,2,...n-1, E[On Zk T 1Zn-1] E[Un z11-1]z: From (3.46), we have, =0. 39 But, E{on^= E[u]Xk/kl= 0. Therefore, we conclude E{un 0,T J- O, for k=1,...,n-1. which completes the proof. 3.2.5.3. Theorem 5 For the dynamical and measurement model equations (3.38) and (3.39), the formulas (3.40)-(3.44) hold. Proof The results (3.40) and (3.41) follow from the application of the definition of Xkik and the dynamical equation, i.e., for (3.40) we have I Xk+1/k E[Xk+i 1Zk^E[Fk Xk^ Zk =F Xkik . A similar argument yields (3.41). Moreover, using the innovation, we can assert that Xk/k Ekk IZ„ Z2^Zk „ Uk 1. Using the fact that X„,„ is also the best linear estimate of Xk, given Z„Z„--- ,Zk_i,uk, we have the formula Xk/k^Xk/k-1 Kk Dk which corresponds to (3.43), where K, is a gain matrix to be determined. It remains to fix Kk, knowing that it minimizes the covariance of the error. The error of the estimated state can be written as = X - Xkik = (X - Xkik - Kk 0k . Hence, the covariance matrix of Ek', E[eki ek+ ^is given by 40 Pk/k Pk, , Kk(fl, P^HkT + Sk )Kk T Kk E[uk(xk —X fr/k — E[(xk — xk,k_i)UkT] Kk T However, E[Uk (Xk Xk/k 1)] E[Zk (Xk E[Hk Xk (Xk X„k I Hk 1k/k-1 • Therefore, we obtain by completing the square Pk/k Pk/k-1 Kk (Hk Plak-1 HkT + SiKkT Kk Hk Pk/1,-1 ^HkT K kT = Pk k + [Kk Pk HkT (HkPkHkT + Sk I ](HkPkHkT + ) X [KIT (HkPkHkT + Sk) 1HkPk — Pk/k-lHkT (11kPkHki. Sk 1HkPk • ^ It follows immediately that the best value of (3.47) Kk is Kk^Pk/k 1 HkT [Hk 'k/k-1 HkT Sk 1, which is (3.42) and (4.47) reduces to (4.44). Thus the proof is complete. 3.3. Extended Kalman Filter Theory The extended Kalman filter (EKF) is a generalized version of the Kalman filter when the dynamical or measurement model equations are nonlinear. The dynamical and measurement model equations are X, = fk(Xk)+ tk,^ (3.48) hk(Xk)+^ (3.49) Zk where Xk is a random n-vector, Zk is a random m-vector, 4k and Tik are sequences of white noise processes of dimension n and m, respectively, fk(Xk) is a nonlinear n-vector function, and hk(Xk) is a nonlinear m-vector function 41 It is possible to derive the EKF in a way similar to the derivation of Lemmas 1 to 6 by Taylor expansion of f„(X„), and h„(X„) about x„,k and neglecting the moments of order higher than two. However, the simplest way to derive the EKF is to linearize the dynamical model equations about X„,, and measurement model equations about Xkik_„ and apply Kalman Filter algorithms to the linearized equations. The nonlinear functions t(X„) and h„(X„) can be expanded in Taylor series about the conditional means Xki, and X„,„ „ respectively, as f(X„) = f(X„„)+F„ (X, — X„,,)+---,^ (3.50) h(X,) = h(X„,„_,)+Hk (Xk —^ (3.51) where and Hk^k/k-1 —^ OX^ Fkafk - (Xk/k ) are nxn and mxn matrices, respectively. Neglecting higher-order terms and assuming we know X„,„ and Xkik_i enables us to approximate the dynamical and measurement models as Xk±i^Fk Xk Zk (3.52) +^+ (3.53) -^Xk + Vk + lk, where u, = fic(X„,,)— F, Xkik , (3.54) Vk^hk (Xkik^Hk Xkik^. (3.55) Now if we apply the KF to the linearized dynamical and measurement model equations (3.52) and (3.53), the EKF will be obtained. The extended Kalman Filter algorithms applied to (3.48) and (3.49) is identical to (3.40)-(3.44) except that F, = Of„(X,„)/ax and Hk = ahk(xkik Nax. 3.4. More About Error Covariance Matrices To apply the EKF or the KF algorithms, derived in previous sections, we have to specify the model error covariance matrix, Qk, the initial state error 42 covariance matrix, Po, and the measurement error covariance matrix, Sk. There are no systematic methods to determine these covariance matrices. In most of the literature, the values of these covariance matrices are determined by trial-and-error based on the experience of the user. However, there is a rule to help us to decide the range of these covariance matrices. Let us start with the model error covariance matrix, Qk. Since the diagonal elements of the model error covariance matrix are the squares of the standard deviations of noises, then, the square roots of the diagonal elements of the model error covariance matrix should be large enough to account for the error between the exact model equations and the system model equations. Similarly, the square roots of the diagonal elements of the initial state error covariance matrix should be large enough to cover the difference between the exact initial values of the state and the initial values of the state specified by the EKF or the KF. Finally, the values of the measurement error covariance matrix are dependent on the instruments that we use to perform the measurements and good values of the measurement error covariance matrix can be found in the manuals for the instruments. 3.5. Exact Filter Since a nonlinear transformation can transform a gaussian variable into a nongaussian variable, the fundamental assumption of the KF theory can be violated. Some other techniques have been proposed to improve the performance of the EKF, for example, the gaussian sum technique [15] [16]. However, if a dynamical system is simple enough, then the conditional probability density function can be computed. In this section, a recursion procedure will be developed to compute the conditional probability. The algorithm is composed of two parts: 1. Evolution of the conditional probability density function. 43 2. Update of the conditional probability density function from the measured data. Consider the dynamical and measurement model equations of a system is given by (3.17) and (3.18). Lemma 8 (Evolution of conditional probability density function) Let pk(xklZk ) denote the probability density function of xk conditioned on measurements Zk. If r1 exists then the probability density function of xk+1 conditioned on measurements Zk is given by Pk+1(xk+11Zk - Pk (fk I (xk+1)1ZIC). cid( afk (xk+1 aXk+1 ) Proof This result is a consequence of Theorem 1. Lemma 9 (Update of conditional probability density function from measured data) Let pk (x.k It - I) denote the probability density function of xk conditioned on measurements Z", then Pk (Xk lZk = ^ (Zk 1Xk ) Pk (Xk iZk 1) (Zk i(P) Pk (PIZ" ) thp Proof From the definition of Zk ={z,,Z2,•••,Zk} Pk (XklZk^Pk (xk ,Zk 1) • From Bayes' rule k(zk I xk ,Zk 1) pk (xk Pk (Xk IZ = )Z" ) p,(zkiZk-1) Since the measurement noise is white noise, i.e., qv, =0 for i#j, 44 p,k(zklxk,Zk 1)=pzk(zk xk). Furthermore, PZk (zk^= fPz, (zk IT)Pk (PlZkld(P, Combining the above results, we get Pk (xk I 13,k (Zk^ j Spzk (zk pk (xklzk-J IT) pk (y1Z/C-1) chp (QED) In the following, we show how to calculate and the probability density function of rik• p,k(zkixk) from Theorem 1 Since zk -,-hk(xj+iik, then by Theorem 1 pzk (Zk^- Pylk (14 hk (Xk )). Since flk is gaussian^N(0,Sj) pzk(zkixk) ((2701/2111 Vdet(Sk)) lexp(--;[zk — hk(Xk)]T^hk (Xk )]). In particular in a one-dimensional problem, rik 1■1(0,(52 ), where a is the standard deviation, (zk - h(x3) 1 = ^ exp pz,(zkixi^ (72 1/27ta2^2^ 2\ The exact filter algorithm will be applied to the enantiomer separation problem in Chapter 5. We conclude this section with the remark that it is not feasible to compute the exact mean and variance except for one- and two-dimensional problems. 45 4. Enzyme Deactivation In a realistic experiment or bioreactor, the activity of the enzyme may decrease as time progresses. However, the concentration of substrate decreases at a slower rate than that predicted by the enzyme-substrate reaction rate equation. It is often difficult or impossible to determine exactly the model equations for the enzyme deactivation process because there are many factors involved. An example of such a process is the enzymatic hydrolysis of cellulose, cellobiose, and glucose [17], [18] cellulose —> cellobiose —> glucose Acquisition of data for the key variables and their control for the above reactions are often impeded by the lack of reliable on-line measurements of the substrate and product concentrations and the difficulties in properly modelling biological activity in processes that are time-varying. One approach to solving this problem involves the development of state estimation techniques that use available measurements to reconstruct the evolution of the state variables (i.e., the variables of the rate equation) or to estimate the parameters of the rate equation. Enzymatic deactivation in the enzymatic hydrolysis of pretreated cellulose was studied by G. Caminal et al. [15] using the extended Kalman filter (EKF) on a model by G. Caminal et al. [16]. Traditionally, the methods used most often to study thermal deactivation of enzymes require an evaulation of activity losses by incubating the enzyme without a substrate for different periods of time and then testing its activity [19]. This methodology requires that a specific set of experiments be carried out; it does not account for the possible effects of the substrate on the enzyme's stability. The power of 46 the EKF is that even if we do not know the exact equation modelling the enzyme deactivation process, we can still apply the rate equation for the substrate obtained from a quasi-steady-state approximation. All other parameters are treated as constants. If we can measure the substrate concentration, then the values of all the other parameters can be predicted by the EKF during the same time period that the measurements were taken. Application of the EKF to identify and track the slowly varying parameters of a multiple enzyme system for pretreated cellulose has the difficulty of being able to choose the system model error covariance matrix appropriately. To study the application of the EKF to an enzyme deactivation problem, a model problem is set up which simulates a single enzyme deactivation in which the parameter that represents the activity of the enzyme is slowly varying. The EKF algorithms is developed in Section 3.3. In the following, the models consist of coupled nonlinear differential equations. The EKF algorithm is applied to Euler discretizations of these models. 4.1. Reactions and Kinetic Equations To simplify the mathematics, we employ a single enzyme reaction in which a substrate A is turned into product P. A+E^C k2 >P+E where A is the substrate, E is the enzyme, P is the product, C is the intermediate enzyme-substrate complex, k_i, k 1, k2 are the rate constants. We use the same notation, i.e., A, E, P. C, for the concentrations of these quantities. The differential equations for the rate of change of concentration of 47 each element in the above chemical reaction can be derived by applying the mass-action law to obtian dA lc, A E + k_, C, dt^ dE kiAE+k ,C+Ic,C, dt^ dC k2C+k1AE—k1C, dt^ dP = k„ C . dt (4.1) If we apply the quasi-steady-state approximation to the exact equation (equivalent to setting dC/dt to zero ), the exact equations can be simplified to dA r A dt^K. + A ' dP dA dt^dt (4.2) where I-. k,, E0, K m (k2 + k_1)/k„ and Eo = total initial enzyme concentration. We simulate the enzyme deactivation by allowing rm to decrease slowly through an exponential law, i.e., rm is modelled by dr. dt ^= —6^ (4.3) where 6 is an artificial small parameter that controls the rate of deactivation. In conclusion, we simulate the enzyme deactivation process by the following system of differential equations dA _^I-. A dt^K+ A' dr ^ = —6 rm dt We will refer to (4.4) as the exact model. (4.4) 48 4.2. On-Line Estimation of Enzyme Activity Using the Extended Kalman Filter To apply the extended Kalman filter to monitor the activity of the enzyme which is governed by (4.4), we have to specify a system model which is a good approximation of (4.4), and then generate measurements which represent the observations of our artificial enzymes system. The measured quantity in this simulation is the substrate concentration, A. 4.2.1. System model In a realistic problem, it is often not possible to obtain an exact mathematical model that describes the physical phenomenon. Therefore, we have to employ a system model which is an approximation to the (real) system that generates the observations. Sometimes, such an approximation is intentional. For example, it may be desirable to employ a system model of lower dimension than the real system in order to gain computational speed and simplicity. In this study, the system model is given by dA^rt. A dt^K+ A' dr. 0 dt dK m O. dt (4.5) 4.2.2. Noiseless Measurements 4.2.2.1. Generation of Measured Data The measured values of substrate concentration, A, are obtained by solving (4.4) using a 4th-order Runge-Kutta method with initial conditions A(0)=3, rm(0)=10, Km(0)=100, and parameter c=0.02.The result is shown in Figure 1. 49 A 2 2. 75 2. 5 2. 25 2 1. 75 1. 5 2 4^6 ^ time Fig. 1 Noiseless measurements of concentration of substrate A. 4.2.2.2. Numerical Results To apply the EKF we also have to specify the error covariance matrix of the system model, Q, the error covariance matrix of the measurement model, R, and the initial state error covariance matrix, Po. Using the guidelines discussed in Chapter 3, one choice of Q, R, Po obtained by trial-and-error is 10 ^0^0^10' 0^0 Q= 0 8 x 10' 0 , R = [101, P0 = 1 0 10' 0 0^0^10'^[0^o 10' -^ - The results of running the EKF on the above choices of error covariance matrices are shown in Figures 2 and 3. Since the estimated values of A and rm are so close to the exact values, it is difficult to observe them in a graph. Therefore, the errors in the estimated values of A and rm are plotted from the exact values. 50 A_err time Fig. 2 Error in estimated values of A from exact measurement of A. A_err = estimated value of A — exact value of A. rrn_err time Fig. 3 Error in estimated values of rm. rm_err = estimated value of rm — exact value of rm. In conclusion, the behavior of rm is being tracked correctly using the EKF with suitable choices of Q, R, and Po. This tracking occurs in spite of the fact that the rm predicted by the system model is not varying with time and that the real rm is exponentially decreasing. 4.2.3. Small Noise Measurements 4.2.3.1. Generation of Measured Data 51 As in the previous case of noiseless measurements, the state that can be measured is the concentration of substrate A. However, exact measurements cannot be obtained in a realistic bioreactor. Thus, the measurements are simulated by adding white noise of zero mean and standard deviation 0.005 to the measurements in Section 4.2.2.1. (see Figure 4). A 3 2. 75 2. 5 2. 25 1. 75 1. 5 Fig. 4 Small noise measurements of concentration substrate A. 4.2.3.2. Numerical Result Using the guidelines discussed in Chapter 3, choices of model error covariance matrix Q. initial state error covariance matrix Po, and measurement error covariance matrix R are Q= 10' 0 0 0 4x102 0 0 10' [0 10 R = {0.0052}, Po = 0 0 0 10" 0 0 0 10' The results of running the EKF on the small noise measurements using the above choices of error covariance matrices are shown in Figures 5 and 6. Again the estimated A is so close to the exact values that it is difficult to observe these differences. Hence, the errors in the estimated values of A from the exact values are plotted. 52 time Fig. 5 Error in estimated A using small noise measurements. A err = estimated value of A — exact value of A. MI exact rm A -estimated rm 10 !Pi^11 3, -9- ■eVi r 2 a 4 time 5 Fig. 6 Estimated rm and exact rm in small noise measurement. As expected, the error in the estimated concentration of substrate A is quite small. The average absolute percentage error of the estimated concentration of substrate A is 0.083%. On the other hand, the average absolute percentage error of the estimated rm is 1.92%. 4.2.4. Large Noise Measurements 53 4.2.4.1. Generation of Measured Data As in Section 4.2.3.1, the measurements are simulated by adding white noise of zero mean and standard deviation 0.05 to the measurements obtained by exact solution of the exact model (see Figure 7). A 2 2. 75 ..• .• e^• .. ' 2. 5 2. 25 2 1. 75 1. 5 2 4^6 ^ time Fig. 7 Large noise measurements of concentration substrate A. 4.2.4.2. Numerical Result By trial-and-error, choices of model error covariance matrix Q, initial state error covariance matrix Po, and measurement error covariance matrix R are found to be 10-7^0^0^10-8^0^0 Q= 0 8 x 10-2^0^R = [0.052],Po^0 10 ^0 0^0^10-9^ 0^0 10' The result of running the EKF for the above choices of error covariance matrices are shown in Figures 8 and 9. For this case, the errors in the estimated concentration of substrate A are large enough to make the difference between the estimated and the exact concentration of substrate A observable. Fig. 8 Estimated A and exact A in large noise measurement. rm solid line - estimated rm dashed line - exact rm 2 4 6 — time Fig. 9 Estimated rm and exact rm in large noise measurement. The average absolute percentage error of the estimated concentration of substrate A is 0.56%. On the other hand, the average absolute percentage error of the estimated rm is 4.69%. 4.2.5. The Choice of Model Error Covariance Matrix 55 So far there is no systematic method to determine the model covariance matrix. However, experience shows that the components of the model error covariance matrix corresponding to the state variable which is governed by a model equation with larger error have to be relatively larger. In the example of enzyme deactivation, the uncertainty in the model equation arises in the state variable rm. In this section, we are going to examine the change in the estimated rm when the model error covariance matrix with a different (rm,rm) component, i.e., the second diagonal element is employed. Before doing the numerical computation, we need to introduce two quantities which indicate the accuracy of the estimated rm. Definition: The absolute percentage error, Abs%, is defined by Abs% =100 x x N n=, nnth^ estimated value — nth exact value n th exact value where N is the total number of measurements. The relative percentage error, Re1%, is defined by I^N, (nthestimated value — nth exact value y Rel% =^x 100 —x^ N^nthexact value Roughly speaking, absolute percentage error measures the amplitude of the fluctuation in the estimated values of a state variable. For example, large absolute percentage error means that the amplitude of the fluctuation in the estimation is large. Relativity percentage error can be positive or negative, and roughly speaking, it determines the best fitted curve of the estimated values. The absolute percentage error and relative percentage error of the estimated rm are calculated with the following sequence of model error covariance matrices, Qn, 56 4 x 10 Q. 0 0 0 0 covariance 0 0 10' where covariance=20 x n x 10-2. With this choice of the model error covariance matrix, we again will examine the cases of noiseless measurement, small noise measurement, and large noise measurement. 4.2.5.1. Noiseless Measurement The sequence of model error covariance matrices for the EKF is specified above and the measurement error covariance matrix, R, and initial state error covariance, Po, are given by R = [101, P0 = 10' 0 0 0 0 10' 0 0 Figures 10 and 11 show the dependence of the absolute percentage error and relative percentage error on covariance defined in Section 4.2.5. 57 Abs % error 0. 26 0. 25 0. 24 0. 22 0. 22 2^4 6^se covariance Fig. 10 Absolute percentage error of estimated rm as function of covariance. Rel^ % error 2 ^ covariance 4^6^0^10 -0. 22 -0. 22 -n. 24 -0. 25 -0. 26 Fig. 11 Relative percentage error of estimated rm as function of covariance. 4.2.5.2. Small noise measurement For the small noise measurement case, the measurement error covariance matrix, R, and initial state error covariance, Po, are given by R = [0.00521Po = 10-8 0 0 0 10-8 0 0 0 10-8 58 Figures 12 and 13 show the dependence of relative percentage error and absolute percentage error on covariance. Abs % error 14 12 10 6 4 ^ 2^4 6 10 covariance Fig. 12 Absolute percentage error of estimated rm as function of covariance. Re! % error a 2 2 4^6^2^10 covariance -0. 2 -0. 4 -0.5 -0. -1 Fig. 13 Relative percentage error of estimated rm as function of covariance. 4.2.5.3. Large Noise Measurement For the large noise measurement case, the measure error covariance matrix, R, and initial state error covariance, Po, are chosen to be 59 10' R = [0.052], P0 = 0 0 0 0 10' 0 0 10' Figures 14 and 15 show the dependence of relative percentage error and absolute percentage error on covariance. Abs % error 25 20 15 10 2 ^ 4^6^ _to covariance Fig. 14 Absolute percentage error of estimated rm as function of covariance. Rel % error 1. 5 1 0. 5 2^4 ^ 6 10 covariance -0. 5 -1 -1. 5 Fig. 15 Relative percentage error of estimated rm as function of covariance. 60 4.2.5.4. Conclusion From the above figures, we discover that as covariance increases the absolute percentage error always increases. This is indicative of the fact that a larger uncertainty in the system model leads to larger uncertainly in our estimation of rm. However, the relative percentage error decreases from positive values to negative values as covariance increases. This result shows that the value of the covariance at which the relative percentage error is zero would be a good choice. For example, for the large noise measurement case the relative percentage error of the estimated rm is zero when the covariance is about 2.25 as shown in Figure 15, and Figure 14 indicates that the absolute percentage error at this value is 17%. To reduce the absolute percentage error of the estimated rm, we can adopt a more sophisticated filter algorthim such as the Jestimator invented by Jazwinski in his 1974 paper (see [20]). However, a simpler way is to sacrifice the relative percentage error. From Figure 14, the covariance decreases to about 1 when the relative percentage error of the estimated rm increases to 1%, and the absolute percentage error of the estimated rm drops to about 11%. Unfortunately, this method cannot be applied to a realistic enzyme system because the exact values of certain state variables are not available. However, this calculation shows that a good choice of the model error covariance matrix exists. 4.3. Iteration of the Extended Kalman Filter In an off-line analysis, we can apply the EKF again to the estimated state obtained by the application of the EKF to a set of measurements. Unfortunately, if the model error covariance matrix for the second application of the EKF is chosen to be the same as the first EKF, then the relative percentage error and the absolute percentage error of rm will become larger 61 than that obtained from the first EKF. On the other hand, if the model error covariance matrix of the second EKF is chosen to smaller than the model error covariance matrix of the first EKF, then the relatively percentage error of rm will become larger. As a result, it is not practical to iterate applications of the EKF to the output obtained from the EKF. However, a regression analysis can be done on the set of estimated values of rm which is obtained from the first EKF. For example, if rm is given in terms of a quadratic function of time, e.g., = a +13 t +7 t2, where a, 13, 7 are constants determined by the estimated values of rm, then the modified model equation for rm is dr ^ dtm =13+7t. Hence, the modified system model is dA^rm A dt^Km +A dr m^t, dt dKm =0 dt Now we can apply the EKF to this modified system model. The following schematic diagram summarizes the iteration procedure. 62 system model Extended Kalman Filter measurements estimated state d Extended Kalman Filter ^ modified system model measurements estimated state 4.3.1. Numerical example We apply the above iteration procedure to the case of the large noise measurements since the estimated rm has both larger absolute percentage error and relative percentage error. Firstly, an estimation by the EKF based on the system model (4.5) is performed. Having obtained a set of estimated values of rm, a linear regression analysis is performed and the best straight line fit is calculated to be = 10.1388— 0.1799 t Hence, the model equation of rm is changed to dr. ^ = 0.1799. dt Consequently, the system model (4.5) is changed to: 63 dA^I-. A dt^K + A dr. ^ =-0.1799. dt dK ^ =0. dt (4.6) Since the system model has been changed, we have to modify the model error covariance matrix. The modified model error covariance matrix. Q, is 4 x 10-5 0 0 80 x 10-4 0 0 0 10-9 0 The absolute percentage error and the relative percentage error of the estimated values of rm based on system model (4.6) are improved tremendously. Figure 16 is a plot of the estimated values of rm for the original system model (4.5) and the modified system model (4.6). The following table summarizes the performance of the first and second iterations of the EKF in terms of the absolute percentage error and the relative percentage error. 64 rm dashed --exact rm solid black ^ estimate rm with modified system model grey --estimated rm with raw system model time 2 4 6 Fig. 16 The estimated values of rm for the original system model (4.5) and the modified system model (4.6). Table 1 Absolute and relative percentage errors of estimated A and rm. 1st Iteration Absolute Percentage Error Relative Percentage Error A 0.6954% 0.3162% 0.5034% 0.2711% I rm 5.8456% 1.9369% 0.9008% 0.3168% 2nd Iteration Absolute Relative Percentage Error Percentage Error 65 5. Separation of Enantiomers If a molecule is not superimposable on its mirror image, then it is chiral. These images are called enantiomers or chiral. Living organisms often show different biological responses to one of a pair of enantiomers in drugs, pesticides or waste compounds. For example, sodium L-(+)glutamate tastes deliciously in food, whereas its mirror image D-(-)glutamate tastes bitter or flat. Furthermore, there are a number of cases where one isomer is largely responsible for the desirable biological action while the other predominates in the toxicity/undesirable action. Therefore, industry develops methods of separating enantiomers. In this section, we are mainly concerned with the enzymatic technique of separation and the application of the EKF to monitor the process of separation of enantiomers. 5.1. Enzyme-Substrate Reaction Rate Equation Consider the following coupled enzyme-substrate chemical reaction A + E^Ci k2 >E + P, ki B +E^C, <^ k2 >E + P2 where A and B are substrates, E is the enzyme, 131 and 132 are the products, C1 and C, are the intermediate enzyme-substrate complexes, k_1, kl, k,, k1, k'1, k'2 are the rate constants. If we apply quasi-steady-state assumption, the following rate equation is obtained [3] ( cf. Chapter 2) 66 dA_ ^VniA K.B A dt^K.„ A + K.A B + dB^Vni„ KT, B dt^K.B A + K.„ B +^KmB (5.1) where A and B are the concentrations of the substrates A and B, respectively, and the parameters of the first-order system differential equations are defined by K= +k k ' mA^k, 2. K^1(11 InA^k; VifiA = Eo k„ Vm], =E0k. where E0 is the initial concentration of the enzyme. A numerical solution of the coupled first-order differential equation obtained from using the 4th-order Runge-Kutta scheme is plotted in Figure 17 for the following initial conditions and parameter values Ao = 0.05, Bo = 0.05, VinA = 1 00 V.B = 5, K niA = 0.001, and Kn. =0.0025. 67 concentration O. 05 0. 04 0. 02 0. 02 0. 01 ^ time 0. 0002^0. 0004^0. 0006^0. 0009^0. 001^0. 0012 Fig. 17 Numerical solution of the rate equation for the enantiomer. 5.2. Enantiomeric Ratio If the rate equation of substrate A is divided by the rate equation of substrate B, then the denominator is canceled. A first-order ordinary differential of A as a function of B is obtained _ V„,A Kn.,B A dB VmB^B • dA (5.2) An exact soultion is given by VmA K InA — ink =^ ( lnB—lnBc), VinB KmA (5.3) where Ao and Bo are the values of A and B at time equal to zero. The solution can be rearranged as In( A E^( °^ B In ^ Bo VmA (5.4) 68 where E is called the enantiomeric ratio and is a constant made up of the parameters in the rate equations. The enantiomeric ratio, E, is a measure of the relative rates of the two completing reactions. The fact that E=constant makes the enantiomeric ratio useful for comparing the effectiveness of different enzyme systems for separating substrate A and B. Moreover, A can be solved in terms of B, E, Ao, and Bo as, B E (5.5) A = Ao(— Bo and B can be solved in terms of A, E, Ao, and Bo as 1 B- T., BO A E —A 0) (5.6) 5.3. Conversion The conversion, C, is a measure of the fraction of how much substrate is left in a chemical reaction and is defined by C1 = A+B Ao +Bo (5.7) If we know the enantiomeric ratio and the conversion, we can predict the concentrations of A and B by solving the following simultaneous algebraic equations A +B =(1-C)(Ao +BO, (5.8) A Ao However, in general, since the second equation has the exponent E, an analytical solution cannot be obtained. 5.4. Perturbation Analysis 69 The graph of the solution of the enzyme reaction equation indicates that the decrease in concentration of substrate B is relativity small compare to the decrease in concentration of substrate A. As a result, the concentration of substrate B can be considered constant to a leading-order approximation. We show this by application of a perturbation method. First, we rescale the enzyme-substrate reaction rate equation and express it in dimensionless form. Define the following dimensionless variables (5.9) If the rate equations are expressed in terms of these dimensionless variables, they become cfA (it^Ä+11+1' dfi= ch^A+B+1' 1 where E = - . (5.10) The lower limit of the enantiomeric ratio is usually between 25 to 50, so c is between 0.02 to 0.04. Therefore, we can apply a perturbation method to (5.10). The leading-order approximation is obtained by setting E equal to 0 and solving the resulting system of first-order differential equations with suitable initial conditions. The leading-order equations become 70 dA A++1^ (5.11) (IT The initial condition^at t = Tic is A = A, B = Bk . The leading order approximation is = f3, , (5.12) (A— Ak) + +^inkA = —(c — Tic). Ak Instead of setting B = ilk, we used (5.6), then the following solution is obtained A ) (A _ A ( \ (5.13) ± Ei in ^ = —(t — tk). In Figures 18 and 19 are plotted A and B as functions of time for both the perturbation solution (5.13) in terms of the original variables A, B. t, and the numerical solution obtained from the 4th-order Runge-Kutta scheme with the initial condition t=0, A0=0.05, B0=0.05. Note that there is excellent agreement between the perturbation and numerical solutions. 71 Fig. 18 Exact and perturbation solutions of A. Fig. 19 Exact and perturbation solutions of B. 5.5. Construction of Difference Equation From Leading-Order Perturbation Solution Having obtained the leading-order perturbation solution of the dimensionless enzyme-substrate rate equations, we can employ them to construct the corresponding difference equations, namely the forward and backward propagating equations. The importance of the forward propagating equations 72 is that they are discrete equations and we can use them to describe the dynamics of the state variables in the EKF algorithms. The importance of the backward propagating equations is that they can be used to derive an exact filter as shown in [12] (see Chapter 2 as well). 5.5.1. Forward Propagating Equations We apply (5.13) to find out the values of -A and 11 at^tk, . The result is a set of coupled difference equations, hic+1 Ak (5.14) (Ak, - Ak) + (1+ hiln A-k+1^—(tk+i tk)Ak ) If tk lc I+1 is small then ikk±i fikk and Ek,iihk will be close to but less than 1. Hence, both (1— ÄkFirikk) and (1— iLiii3k) are small and positive. Since both 1n(Ak±ii-Ak) and 1414k±irBk) can be rewritten as ln(1 — (1 — Ak±lAk)) and 111(1 - (1 - hk+iiiik)) , respectively, the first-order Taylor expansions of them are given by in(AkjAk) (-I A, +,/A in(hk+iihk) (-1 + hk+1/k) From the above two equations, (5.14) is reduced to Ak+1 = Ak eCk+1^Ak Ak ijk +1 (5.15) tk) hk Bk, =^6 erk+1 Ak Bk +1 Note that (5.15), is equivalent to Euler discretization of (5.10), and is called the forward propagating equations. 73 When (T1,1— -c,) is small, the numerical solutions of the above difference equations are very close to the numerical solutions obtained using the 4thorder Runge-Kutta scheme (see Figures 20 and 21). Fig. 20 Numerical solutions of substrate A from the forward propagating equations and the 4th-order Runge-Kutta scheme. Fig. 21 Numerical solutions of substrate B from the forward propagating equations and the 4th-order Runge-Kutta scheme. 5.5.2. Backward Propagating Equations 74 A, andB k+lare expressed in terms of . Ak and -Ek in the forward propagating equations On the other hand, equations expressing A, and Bk in terms of it„ and ifl„ are the backward propagating equations. To find the backward propagating equations, we can treat the forward propagating equations as simultaneous algebraic equations and consider A, and fi„ as the unknowns. However, the algebraic equations may not have an analytic solution if it is non-linear in A, and fi,, so we have to find an approximate solution of the simultaneous algebraic equations. A general method is discussed in Appendix 2, but in this section, we will employ the perturbation solution of the dimensionless enzyme-substrate rate equation to derive the backward propagating equations. It turns out that the equations obtained using these two approaches are the same. We apply (5.13) to determine the values of A and 14 at '1^and rename the subscripts k and (k-1) as (k+1) and k, respectively. Then the following system of coupled difference equations is obtained Bk+l = fi A.-,z 1) k ( : 6,^ (5.16) (Ak — A„)+ (1+ ii„,)14^ I = A, If tki-i tk -(tk tk+, is small then ikk ritk, and ik ii3-„ will be close but large than 1. Hence, ( Ak/Ak+i 1) and (Bk/Bk^l-1) are small and positive. Since both ln(i1k rk„) and In(fik/f3,,) can be written as MO +(AkiAk —1)) and + (11„M„ —1)), respectively, the first-order Taylor expansions of them are given by 1n(A„/A„) —1+ Ak/A„, 75 From the above two equations, (5.16) is reduced to Ak+.1 ^tk+1)Akk+1 Ak4 hk = ttk-fi 13,„ + 1 (5.17) (t^)i36 -"2"+1 • A, + B„ +1 Equation (5.17) is called the backward propagating equations. To determine the initial conditions for the backward propagating equations, we solve the forward propagating equations (the results were shown in Figures 20 and 21). This determines the values of À and 11 at the initial point for the backward propagating equations. The result is shown in Figures 22 and 23, A 0.06 forward propagating 0.04 0.03 backwar propagating 0.02 ^ time 0.0002 0.0004 0.0006 0.0008 0.001 0.0012 Fig. 22 Numerical solutions of the forward and backward propagating equations for substrate A. 76 I^ ^ 0.0004 0.0006 0.0008 0.0495 0.049 0.0485 ackward propagating thne 0.001 0.0012 orward propagating 0.048 0.0475 0.047 Fig. 23 Numerical solutions of the forward and backward propagating equations for substrate B. We observe a small gap at the beginning of Figures 22 and 23, which occur because the backward propagating equations are not the exact inverses of the forward propagating equations. The error obtained in the general theory in Appendix 2 is 0(At2). 5.6. On-Line Estimation of Substrate Concentration by Extended Kalman Filter 5.6.1. Defining the problem In a realistic chemical reaction separating enantiomers, we have to determine the individual concentration of the enantiomers so that we can monitor the relativity purity of the product and decide when to stop the reaction. However, it is very difficult to measure the concentration of individual substrates, and it is relatively easier and much cheaper to measure the total concentration of the enantiomers. Thus, the problem is to estimate the individual concentration of each enantiomer from the sum of concentrations of the enantiomers. To study the problem in more detail, we simulate measured 77 data and apply the EKF technique to filter the measurement noise and estimate the concentration of each substrate. 5.6.2. Generation of Measured data As discussed in last section, the quantity that can be measured easily is the sum of concentrations of substrate A and B. For the purpose of this study, we generate the measured data by solving the enzyme-substrate rate equation (5.1) using the 4th-order Runge-Kutta scheme with initial condition Ao = 0.05, Bo = 0.05 and parameters value \CA = 100, K 0.001, V., = 5 K.„ = 0.0025. (With this set of parameter values, the enantiomeric ratio, E, is equal to 50.) Since a set of exact measurements cannot be obtained in a realistic reactor, we simulate the noise with white noise. The white noise assumption is only an approximation; however, the Central Limit Theorem states that the limit distribution of a sequence of random process tends to a gaussian distribution, which motivates the gaussian distribution in the white noise assumption. The gaussian distribution employed in this simulation has zero mean and standard deviation 0.001. 5.6.3. Numerical Results To apply the EKF, we have to specify the dynamical model of the state variables, the measurement model, the error covariance matrices of the system model, and the measurement model. The system model is given by (5.1). The measurement model is determined by the algebraic relationship between state variables and the quantities which we measure. Mathematically, the measurement model for this simulation is S=A+B where S is the sum of the concentrations of the substrate enantiomers A and B. 78 So far there is no systematic method to determine the error covariance matrix of the dynamic model and the measurement model, all we can say is that the larger the error in the model the larger the entries in the covariance matrix should be. The determination of a good choices depends on experience and using a trial-and-error method. Choices of the model error covariance matrix. Q, the measurement error covariance, R, and the initial state error covariance matrix, P0, is Q= [1.610-9^0 10 -9^0 R = [0.0012],^= [ 8.510-ill^ 0 ^ 0^10-12]• The result of running the EKF, for the above choices of covariance matrices, are shown in Figures 24, 25, and 26. Enantiomeric Ratio et*i. 50.6 : 50.4 50 2 0.1^0.2^0.3^0.4^0.5 Conversion Fig. 24 Enantiomeric ratio calculated from the estimated concentrations of substrates A and B. 79 A 0.05 0 .04 0.03 0.02 0.01 0.0002 0.0004 0.0006 0.0008 0.001 0.0012 time Fig. 25 Exact and estimated concentrations of substrate A. 0.0495 0.049 --At-A.09..20.00040.00060.0008 0.0010.0012 time • 0.0485 0.048 0.0475 0.047 Fig. 26 Exact and estimated concentrations of substrate B. 5.6.4. Filter Divergence The estimated values of the substrate concentrations are extremely good compared to the exact values of substrate concentrations obtained by solving the rate equation by the 4th-order Runge-Kutta scheme. However, the statistics obtained from the EKF are inconsistent with the actual statistics when the residues (the residue of a state variable is the difference between the exact value and the estimated value of the state variable) of A and B are plotted. 80 residue and s.d. of A residue of A 0.0001 0.00005 60020.00040.00060.ms o. 111-19-reol2 time —0.00005 —0.0001 standard deviation predicted from EKF Fig. 27 Residue of A and predicted standard deviation from the EKF. residue and s.d. of B 0.00003 0.00002 0.00001 °:1112'2 -"I:1'7M 4 0.0JTtl00t1.001 0 . 0012time —0.00001 residue of B —0.00002 —0.00003 predicted standard deviation from EKF Fig. 28 Residue of B and predicted standard deviation from the EKF. In fact, as shown in Figure 27, a significant portion of the residue of substrate A lie outside the two curves representing the standard deviation. This fact causes the gain matrix of the filter to become too small to sufficiently amplify the subsequent data to correct the values of the state predicted by the system 81 model. If the system model is not a good approximation, it could lead to a divergence between the subsequently estimated state and the measurement. In the next section, we will show how to correct this problem. 5.6.5. Correction of Filter Divergence A number of approaches to deal with the divergence problem have been suggested. In fact, these approaches are useful when there is an inconsistency between the actual statistics and the predicted statistics from the EKF. One of the simplest approaches is to increase the model error covariance (for a motivation of this method, see the next section). Roughly speaking, the increased model error covariance would hopefully keep the gain matrix away from zero and the subsequent measured data will not be ignored. Therefore, the model error covariance matrix, Q, is increased to [4 -10-2 0 0 9.10' and we recalcuate the estimated concentrations of the substrates A and B with this new model error covariance matrix. Then we calculate the enantiomeric ratio and the conversion. The results are shown in Figures 29 to 33. The statistics of substrate A has been improved as shown in Figure 31, but, unfortunately, the accuracy of the enantiomeric ratio has to be sacrificed. However, the changes in the estimated values of the concentration of substrate A and B are so small that the changes are not observed. The reason the enantiomeric ratio has been changed is that E is dependent on 1nB0ilnB and a small change in B will produce a large change in E. 82 0.05 A 0.04 0.03 0.02 0.01 0.1 0.2^0.3 0.4^0.5 conversion Fig. 29 Exact and estimated concentrations of substrate A. 0.3^0.4 0.5 conversion 0.0495 0.049 0.0485 0.048 0.0475 0.047 0.0465 Fig. 30 Exact and estimated concentrations of substrate B. 83 residue of A residue and s.d. of A standard deviation of A 0.0006 0.0004 • 0.0002 .D.9.4020.00040. pan 0.0008 10.001 0.0012 time —0.0002 —0.0004 -1^• Fig. 31 Residue of substrate A and predicted standard deviation from the EKF. residue of B residue and s.d. of B standard deviation 0.0002 0 0001 0.00020.00040.%fkgN..0008 1:061 0.0012 time —0. 0001 —0.0002 Fig. 32 Residue of substrate B and predicted standard deviation from the EKF. 84 Enantiomeric ratio 52 51.5 51 50.5 ......^•... 0.1^0.2 0.3^0. 4 0.5 conversion Fig. 33 Estimated enantiomeric ratio vs. conversion. 5.7. An Example Of Filter Divergence The fundamental problem in filter divergence is that the dynamics of the state are subject to error and if the EKF is constructed on the basis of an erroneous model and the system noise which is much smaller than the model noise, then the gain matrix of the filter will become very small when it operates on many data and the subsequent measurements will be ignored. As a result, the filter will follow the wrong state predicted by the dynamics of the system. We present an example from Schlee et al. [21] to illustrate this. A more fundamental illustration of the divergence problem is given by Fitzgerald [22]. 5.7.1. Example Consider the problem of estimating the altitude h of a space vehicle from altimeter data. Suppose that the filter is based on the assumption that the altitude of the space vehicle is a constant, whereas, in fact the space vehicle is climbing at a constant speed. Define the following notation y: measurement at instant tk. hk: actual altitude of space vehicle at instant tk. 85 hk: the altitude variable of the space vehicle in the dynamical model of the EKF at instant Pk/k: tk the update estimated error covariance matrix at instant The measurement model is yk hk + vk , where vk is the white noise at instant tk. The filter dynamical model is hk+i^hk• The actual dynamical model is h, = hk +S. Substituting the above equation into the EKF, a system of difference equations for kik and Pk/k are obtained Pk-14-1^( fikik^ p^ + R k-l/k-1 Yk fik-1/k-1)) (Pk-1/k-1 )2 Pk/k^k-lik-1 = P (5.18) (5.19) k-lik-i + R. Equation (5.19) has the solution PoR PkIk k Po + R (5.20) • Substituting (5.20) into (5.18) yields ii k+1/k^ 11 kik^P° (k^1)P0 + R (Yk kik ). (5.21) The actual estimation error, Fik, is defined by ' hk^hkik • Thus, the difference equation for the actual estimation error is kP+R fik-11/k^(kik ± ' +1)P0 ° (k +R) P0 yk+1 (k + OP0 + R (5.22) 86 The solution of (5.22) is kik Rk — Ok/2]Po + k R^Po^k Ey s^ kPo+R^kPo +R^kPo +R (5.23) R = ^ 110/0 The first term of (5.23) goes to zero as Ic-->00. Furthermore, the last term of (5.23)^goes to^zero^in 1.i.m.xn = x^if limE[Ix„ — n__00 the^mean square sense, i.e., = 0. However,^the second term of (5.23) involving the rate s is unbounded in k and goes to infinity as (k-1)/2. As a result, the estimation error is unbounded and the filter is diverging. Finally, if the speed s is zero then the middle term vanishes and the filter becomes stable. 5.8. Estimation by Exact Filter Algorithms We now apply the exact filter algorithms developed in Chapter 3 to estimate the substrate concentration from the rate equation for the enantiomer in one dimension. This will be compared to the results obtained by the extended Kalman filter when there is error in the initial value of the concentrations of substrates A and B. 5.8.1. Dynamical and Measurement Models in One Dimension Substituting (5.5) into the rate equation (5.1) reduces the system of two coupled first-order ordinary differential equations into a first-order ordinary differential equation dB ^Vm„ Km,„ B dt^K ^+Km, B +Km, Km, • (5.24) To apply the discrete EKF, we have to discretize the rate equation (5.24). If Euler's discretization scheme is employed, the following difference equation is obtained Vrni, K Bk Etic+i = Bic At ^ Km, Ao(Bk/Bo)E +Km, B, + KmA Km, • (5.25) 87 However, according to Lemma 8 in Section 3.4, we need to solve Bk in terms of Bk+i in order to calculate the evolving conditional probability density. A leading-order approximation of the solution is given by the backward propagating equation Bk = Bk±i ± VrnB K At ^ Kinn Ao(Bk±i /Bo )E^Bk±i K ^• (5.26) The measurement model is Sk = Ak Bic + v„ where vk is white noise with standard deviation a. Combining the above equation with (5.5) results in a measurement equation with one variable Sk = Ao ^ Bo Bk , +B k V k • From Theorem 1 in Section 3.1 and the above equation, the propability density of Sk conditioned on Bk is given by p(SklBk) = )2 cs2^ 1 (S^Ao (Bk /B0 — Bk where Ao, Bo are the initial concentrations of the substrates A and B, and a is the standard deviation of the noise of the measurements. 5.8.2. Numerical Results We employ the dynamical and measurement models derived in the last section and apply the exact filter algorithms (Lemmas 8 and 9 in Section 3.4) recursively to calculate the conditional probability density at the times when the measurements are taken. The estimated state of the dynamical model is the 88 mean of the conditional probability density which can be calculated by numerical intergation. In conclusion, if there is an error in the initial value of the state, the estimated state obtained using the exact filter will converge to the exact curve of the state at a faster rate than the estimated state obtained using the extended Kalman filter. The following choice of parameters were chosen for the exact filter Initial state (A0,B0) = (0.0502,0.0502). Initial state error covariance matrix =0.01. The covariance matrix of noise-0.0012. The estimated concentractions of substrate B obtained by the exact filter and the EKF are shown in Figure 34. Fig. 34 The estimated concentractions of substrate B obtained by the exact filter and the EKF. 89 6. Discussion and Future Research 6.1. Discussion The enzyme kinetic rate equations are developed in Chapter 2 with special emphasis on the quasi-steady-state assumption. We demonstrate that the kinetic rate equations obtained using the quasi-steady-state assumption are the leading-order equations of a system of coupled singularly perturbed first-order differential equations. Thus, the kinetic equations contain error terms that may give some guidelines for choosing the order of magnitude of the model error covariance matrix in a real enzyme system. However, the method using the differential equation is not an efficient way to get the quasi-steady-state rate equation for a complicated multi-enzyme system. Another method is the King Altman algorithm (see [1]). One common feature of the two simulated enzyme systems in Chapters 4 and 5 is that the evolving states of a partially observable system are reconstructed using the EKF algorithm. Estimations of the states of the enzyme system without and with proper model equations are studied in Chapter 4 and 5, respectively. In both cases, different amounts of noise were added to the measurements. The accuracy of the estimation using the EKF algorithm decreases when the noise in the measurements increases. There is a lower bound in the error of the estimated states that cannot be further reduced by changing the parameters of the EKF, e.g., the model error covariance matrix, the measurement error covariance matrix, and the initial value error covariance matrix. So far there is no systematic method to determine the values of these covariance matrices. In most of the literature, the values of these covariance matrices are determined by trial-and-error based on the experience of the user. However, there is a rule 90 to help us decide the ranges of these covariance matrices. The square roots of the diagonal elements of the model error covariance matrix should be large enough to cover the error between the exact model equations and the system model equations. Also the square roots of the diagonal elements of the initial state error covariance matrix should be large enough to cover the difference between the exact initial values of the state and the initial values of the state specified in the EKF or the KF. The largest error in the estimated state of the simulated enzyme deactivation system in Chapter 4 is about 6% (cf. Table 1). In the off-line estimation, this error can be reduced to 0.9% by iteration on the dynamical model equations (cf. Chapter 4). To drop this error from 6% to 0.9%, only one iteration is needed. However, it is difficult to generalize this method to on-line estimation. To reduce the error in an on-line estimation, we need to employ more sophisticated techniques such as the Lestimator developed by Jazwinski in 1974 [20]. In the simulated enzyme system in Chapter 5, we found that the statistics obtained from the EKF are inconsistent with the actual statistics obtained when the residues of the states are plotted. This phenomenon is known as "filter divergence". The explanation most often offered for this phenomenon is that the estimated error covariance matrix becomes unrealistically small, so that too much emphasis is put on the value of the current estimated state and subsequent measurements are effectively ignored. Fortunately, this divergence in our simulated system can be corrected by increasing the values of the entries of the model covariance matrix. However, in general, it is not true that the divergence can always be removed by increasing the intensity of the process noise assumed in the filtering model. Indeed, a more fundamental illustration of the divergence problem is given by 91 Fitzgerald [21]. He suggests that there are two types of divergence, true divergence and apparent divergence. True divergence occurs when the meansquare error can actually be shown to approach infinity with increasing time. Apparent divergence occurs when a steady state is reached but the associated errors are too large to allow the estimated state to be useful. The divergence in our problem is apparent divergence. Indeed, the apparent divergence can be removed by more sophisticated algorithms, for example, by the limited memory optimal filter [23] which computes the best estimate based on a "moving window" over the most recent data. Finally, an exact filter algorithm is used to compute the conditional probability density function of the substrate concentration in Chapter 5. The mean of this probability density functions, computed by numerical integration, is compared to the values of the estimated concentration obtained from the EKF when there is error in the initial values of the concentrations of the substrates. We find that the exact filter has a faster rate of convergence than the EKF. However, the computational time taken by the exact filter is much longer. Thus, it is not feasible to apply the exact filter algorithms in on-line estimation. 6.2. Future Research Central to the estimation of the state variables of a dynamical system is the determination of the probability density function of the state variables conditioned on the set of measurements. If the dynamical system is linear, the state variables that are gaussian initially will evolve to a set of gaussian random variables. To determine the probability density function of a gaussian random vector, it is sufficient to determine their mean and covariance matrix. The Kalman filter gives an exact solution for the linear system. As shown in the derivation of the EKF algorithms in Chapter 3, the dynamical equations are 92 linearized at the estimated state variables. This will introduce a significant error when the intervals between measurements are large. In the simulated enzyme system in Chapters 4 and 5, the intervals between the measurements are chosen to be sufficiently small to minimize this error. This requirement can be realized for the measurements taken by a radar, but may be not be realized for some enzyme systems. Although the probability density function conditioned on the set of measurements evolving in time can be described with a system of differential or difference equations, an exact solution is difficult to obtain. Moreover, the state variables in a nonlinear estimation problem are non-gaussian. The nongaussian conditional probability density function cannot be characterized by a finite set of parameters, e.g., mean and moments. Thus, the unknowns of the exact nonlinear filter are infinite (cf. the unknowns in the linear Kalman filter are the means and the first-order moments). Therefore, it is natural to investigate the approximation schemes for the conditional probability density function. Motivated by Fourier series, an approximation scheme based on expanding the conditional probability density function using orthogonal series (e.g., the Edgeworth series [24]) is proposed [25]. However, this approach has a distinct disadvantage in that, when truncated, the resulting series is not a valid density function. Therefore, it is necessary to keep a large number of terms in the series to minimize the error. This will make the scheme computationally unattractive. However, there are other expansion schemes to overcome this difficulty. One of them is the gaussian sum approximation scheme proposed by Sorenson and Alspach [15], [16]. Roughly speaking, this approximation scheme is to expand the conditional probability density function by sequences of weighted gaussian distributions of different mean and covariance. The resulting filter equations are a set of extended Kalman filters, 93 with each filter tracking the evolving mean and covariance matrix of its assigned gaussian density. In a low noise set of measurements, the gaussian sum filter can be very nearly optimal. In a high noise set of measurements, we have to frequently reinitialize the convariance matrices of the individual EKF with a smaller convariance [13]. Apart from the probabilistic approach discussed in this thesis, there is another technique using neural networks to estimate and predict the future state of a dynamical system. In the probabilistic approach, a dynamical model for the system process is essential. On the other hand, in certain circumstances, even a first-order approximation to the exact equations of a process is not known, e.g., the value of a stock in the stock market. In this situation neural networks are an excellent tool to learn relationships without requiring knowledge of the dynamical equations and then to predict future states. The theory of neural networks was developed back in the 1950s and has become one of the hottest topics in artificial intelligence over the past 10 years. A network operates in analogy to a human brain. The storage and processing of the information are facilitated by modifying the synaptic strengths of the junctions between interacting neurons. Neural network technology has the greatest potential in areas such as speech and image processing, but applications to areas such as process control, dynamic modeling, and robotics are currently being developed. Indeed, Lapedes and Farver [26] studied the ability of layered feed-forward networks to forecast time series, which is an important problem in economics, meteorology, and many other areas. The advantage of the neural network technology is that, no information about the dynamical equatins is required. All that is required is to 94 design the architecture of the neural network and to feed the neural network model the necessary information to let it learn by representative examples. In off-line estimation of a dynamical system that is partially observable and has no proper dynamical equations (e.g., the simulated enzyme deactivation system in Chapter 4), we can use a neural network to process the estimated state obtained using the EKF, and then develop the dynamical equations for the system. However, more research still has to be done on developing a neural network type on-line estimation scheme in analogy to the EKF. 95 Appendix 1 Matching Slope Technique for a Singularly Perturbed First-Order System In this appendix, we will prove that the leading-order terms in the rate of change of the substrate concentration for the single enzyme reaction described in Chapter 2 is always negative, We determine the initial values for the outer equations of a singularly perturbed first-order system by matching the slopes of the inner and outer solutions when the inner solutions are expressed in terms of the outer variable which in turn is set equal to zero. The matching method will be applied to find the initial condition for the outer equation of the singularly perturbed first-order dimensionless enzyme system derived in the Heineken et al. 1967 paper [2]. The advantage of this method is that the explicit solutions of the outer equations are not required. This method is especially convenient when the outer solution is an implicit function of the outer independent variable, where it is difficult to apply the method of matched asymptotic expansions. Inner and Outer Equations In the following we summarize the inner and outer equations in the Heineken et al. paper. Inner equations are given by dy = -6 [y - (y^- X) z],^ dc (A.1.1a) dz =y—(y+Oz. da (A.1.1b) The perturbation expansions are assumed to take the form ^ (A.2.2a) z(a) = zo(u)+ c z,(a)+ c 2 z, (a) F •^ (A.2.2b) Y(G) = Y0(175)+ 6 Yi (CS) + 62 y2 (G)+-, - - - The corresponding initial conditions are given by 96 yo (0) = 1, zo (a) = 0, y.(0) 7-- 0, z.(0) = 0, for n=1,2,3,.... The leading-order equations are obtained by substituting (2) into (1) and setting 6 =0, thus dyo 0 da (A.1.4a) dz ° =y0 —(3'0^zo. da (A.1.4b) The equation for the first-order and second-order terms are given by dy, = +(p.— X); + yo zo,^ da^ ZI (A.1.5a) =^+(p^+ yo z, +^zo,^ (A.1.5b) dy, =^+(p.— X) z, + y, zo + yo z„^ da (A.1.6a) dz, da (A.1.6b) da = y2 -y2 ZO Y1 Z1 -y0 z., ^ Outer equation For the outer region, we rescale time by t=-66 and let y(r)=y(c) and (-r)•= z(r7), thus dt =—Y“Y-1-11-20i, _ 6^ y dt _ _ —( y + p) z. (A.1.7a) (A. 1 7b The perturbation expansion for the outer solution is given by y(t) = yo(T)+ 6 yl(T) +62 y2(T)+--, (A.1.8a) z(T) = zo(T) + c Zi (t) ± (A.1.8b) 62 Z2 (t)-1-• • The equations for the leading-order and first-order terms in the solution are given by ^ 97 y0 + (y0 + i - x.) z0, dY°^ at 0= Yo -670 + 2,) + ^Y,^— x) dY1^ crio _ (IT yi yo z,^Y, zo^zi•^ (A.1.10b) Computation of the Initial Conditions The leading-order and first-order inner solutions obtained by Heineken et al. [2] are Yo(a) = 1,^ 1 — exp(—(1+ u) zo(a)= ^ 1+ u Yi(a)= (A.1.11a) (A.1. 11b ) Xa^(1+11—X) ^,^ [1 exp(—^(Y)] (1+ P)^(1+ 1-1)- L (A.1.12a) (1+ — 22) (1+u-2X) 2((l+ ) + )^0+03^(1+1)4^0+04 exp[—a)] zi ,e^2■, u^u ex_((1 + t) )] (1+04 {^(1+ 04 (52 + (1+ )(1+ —(-1)(3' + (1+2u—X-2Xu+ull.^(A.1.12b) Therefore, the inner expansion written in terms of the outer independent variable t is (A.1.13) -11^X^)+E (I +11-X) +0^ 6/^0+^0 +02 From (A.1.9), we get dyo^x yo dt^y0 + t (A.1.14) ^ ^Yo ^Yo 98 so that as T--->0 = _Y°(°,) x . Yo (0) + This relationship between the slope and the value of the function at T=0 is valid since this corresponds to an interior point in the domain. We require the above slope to be equal to the leading-order slope of the inner solution expressed in terms of the outer variable "C, as T—>0 (0)^X n (A.1.16) P^+ P, The unique solution is given by yo(o) = 1 From (A.1.9b) the initial value of Zo is (A.1 .1 7) io (0) = 1+14^ Therefore, the initial condition for the outer leading-order expansion is 37-0(0) =1. For this initial condition, we see that dyoidt o. Combining this result with the inner solution, we conclude that the leading-order expansion for the rate of change of the substrate of the single enzyme reaction is negative. Furthermore, if the perturbation series expansion is uniformly convergent, then the rate of change of substrate is always negative. Next, we calculate the initial condition for the first-order outer solutions. If we differentiate (A.1.9b) w.r.t. 0 = Ner) —^N(r) T, we get (t) — (t) yo(t) —^(t) .^ (A.1.18) Evaluating (A.1.18) at T=0 and substituting y(o), yo(o), io(o) into it, we get Z(0) =^X^ . (1+)3 Evaluating (A.1.10b) at T=0 and substituting Z,;(0) into it, we get (A.1.19) 99 X u, 311(°)(0) 0 +^= (0) Zi (0)^(1+ pt) (A.1.20) Evaluating (A.1.10a) at T=0, we get ^(o) ^ y;(o) =^+ ii(o) +^+ (p. -^(o)^ 1+ (A.1.21) If 37;(0) is known, then (A.1.20) and (A.1.21) form a set of simultaneous algebraic equations for the unknowns But y,(o) and MO). --37(0) can be calculated from the inner solution using the matching slope technique. Since y(o) is a first-order term, we need to determine the first-order inner expansion. When the first-order perturbation expansion of the inner solution is expressed in terms of the outer independent variable c, there is a contribution to it from the second-order solution y2(a). An analytic solution of y2(cy) can be obtained from Izt + ?Of zi^ct^Yi^zo (Oct + f Yo^zi (Oct•^(A-1 -22) 0^0^0^0 Y2 (a) -1311(0 — The complete solution is tedious; however, we just need the first-order contribution when 82 y2(a) is written in terms of t, namely 62 Y2(a) -= 6° term + 6 'I 2 A' g" 0(62). +^+ The first-order expansion of y(T/c) is Y(T/g) — Yo^+ 6 y,(T/g)+ c2 y2(T/c)+--6° term + [j^— X^2X la (1 +^X) ± ow] ± 0(62). ^( 1+t) (1 + (A.1.23) Hence, the first-order slope is the derivative of the coefficient of c with respect to t, i.e., 2X 4(1+ — 20/0 +^. The matching condition implies 100 1_1(1+11-2) F(0) =^ (1+ 04^• (A.1.24) Substituting (A.1.24) into (A.1.21), and solving (A.1.20) and (A.1.21) simultaneously, we obtain yi(o) = ( 1+11+ x) (0) p (-1 + + 2 X) +11)4 The yi(o) condition we obtain here using the matching slope technique agrees with the one obtained by Heineken et al. in 1967 [2]. 101 Appendix 2 General Theory of Propagating Forward and Backward Equations Consider the following equation dx f( dt^`x). (A.2.1) If we apply the mean-value theorem to x(t+h), we get x(t + = x(t)+ h^p)t + p(t + h))^ (A.2.2) for In order to apply (A.2.2) to (A.2.1) we prove a lemma Lemma 1 For a monotone increasing or decreasing function, x(t), which can be expanded in power series, then x(t + p h) = (1 — p) x(t) + p x(t + h) + 0(h2) ^ (A.2.3) Proof Since x(t) can be expanded in terms of power series, it suffices to prove the lemma for x(t)t', namely (t + p^=^+ p (t + h)n + 0(h2)^ (A.2.4) where p is independent of n. This task can be achieved by expanding both sides of (A.2.4) by using the binomial theorem and equating first-order terms in h. (QED) Using (A.2.1) in (A.2.2) gives x(t + h) x(t) + h f(x((1 — p)t + p(t + h))) ^ = x(t) + h f(x(t + p h)). (A.2.5) Apply Lemma 1 to (A.2.5), which simplifies to x(t +^x(t) + h f((1 — p) x(t) + p x(t + h))^ (A.2.6) 102 If we rename x(t+h) as xk_,1, x(t) as xk, and h as At then (A.2.6) can be rewritten as a difference equation, i.e., Xk^h 1 —^+ p^)• (A.2.7) Now we consider two cases I. p=0 The difference equation is xk+i^xk h xk ). (A.2.8) which we call the forward propagating equation. II. p=1 The difference equation is x, = x, + h xk X k+1 x„), or h (A.2.9) which we call the backward propagating equation. Next we are going to examine whether (A.2.9) is an inverse of (A.2.8). In other to do this, substitute (A.2.9) into (A.2.8) and get 0 = At f(x,)— At f( xk + At f(xk )).^ (A.2.10) If the Taylor series expansion is applied to (A.2.10), the right hand side of (A.2.10) becomes 0(At2). Thus, (A.2.9) is an inverse of (A.2.8) with an error of 0(At2). We illustrate this by an example of f(x)=-x. Example When f(x)=-x, the forward propagating equation becomes x k+1 =^At x k^ (A.2.11) the backward propagating equation becomes, xk x, + At xk „ . ^ (A.2.12) 103 Solve the forward equation numerically for At=0.01, x=1 at t=0 and then solve the backward equation numerically for At=0.01 and initial condition obtained from the solution of the forward propagation equation. The result is plotted in Figure A.1 0.1 0 . 2^0 . 3 0 . 4^0.5 0.9 0.8 0.7 0.6 Fig. Al Numerical solutions of the forward and backward equations using (A.2.11) and (A.2.12). We can see a deviation at the beginning of the graph which represents an error of 0(0.012). Since, in our example, f(x) is a linear function, the forward propagation can be inverted to obtain the backward propagating equation, which is Xk = xk „ 1- At (A.2.13) Similarly, if we solve (A.2.11) and (A.2.13) numerically with the same initial condition, Figure A.2 is obtained. 104 Fig. A2 Numerical solutions of the forward and backward equations using (A.2.11) and (A.2.13). This time the deviation at the beginning of the graph disappears. 105 Bibliography [1] A. Fersht, "Enzyme Structure and Mechanism", W.H. Freeman And Company, New York, 1985. [2] F. Heineken, H. Tsuchiya, and R, Aris, "On the mathematical status of the pseudo-steady state hypothesis of biochemical kinetics", Math. Biosci. 1, 95113, 1967. [3] L.A. Segel, "On the validity of steady state assumption of enzyme kinetics", Bulletin of Mathematical Biology 50, 579-593, 1988. [4] L. Edelstein-Keshet, "Mathematical Model in Biology", Random House, New York, 1987. [5] A. H. Nayfeh, " Introduction to Perturbation Techniques", Wiley, New York, 1981. [6] N. Wiener, "Extrapolation, Interpolation & Smoothing of Stationary Time Series", The M.I.T. Press, Cambridge, Mass., 1949. A.N. Kolmogorov, "Interpolation and Extrapolation", Bull de l'academie des sciences de U.S.S.R., Ser Math. 5, 3-14, 1941. [7] R.E. Kalman, "A new approach to linear filtering and prediction [8] problems", J. Basic Eng., Trans. ASME. Series D 82, 35-45, 1960. [9] T.P. McGraty "Stochastic Systems and State Estimation", J.Wiley & Sons, New York, 1974. [10] A. Papoulis, "Probability, Random Variables, and Stochastic Processes", Third Edition, Mc-Graw Hill, New York, 1991. [11] G.R. Grimmett and D. R. Stirzaker, "Probability and Random Processes", Oxford, University Press, Oxford, England, 1990. [12] A.H. Jazwinski, "Stochastic Processes and Filtering Theory", Academic Press, San Diego, 1970. [13]^B. Anderson and J. Moore, "Optimal Filtering", Prentice-Hall, Englewood Cliffs, NJ, 1979. 106 [14] A. Bensoussan, "Stochastic Control of Partially Observable Systems", Cambridge University Press, Cambridge, 1992. [15] H.W. Sorenson and D.L. Alspach, "Recursive bayesian estimation using gaussian sums", Automatica, AC4, 465-479, 1971. [16] D.L. Alspach and H.W. Soreson, "Nonlinear bayesian estimation using gaussian sum approximations", IEEE Trans. Automatic Control, AC-17, 439448, 1972. [17] G. Caminal, J. Lafuente, J. LOpez-Santin, M. Poch, and C. Solk "Application of extended kalman filter to identification of enzymatic deactivation.", Biotechnol. Bioeng., 29, 366-369,1987. [18] G. Caminal, J. LOpez-Santin, and C. Solk "Kinetic modeling of the enzymatic hydrolysis of pretreated cellulose.", Biotechnol. Bioeng., 27,12821290,1985. [19] G. Greco, G. Marrucci, N, Grizzuti, and L. Gianfreda, Ann. NY Acad. Sci. 434, 7, 1984. [20] A. Jazwinski, "Brief paper". Automatica 10, 203-207, 1974. [21] F. H. Schlee, C.J. Standish, and N.F. Toda, Divergence in the Kalman Filter, AIAA J. 5, 1114-1120, 1967. [22] R.J. Fitzgerald, Divergence of the Kalman filter, IEEE Trans. Automatic Control, AC16, 736-747, 1971. [23] A. Jazwinzki, "Limited memory optimal filtering", IEEE Trans Aut. Control 13, 588, 1968. [24] H. Cramer, "Mathematical Method of Statistics." Princeton Univ. Press Princeton, New Jersey, 1961. [25] H.W. Sorenson and A.R. Stubberud, "Non-linear filtering by approximation of a a posteriori density", Inter J. Control 8, 33-51, 1968. [26] A.S. Lapedes and R.M. Farber, "Non-linear signal processing using neural networks : prediction and system modeling", Los Alamos preprint LAUR-87-2662, 1987.
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Application of the extended Kalman filter to enzyme...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Application of the extended Kalman filter to enzyme reactions Lai, Henry K. H. 1993
pdf
Page Metadata
Item Metadata
Title | Application of the extended Kalman filter to enzyme reactions |
Creator |
Lai, Henry K. H. |
Date Issued | 1993 |
Description | In monitoring biological processes, measurement of key variables is often impeded by the lack of reliable on-line measurements of biomass, substrate, and product concentrations, and the difficulty to properly model biological activity in processes that are nonlinear and time-varying. One approach to solving this problem involves the development of state estimation techniques that use available measurements to reconstruct the evolution of state variables or to estimate the bioprocess parameters. The use of filtering theory for state estimation provides a means of incorporating a deterministic model into a method of forecasting future states which includes the probabilistic uncertainties of both the system and the measuring devices. The state estimation technique we use is the discrete extended Kalman filter. This method allows the use of an approximate model and partial measurements of process variables. To study the application of the discrete extended Kalman filter to a biological problem, we used two different enzyme reactions as model systems. The first model problem is a simulation of cellulose hydrolysis with enzyme inactivation. The second model problem is a simulation of the separation of a chiral substance into its respective enantiomers. In the first problem, the hydrolysis of cellulose, a slowly varying parameter is tracked after correctly choosing the model error covariance matrix. In the second problem, we reconstruct the complete composition of the system from a partial set of measurements. Although these model problems are simple cases involving only single enzymes, the extended Kalman filter can be applied to more complex systems of enzymes. |
Extent | 3893416 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2008-09-16 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0079750 |
URI | http://hdl.handle.net/2429/2036 |
Degree |
Master of Science - MSc |
Program |
Mathematics |
Affiliation |
Science, Faculty of Mathematics, Department of |
Degree Grantor | University of British Columbia |
Graduation Date | 1993-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1993_fall_lai_henry.pdf [ 3.71MB ]
- Metadata
- JSON: 831-1.0079750.json
- JSON-LD: 831-1.0079750-ld.json
- RDF/XML (Pretty): 831-1.0079750-rdf.xml
- RDF/JSON: 831-1.0079750-rdf.json
- Turtle: 831-1.0079750-turtle.txt
- N-Triples: 831-1.0079750-rdf-ntriples.txt
- Original Record: 831-1.0079750-source.json
- Full Text
- 831-1.0079750-fulltext.txt
- Citation
- 831-1.0079750.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0079750/manifest