- Email: [email protected]
[25] Simulation of hemoglobin kinetics using finite element numerical methods
[25]
FINITE ELEMENT NUMERICAL METHODS
517
[25] S i m u l a t i o n o f H e m o g l o b i n K i n e t i c s U s i n g F i n i t e Element Numerical Methods
By
ROBERT L. BERGER, NORMAN DAVIDS,
and
M I C H E L E PERRELLA
Introduction Ordinary differential equations (ODEs) have been used for centuries as mathematical models of physical and chemical processes. The laws of mass action, for example, can be used to derive an ODE or system of ODEs describing a kinetic process. The advantage of using ODEs to express complex models arising from the laws of mass action is that they are easy to derive, even if they are difficult or impossible to solve analytically. Numerical solutions to such equations have been one of the great advantages of the computer age. The numerical analysis literature is full of solution methods, and a survey of that literature is beyond the scope of this chapter. The interested reader should see such topics as ordinary differential equations, stiff differential equations, initial value problems, and boundary value problems. In this chapter, we discuss one elementary method of analyzing biochemical kinetic problems to show what some of the issues are, including a discussion of how such methods appear to a user, in terms of ease of use. The general solution of a complex reaction, in this case hemoglobin (Hb) binding to carbon monoxide (CO), is then discussed in detail, including a global fit using a curve fitter. Finite Element Approach In simulating a rate process, we start from the law of mass action, which states that the rate of the forward reaction is proportional to the concentrations of the reactants, say A and B. If these are written as [A] and [B], Of = kf[A][B] and in a similar manner, vb = kb[C][D] if C and D are the products) In the simplest case we might have reactant A going to B, with the back reaction going so slowly that it can be neglected. Then we have
dA/dt = -kf[A]
(1)
where kf is the rate constant for the reaction and [A] is the concentration of the particular component. I W. J. Moore, "Physical Chemistry," 3rd ed., p. 169. Prentice-Hall, Englewood Cliffs, NJ, 1962.
METHODS IN ENZYMOLOGY,VOL. 232
Copyright © 1994by AcademicPress, inc. All rights of reproduction in any formreserved.
518
MATHEMATICAL ANALYSIS AND MODELING
[25]
Equation (1) may be interpreted in two ways: (1) that it is an ODE, and (2) that it is a statement about finite steps dA and dt. In the first case the ODE must be solved, which in this elementary case is easy, the solution in closed form being A = A0 e x p ( - k f t ) . We can, however, ignore this solution completely by refraining from having dA, dt pass to the limit, which is zero. In this way we bypass having to solve the differential equation either in closed form or numerically, because it is not generated in the first place. Thus Eq. (1) can be rewritten as dA = -kf[A] dt and used repetitively to provide an accounting of the increments of reaction products being created during a succession of time intervals dt. Following this step, we adjust or update the value of A by the relation [A]' = [A] + d[A]
(2)
These operations develop into a computer program when these steps are carried out for each reaction and each set of reactants, and include the p r o p e r i n p u t - o u t p u t statements and control statements. The particular scheme for doing this was developed by Davids and Berger 2,3 and is called the finite element or F E method. The above steps are inherent in the use of any digital procedure for simulating continuous processes: (I) The time must be divided into a discrete set of intervals so that the state of the system, stored as numerical values, can be updated repetitively, as described above; (2) the value of the time step dt and other numerical parameters must be chosen with care so that the computations are not only accurate but, even more fundamentally, remain stable; that is, do not cause a rapidly increasing error to build up. This point is discussed in more detail in the next section.
Stability Problem-Multiple Reaction Systems An inherent problem with any procedure that describes a physical system with finite sequential steps is that instability occurs when initially small errors grow, tending to build up large fluctuations in a system. Instability problems may arise when multiple reactions exist in a system having widely varying orders of magnitude in rate. As an illustration consider the following system: A + B C + D
kl k2
) C
(3)
, A + E
(4)
2 N. Davids and R. L. Berger, Commun. ACM 7, 547 (1964). 3 R. L. Berger, Discuss. Biophys. Soc., Biophys. J. 24, 2 (1978).
[25]
FINITE ELEMENT NUMERICAL METHODS
519
TABLE I COMPARISON OF EXACT SOLUTION VERSUS FINITE ELEMENT SOLUTION a
[A]I[Ao] kit
Analytical
Simulation
[B]/[B0]
0 0.2 0.4 0.6 0.8 1.0 2.0 3.0 4.0 5.0
1.000 0.819 0.670 0.549 0.449 0.368 0.135 0.050 0.018 0.007
1.000 0.818 0.669 0.547 0.448 0.366 0.134 0.049 0.018 0.007
0.000 0.164 0.268 0.331 0.363 0.374 0.302 0.204 0.144 0.114
a Comparison of the solutions of the kinetics of Eq. (5) according to the exact mathematical solution of Eqs. (6)-(8) and the finite element method simulation. Conditions: [A]0 = 1, kl = 1, k 3 = 1, k 2 = 0.1, and dt = 0.01 sec.
in which kl is several orders of magnitude faster than k2. In such a case the stability of the numerical model is maintained by choosing dtd for the fast reaction, say dt d = (1/100) dt, and then taking 100 steps for this reaction every time a single step of the slow reaction, which uses the time step dt, is calculated. We can think of this procedure as "synchronizing the clocks" associated with each of the reactions. Several kinetic cases have been taken from the literature to validate the procedure. As an example, consider the following case of two consecutive reactions: A
kl
> B
.
k3
" C
(5)
k2 A first, irreversible unimolecular reaction is followed by a reversible one. A general treatment of such reactions as given by Benson 4 and Johnson has given a two-exponential solution5 (extensive tables for this reaction 4 S. W. B e n s o n , " T h e Foundations of Chemical Kinetics." McGraw-Hill, New York, 1960 (see esp. Chapter III). 5 K. A. Johnson, this series, Vol. 134, p. 677.
520
MATHEMATICAL ANALYSISAND MODELING
[25]
based on exact m a t h e m a t i c a l formulas are available6). The rate equations for this reaction scheme are dA/dt = -kl[A] dB/dt = k~[A] + k2[C] - k3[B] dC/dt = -k2[C] + k3[B]
(6) (7) (8)
Table I gives a c o m p a r i s o n with the result c o m p u t e d f r o m the closed f o r m solution, taken f r o m K u b a and K o n o w a l o w 6 and from the F E method, assuming [A]0 = 1, kl = 1, k 3 = 1, and k2 = 0.1. The following is the F E p r o g r a m used, written in BASIC. The reaction of bicarbonate and hydrochloric acid is a well-known e x a m p l e o f this type of reaction scheme.7 On mixing, the reactants combine quickly (108 M -I sec -t) to f o r m carbonic acid. The b r e a k d o w n to carbon dioxide and w a t e r is also fast, 25 per second at 25 ° , and the reverse reaction is slow, e v e n at p H 7.4. A p p e n d i x A contains the Fortran program for solving this case b y the finite element method, but without the curve fitter. This requires making a first estimate o f the values of the rate constants. The calculation is then m a d e and a visual appraisal is made, comparing the F E output with the experimental data. F r o m the rates calculated, a new estimate can be m a d e and the process repeated until satisfactory a g r e e m e n t is reached. Appendix B shows a m o r e complicated case with the automatic c u r v e fitter added. This point is discussed in greater detail later on page 523.
Association Reaction between Hemoglobin and Carbon Monoxide Finite element analysis can be used to analyze the interaction between hemoglobin and CO in the general case, in which both the CO association and dissociation reactions occur, along with structural transition equilibria, during the course of the reaction. We consider here a special case in which substoichiometric a m o u n t s of CO are quantitatively bound by hemoglobin in a period of time during which the dissociation of the ligand is negligible. At the end of this preset reaction time, the liganded intermediates are q u e n c h e d and separated by cryogenic techniques. 8,9 The experimental data are in the f o r m of concentrations of the intermediates as a 6 D. W. Kuba and D. D. Konowalow, "Report ARL 63-44 (ASTIA)." University of Wisconsin, Madison, 1963. 7 L. Rossi-Bernardi and R. L. Berger, J. Biol. Chem. 243, 1297 (1968). 8 M. Perrella, L. Benazzi, L. Cremonesi, S. Vesely, G. Viggiano, and R. L. Berger, J. Biochem. Biophys. Methods 7, 187 (1983). 9 M. Perrella and L. Rossi-Bernardi, this volume [21].
[25]
FINITE ELEMENT NUMERICAL METHODS
521
30
'INITIALIZE:
40
A=I
'Initial
50
B=0
'
"
"
B
60
C=0
'
"
"
C
70
KI=I
'Rate const
for A - - - > B
80
K3=I
'Rate const
for B - - - > C
90
K2=.I
'Rate const
for b a c k reac B < - - - C
i00
T=0
ii0
DT=.01
'Set t i m e
120
KM=I01
'Set n u m b e r
130
'
conc
of A
step of
intervals
140
'START
REACTION
150
FOR K=I
TO KM
160
LPRINT
170
D=KI*A*DT
'Law o f M a s s
180
A=A-D
'Update
A
190
B=B+D
'Update
B
200
D=K3*B*DT-K2*C*DT
'Law o f M a s s
210
B=B-D
'Update
of B again
220
C=C+D
'Update
of C
230
T=T+DT
'Update
of time
240
NEXT
'End l o o p
250
STOP
K
USING
"#.##
###.###
#.###
#.###";;T,A,B,C Action
for A--->B
Action
for B<--->C
SCHEME I
function of fractional CO bound. Time information is not provided because at the high values of the hemoglobin concentration required by the cryogenic technique, the CO-binding reaction is so fast that sophisticated quench-flow methods are required to study the reaction on a time-dependent basis. Program Figure la is a general schematic showing the 10 ligation states and the pathways (16 in all) that can be taken from deoxyhemoglobin (Hb) to
522
MATHEMATICAL ANALYSIS AND MODELING
[25]
24
b
r"'-'-I k3
- 1231
k7
FIG. 1. (a) General schematic showing the 10 ligation states of hemoglobin and the 16 forward pathways leading from deoxyhemoglobin to fully liganded hemoglobin. Each ligation state is assigned indexes U, w h e r e j denotes the code of the isomer in ligation state i. Each rate constant k is assigned subscripts ijk: i, matrix from which transition originates;j, matrix element where transition terminates; k, number of sites liganded. (b) Simplified 12-rate constant path. Species 21 and 22 are grouped together because they cannot be resolved by the cryogenic technique for the isolation of the intermediates.
[25]
FINITE ELEMENT NUMERICAL METHODS
523
carbonmonoxyhemoglobin (HbCO). Each ligation state is located by the number of liganded sites (horizontal position). In a previous notation each state was also located by its position from the top of the diagram (columnwide). 3 We have changed this notation for the double-liganded species to conform to another such scheme now commonly used. l° A variable defining the rate constants is now introduced as (KON)iik, where i identifies the matrix element from which transition originates; j identifies the matrix element where transition terminates; and k is the number of sites liganded. Care must be taken in assigning statistical factors, 1 or 2, for the type of heme available for any given path. A value of zero for the variable corresponding to a particular transition in Fig. la prevents that transition from taking place. By putting in enough zeros, one can narrow the process to a single pathway or group of pathways. The general scheme for doing the simulation and optimizing the rate constants by global curve fitting is presented below, with applications to a simple fourstep model of CO hemoglobin binding and then using a 12-rate constant model (Fig. lb). Outline of Kinetics Curve-Fit Program Standard methods exist for determining the best fit by least squares to a set of data points, using a particular type of function, such as a polynomial or a series of exponentials, l~ The method of finite elements starts with the differential equations for the system in finite form. One then supplies the kinetics program with a set of trial parameters, that is, starting values or values determined from a previous iteration. We then go through the kinetics program and determine the error between the predicted curve and the data. These errors are then used to calculate the next mean-square error, which is compared with the previous one. A decision is then made as to whether to repeat the process with a disturbed set of parameters, until an optimum is reached. In this manner, it is possible to curve-fit a set of data by using kinetic models in place of numerical functions altogether. Figure 2 shows the scheme diagrammatically. 1. Read in and store the data, as identified by file name. 2. Enter the initial concentrations of the solutions, [Hb] and [CO]. 3. Enter the assumed transition sequence case. 4. Enter the assumed starting values for the rate constants. 5. Enter the number of iterations desired. l0 G. K. Ackers and F. R. Smith, Annu. Rev. Biophys. Biophys. Chem. 16, 583 (1987). 11 D. W. Marquardt, S I A M ] . 2, 431 (1963).
524
[25]
MATHEMATICAL ANALYSIS AND MODELING INITIAL VALUES of Rate Constants Bol, etc. Hb, CO for
~ 1
[
CURVE - FIT
I
[ [
~
Are these the best parameters?
Nfo-gl
I KINETICSl
FIG. 2. Outline of the finite element solution for the h e m o g l o b i n - C O reaction.
6. The c o m p u t e r then begins by passing the values of the given rate constants to the subroutine CFIT, where they are stored as the vector: parm;, i = 1, nparm. 7. These values are then passed to subroutine K I N E T , which works out a prediction function y(t) according to the transition sequence described there. The error is then calculated from e(t) = ydat(t) - y(t), where ydat(t) equals the observed data. 8. We then return to C F I T to obtain the sum of squares: s = X[e(t)] 2. 9. Each parameter is then perturbed in turn by the relation parm; = parm;(1 + 0.001). 10. With the above steps, the gradient vector can be set up normal to the surface o f a constant sum of squares in parameter space. The Jacobean matrix of the system can then be constructed. The inversion and diagonalization of this matrix, according to the G a u s s - N e w t o n procedure, provides an updated set of parameters and a new sum o f squares o f error, hopefully lower in magnitude. 11. The cycle of iterations terminates when the mean-square error converges. This can be determined by specifying a preset tolerance or the n u m b e r of iterations to be made. The user should review the display o f the successive sum of squares and then, by typing an integer, select additional iterations as desired.
Global Fit In general, any curve-fit scheme will yield more precise values for the parameters when more data points are available. In fact, for any kind of precision we must have the condition ndata >> nparams
(9)
[25]
FINITEELEMENTNUMERICALMETHODS
525
In the simple four-step scheme of the reaction, in which the subunit heterogeneity is ignored, there are five values of the concentration of the intermediates, B(i), available for determining four parameters, if time information is provided. If time information is not provided, there are three parameters relative to a fourth fixed parameter for making the determination. Because the sum of the concentrations of the intermediates at a certain value of the fractional CO bound equals the total hemoglobin concentration, only four values of the concentrations of the intermediates are independent. Thus for the simple four-step scheme it is possible to curve-fit the parameters by using the data from a single experiment at some fractional value of CO bound to hemoglobin. If data from experiments carried out at many different values of the fractional CO are available, a more reliable estimate of the error in the curve-fitted parameters is obtained by combining all of the data into what is called a "global" fit. This is done by simply "stacking the data" until we run out of data, as shown in Table II. Thus, if we have n = 10 experiments to work with in this fashion, we obtain 90 independent data points for the adjustable parameters relative to the reference parameter and condition (9) is fulfilled. If subunit functional heterogeneity is considered, I0 different ligation states are possible according to Fig. la. Because the cryogenic technique
TABLEII STACKINGMETHOD [BH] [B21] [B°~]t experiment1 [B31] [B41] [BI2] t [Bo2] [B22] experiment2
[B32] [B42]
[B.] I experimentn
526
MATHEMATICAL ANALYSIS AND MODELING
[25]
cannot determine the concentrations of intermediates B2~ and B22 separately, 9 values of the concentrations of the intermediates, B(i), are determined, of which only 8 values are independent, against a total of 12 parameters to curve-fit for, as shown in Fig. lb. Thus, even if k~l and k12 are assumed to be equal reference parameters, the number of adjustable parameters (10 values of the k) is greater than the number of data points, that is, 8 concentration values, and condition (9) is not fulfilled. In such a case the global fit using the data from several experiments at different values of fractional CO bound is mandatory. The Fortran program used to do this is shown in Appendix B.
Application to Simple Four-Step Scheme Because the functional heterogeneity of the subunits in the association reaction with CO is slight, 12the general process in Fig. la can be simplified by reducing it to a simple four-step process by putting kll I = kl21 k132 = k242
(10) (11)
kll 2 = k212 = k222 = k122
(12)
k313 = k423
(13) (14)
kll 3 -- k123 = k213 = k223 kll 4 = k214
(15)
Furthermore, when no time information is provided by the data, either the rate constants are calculated relative to a reference rate in the process or such a reference rate is assigned a value taken from the literature. Table III shows the results of the analyses of the experimental data, obtained by quench-flow and cryogenic isoelectric focusing methods s,9 to give the intermediates, to obtain the rate constants in the association reaction 12 between Hb and CO at 20°, in 0.1 M KC1, pH 7, assuming k4 = 6.0/zM -1 sec -l, as obtained by Gibson using optical stopped-flow techniques ~3 under similar conditions. The optimum time step dt was determined to be 0.01 msec. The analyses, starting with the data at the lowest fractional value of CO bound, were carried out by stepwise addition of the data at higher CO saturations to obtain the final global fit to all the
12 M. Perrella, N. Davids, and L. Rossi-Bernardi, J. Biol. Chem. 267, 8744 (1992). 13 Q. H. Gibson, J. Biol. Chem. 248, 1281 (1973).
[25]
FINITE ELEMENT NUMERICAL METHODS
527
TABLE III RATE CONSTANTS FOR REACTION OF HEMOGLOBIN PLUS
CO: FIVE-SPECIES MODELa
Fitted parameters (p.M -I sec -l)
Global fit up to experiment no.
Sco (%)
kl
k2
k3
SSE b (× 104)
I 2 3 4 5 6 7 8 9 10c
17.6 18.5 21.8 30.1 39.9 45.3 54.9 59.8 72.4 82.7
0.0777 0.0952 0.0944 0.0975 0.103 0.108 0.110 0.112 0.113 0.113
0.199 0.227 0.226 0.231 0.241 0.253 0.256 0.261 0.264 0.265
1.081 1.302 1.333 1.336 1.412 1.472 1.474 1.465 1.474 1.472
11.403 6.231 6.618 5.974 5.638 5.445 5.436 5.376 5.371 5.368
0.1024
0.242
1.382
6.286
Average valuesd:
a Rate constants for the association reaction between Hb and CO at 20°, pH 7, from data on the reaction intermediates as studied by cryogenic techniques. ~2Parameters (k 1, k2, k3) (intrinsic rate constants, i.e., apparent rate constants divided by the statistical factors 4, 3, and 2, respectively) were fitted by the finite element method as a function of the cumulative number of experiments at different CO saturations Sco, (1, 1 + 2 . . . . . 1 + 2 + -.- + 10), using k4 = 6.0/xM -1 sec-l. 12 b SSE, Sum of square errors (see text). c The error in the rate constants is - 10%. d Rate constants obtained by averaging the values of the parameters fitted to the data at each CO saturation value.
d a t a . T h e c a l c u l a t e d r a t e c o n s t a n t s at e a c h s t e p w e r e u s e d to fit t h e d a t a at all C O s a t u r a t i o n s ; t h e s u m o f s q u a r e e r r o r s f o r e a c h fit ( S S E ) a r e s h o w n in T a b l e I I I . Clearly, limiting values for the rate constants are reached inthe stepw i s e p r o c e s s a n d t h e m i n i m u m v a l u e f o r S S E is a t t a i n e d w i t h t h e p a r a m e t e r s f r o m t h e g l o b a l fit. F i g u r e 3 s h o w s t h e r e l a t i v e c o n c e n t r a t i o n s o f t h e s p e c i e s in d i f f e r e n t s t a t e s o f l i g a t i o n v e r s u s C O s a t u r a t i o n as c a l c u l a t e d b y t h e finite e l e m e n t method, assuming negligible functional heterogeneity of the subunits and u s i n g t h e v a l u e s o f t h e r a t e c o n s t a n t s o b t a i n e d b y t h e g l o b a l fit. If we now expand the possible pathways by utilizing the scheme of F i g . l b , f i x i n g t h e l a s t 2 r a t e c o n s t a n t s , k u = kl2 = 3.0 ~ M -1 s e c -1, a n d a g a i n d o t h e g l o b a l fit u s i n g all 10 e x p e r i m e n t s , w e o b t a i n t h e r e s u l t s s h o w n in T a b l e I V .
528
M A T H E M A T I C ANALYSIS AL AND MODELING
[25]
0.8
0
4
0.6
fi 0.4
0.2
0
0.2
0.4
0.6
0.8
Sco
FIc. 3. Fraction f/of liganded intermediate i (i = 0-4) versus CO saturation Sco in the association reaction between Hb and CO at 20° in 0.1 M KCI, pH 7) 2 The curves were calculated using the rate constants obtained by the global fit to the experimental data on the concentrations of the intermediates (see Table III), assuming negligiblefunctional heterogeneity of the subunits. T h e results o f this global fit to the CO-binding data, according to the 12-rate constants s c h e m e in Fig. Ib, confirm the finding of Perrella et a l ) 2 that the fl subunit reacts faster than the a subunit in the first CO-binding step (k2/k~ = 1.91 - 0.28). A similar subunit heterogeneity is o b s e r v e d b y comparing the rate constants for the second CO-binding step along the a a - ~ f l p a t h w a y (upper p a t h w a y s in Fig. Ib) and along the flfl-aa p a t h w a y (lower p a t h w a y s in Fig. lb), which yield the ratio k6/k 3 -- 2.98 - 1.76. Because species 21 and 22, both having one ct- and one fl-liganded subunit, are not equivalent and the cryogenic technique allows only the determination o f the sum of their concentrations, the rate constants k4, ks, k8, and k 9 are c o m p o s i t e constants. H o w e v e r , within experimental error, k 4 = k 5 and k 8 = k9. Thus the global fit confirms a previous conclusion that the CO-binding properties o f species 21 and 22 are similar ~2and that all binding steps are c o o p e r a t i v e with an acceleration after the binding of the second CO molecule. The simple model s h o w n above, and the m o r e complex 12-constant model, can be solved by analytical methods and the results are in close agreement. H o w e v e r , w h e n one goes to the reactions with o x y g e n and hemoglobin, analytical models of the 12 rate constants, including the b a c k reactions, which cannot be neglected as in the CO-binding case, b e c o m e
[25l
FINITE ELEMENT NUMERICAL METHODS
529
T A B L E IV RATE CONSTANTS FOR REACTION OF CO PLUS HEMOGLOBIN: NINE-SPECIES MODELa Ligation of ct subunits
Ligation of/3 subunits
kI k3 k5 k8 kj0
k2 = k0 = k6 = k7 = k9 =
= = = = =
0.0410 -+ 0.0042 0.0438 - 0.0215 0.153 -+ 0.054 0.454-+ 0.117 0.729 -+ 0.363
0.0782 -+ 0.0080 0.162 --- 0.03I 0.131 -+ 0.0428 0.133 +- 0.220 0.713 -+ 0.159
a Rate constants, k. (n = 1-12) b, for the binding of CO to deoxyhemoglobin at 20 ° in 0.1 M KCI, pH 7, calculated from the concentrations of the intermediates according to the reaction scheme of Fig. lb, assuming kll = k12 = 3.0 tzM -~ sec -~, that is, one-half the value measured by optical stopped-flow kinetics for the binding of the last CO moleculeJ 3 b Intrinsic, which is apparent constant k', divided by the statistical factor (in ~ M -j sec-1). The concentration is the heine concentration, kj = k l / 2 ; k 2 = kr/2; k 4 = k ' J 2 ; k 5 = k'5/2; ks = k~/2; k 9 = k ¥ 2 ; in all other cases k, = k,~.
intractable in general. We thus believe that the finite element method offers an important technique for solving or simulating large kinet~ problems. General
Conclusions
The FE method has the advantage that its steps are form,u'lated easily, as they deal directly with the physics or chemistry o f the problem. It simulates the problem in the sense that if the computation is stopped for any reason at an earlier time, it shows the state of the system at that time. Furthermore, it can easily deal with an additional boundary condition or the presence of additional effects, such as electric or magnetic fields, without invalidating the work up to that point. On the other hand, instabilities are difficult to deal with, just as with other numerical methods. In a more general sense, solving kinetic equations requires a number of methods. We have found that the program MATLAB 14 serves very well. We have presented a method of solving kinetic equations using any computer and a language such as Fortran, but we have found that the 14 MATLAB, PC or 386 version available from Mathworks, Inc., 24 Prime Park Way, Natick, MA 01760.
530
MATHEMATICAL ANALYSIS AND MODELING
[25]
utilization of MATLAB allows any of the numerical methods to be written down without concern of portability, as these programs will run on any version of MATLAB. These programs are still under development. It should be pointed out, however, that approaching a problem by numerical analysis of ODEs gives mathematical analysts the advantage of having at their disposal many established numerical and error analysis methods. They also have at their disposal such computer procedures as MATLAB, and so on, that contain within them the broad experience of many other analysts and that are universal in scope as well as convenience. Initial efforts to put the FE method into MATLAB and to utilize stiff differential equation solvers as well as curve fitters seem to be working and it is hoped these efforts will be reported in detail in the near future. 15 Appendix A: Program without Curve Fitter C
PROG:
C02.F4,
CARBON
DIOXIDE
KINETIC
MODELS
C C
REACTIONS
IN ACID
SOLUTIONS
C
COMPONENTS:
C
C02
= CARBON
C
H2CO3
= CARBONIC
C
HCO3-
= BICARBONATE
C
C03--
= CARBONATE
C
TOT
= TOTAL
C
PC02
= GASEOUS
C
KC02
= HYDRATION
C
KH2CO3=DEHYDRATION
C
HEATS:
DIOXIDE
(DISSOLVED)
ACID ION(-) ION
CARB.
(--),
DIOX.
HCO3=H+C03
CONTENT
F O R M O F C02, VELOCITY
OF SYSTEM
MEASURED
CONSTANT
VELOCITY
AS PARTIAL
PRES.
(THERMAL)
CONST.
C
DELHI
= HEAT
OF DE-IONIZATION
CAL/MOL
C
DELH2
= HEAT
OF DEHYDRATION
CAL/MOL
C
PROG.
CONTROL
C
KS=0,
TERMINATE
C
KS=I,
RECYCLE
DIMENSION
(=CO2+H2CO3+HC03-)
SWITCHES EXECUTION
PROGRAM
Y(500),TAU(500),Q(500)
15 R. I. Shrager, Division of Computer Research and Technology, National Institutes of
Health, Bethesda, MD (personal communication).
[25]
FINITE ELEMENT NUMERICAL METHODS
REAL*4 N A , N A H C O 3 , K C O 2 , K H 2 C O 3 , K C L , K H , K H 2 0 , K H C O 3 , K L 2 C 0 3 REAL*4 I S T R , M U C O E F F IST D I S S O C H2CO3
(KI)
A=3404.71 D=14.8435 C=0.032786 CONC OF H20
(MOLS/L)
H20=55.5 P R I N T 804 804
FORMAT
(' T E M P E R A T U R E O F E X P E R
(DEG-C):
R E A D 3, T E M P TABS=TEMP+273.16 H2CO3 DISSOC CONST PKI=A/TABS-D+C*TABS EI=I0**(-PKI) WRITE(21,911)TEMP,PKI,EI IF
(TEMP.EQ.37.) GO TO i0
ASSUMED RATE CONSTANTS
(/SEC):
D I S S O C D A T A F R O M LONG, BIOCHEM.HANDB. P43 D E H Y D R A T I O N C A R B O N I C ACID: H 2 C O 3 = C O 2 + H 2 0 VELOCITY CONST
(AT 25 DEG)
KH2C03=25.5 BACK REACTION RATE
(HYDRATION)
KCO2=.0257/H20 2ND D I S S O C C O N S T HCO3 PK=I0.329 A=2902.39 D=6.4980 C=0.02379 GO T O 20 C i0
CONTINUE
(K2), 25 D E G
'$)
531
532
MATHEMATICAL ANALYSIS AND MODELING VELOCITY
CONST
AT
37 D E G
KH2C03=72.00 KCO2=.0620/H20 20 C C
CONTINUE DISSOC
CONSTS
IONIZATION
C
EQUILIBR
C
PKWHEN
VS TEMP & FORMATION
CONST
H2CO3
CORRECTED
AT
FOR
OF CARBONIC 25 D E G
SALT
ACID:
H+HCO3=H2C03
IS 3 . 7 6 5
EFFECT
PK=3.5 PK AT
37 DEG:
PK=3.6 EH2CO3=I0**(-PK) ASSUMED
FORMATION
RATE
OF H2CO3
FROM
H+HCO3-
KHCO3=SEI0 DISSOCIATION
RATE
(IONIZATION)
OF H2CO3
KL2CO3=EH2CO3*KHCO3
FORMATION
OF H20:H+OH=H20
KH=I. 4Eli KH20=2.5E-5 DELHI=I205. DELH2=I455. SOLUBILITY
CONSTANTS
A=0. 0301 COMPUTER NFIL=23 KM=I00 K2=0 KPR = 1 VOL=I. 0 T=Oo 0 TI=0.0
INITIALIZATION
& CONTROL
CONSTS
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS TW=0.0 K=0 TW=0 QI=0. Q2=0. QMAX=0.
C
P R E P A R A T I O N OF S O L U T I O N
i: HCL
CALL A N M O D E WRITE(NFIL,912) 912
FORMAT
(' S O L U T I O N
1 (HCL)')
PRINT 801 801
FORMAT
(' M O L A R CONC SOLUTION 1 (HCL):
)
R E A D 3,HCL 3
FORMAT
(2F)
PRINT 802 802
FORMAT
(' CONC SALT ADDED:
)
R E A D 3,KCL C
H C L COMES TO E Q U I L I B R I U M
C
E Q U I L I B R I U M CONST H C L
(IONIZES):
HCL=H+CL
HI=HCL CLI=HCL PH0=-ALOGI0(H1) H=H1
PH=PH0 C
CHARGE B A L A N C E CHG=H-CL WRITE(NFIL,901)
C
PREP.
SOLUTION
K,K2,T,PH,HUM,OHUM,HCO3M,H2CO3M,CO2M,CHG
2:NAHCO3
PRINT 803 803
FORMAT
(' M O L A R CONC S O L U T I O N 2 ( N A H C O 3 ) : )
READ 3,NAHCO3 C
NAHCO3
COMES TO EQUILIBRIUM:
NAHCO3=NA+HCO3
533
534
MATHEMATICAL ANALYSIS AND MODELING NA=NAHCO3 HCO3=NAHC03 HCO3M=HCO3*IE3 WRITE(NFIL, 9 0 1 ) K , K 2 , T , P H , H U M
C 1 C
R E S T A R T OF C O M P U T A T I O N
HERE
CONTINUE M I X I N G A F T E R 1 MILLISEC,
NAHCO3+HCL=H2CO3+NACL
WRITE(NFIL, 903) 903
FORMAT
(IH ,'SOLUTIONS
1 & 2 B E C O M E MIXED')
CL=CL1/2.0 NA=HCO3/2.0 HCO3=HC03/2.0 H=HI/2.0 H2CO3=0.0 H2CO3M=H2CO3*IE3 HUM=H*I.E6*VOL OHUM=OH*I.E6*VOL C02M=CO2*IE3 CALL MOVABS(IX0,IY0) CO2=0.0 C
C A R B O N I C A C I D FORMED:
N+HCO3=H2C03
DT=5E-10 DO i01 K 2 = 1 , 4 0 D=KHCO3*H*HCO3*DT-KL2CO3*H2CO3*DT H2CO3=H2CO3+D H=H-D HCO3=HCO3-D C
EXIT AT E Q U I L I B R I U M TOL=0.1E-5 IF
(ABS(D).LT.TOL)
T=T+DT i01
CONTINUE
GO TO 13
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
13
CONTINUE
C
CHG BALANCE CHG=H+NA-CL-HCO3 HI=H
C
CONSTANTS
IN R O S S I - B E R N A R D I - B E R G E R
PAPER
A=.07 B=.I
CI=2.303*(I.+(EH2CO3/(B-A))) C2=2.303*(I.-(EH2CO3/(B-A))) C
IONIC STRENGTH AT START OF REAC
(M/L)
ISTR=0.1 C
DEBYE HUCKEL FORMULA MU=ISTR DB=0.505*SQRT(MU)/(I.+I.6*SQRT(MU))
C
MEAN IONIC ACTIVITY
COEFF
F=I0**(-DB) CI=2.303*(I.+(EH2CO3/F/(B-A))) C2=2.303*(I.-(EH2CO3/F/(B-A))) WRITE(21,911)CHG,DB,F,EH2CO3 911 C
FORMAT
(IH ,4E13.3)
KINETICS
LOOP HERE
DTI=.001 TMAX=KM*DTI DO 102 K=I,KM PH=-ALOGI0(N) CHG=H-HCO3-OH-CL+NA HCO3M=HCO3*IE3 N2CO3M=H2CO3*IE3 CO2M=CO2*IE3 IF((K-I)/EPR*KPR-(E-I))II,12,11 12
CONTINUE
535
536
MATHEMATICAL ANALYSIS AND MODELING P A R T I A L PRES.
CO2
(MMHG)
PCO2=(CO2M+H2CO3M) IA HUM=H*I.E6*VOL OHUM=OH*I.E6*VOL WRITE
(NFIL, 901) K , K 2 , T I , P H , H U M , O H U M , H C O 3 M , H 2 C O 3 M , C O 2 M , C H G
901
FORMAT(IH
ii
CONTINUE
,213,E10.2,F6.2,2X,2EI0.3,2X,3E10.3,2X,2EI0.3)
C
KINETICS
PLOT
C
PLOT OF PH VS T I M E IX=T1/TMAX*IXINT+IXO IY=(PH-2)/5.*IYINT+IY0 CALL DRWABS(IX,IY)
D E H Y D R A T I O N OF C A R B O N I C ACID:
H2CO3=H20+CO2
D=KH2CO3*H2CO3*DTI-KCO2*H20*CO2*DTI H2CO3=H2CO3-D C02=CO2+D H20=H20+D I N C R E M E N T A L HEAT P R O D U C E D BY D E H Y D R A T I O N DQ2=D*DELH2 Q2=Q2+DQ2 C
I O N I Z A T I O N OF W A T E R
C
H+OH=H20
C
CALL WATER(H,OH,H20,KH20,KH,PH)
C C 15
I O N I Z A T I O N OF C A R B O N I C A C I D
(H+)+(HCO3-)=H2CO3
DO 103 K 2 = I , 4 0 HUM=H*I.E6*VOL
C 902
WRITE
(NFIL,902)
FORMAT(IH
K2,H,HCO3,H2CO3,D
,I3,4E10.3)
D=KHCO3*H*HCO3*DT~KL2CO3*H2CO3*DT H2CO3=H2CO3+D
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
537
H=H-D HCO3=HCO3-D K2M=K2
EXIT TEST IF E Q U I L I B R I U M R E A C H E D IF
(ABS(D).LT.TOL)
GO TO 18
INCR HEAT PRODUCED BY D E - I O N I Z A T I O N DQI=D*DELHI QI=QI+DQI T=T+DT 103
CONTINUE
18
CONTINUE
C
TOTAL HEAT Q(K)=QI+Q2 IF
(Q(K).GE.QMAX)
IF
(T.GE.0.2)
QMAX=Q(K)
GO TO 2
TI=TI+DTI 102
CONTINUE
2
CONTINUE
C C
EQUATION ROSSI-BERGER
PAPER J. Biol.
Chem.
243:1297-1302,1968.
P H = - A L O G I 0 (HI) DO 104 K = I , K M H = 1 0 * * (-PH) Y (K)=CI*ALOGI0 (HI/H)+C2*ALOGI0 ((HI+B-A) / (H+B-A)) T A U (K) =Y (K)/KH2CO3 I X = T A U (K)/TMAX*IXINT+IX0 IY= (PH-2)/5. *IYINT+IY0 PH=PH+. i C A L L DRWABS (IX, IY) IF 104
(TAU(K).GT.TMAX)
CONTINUE
GO TO 17
538 17
M A T H E M A T I C A L ANALYSIS AND MODELING
CONTINUE WRITE(5,921)TMAX
921
FORMAT('
',F5.3)
C TI=0. DO 121 K=I,KM IX=T1/TMAX*IXINT+IX0 IY=Q(K)/100.*IYINT+IY0 TI=TI+DT1 121
CONTINUE WRITE
922
(5,922)QMAX
FORMAT
(1H ,E13.3)
READ 501,KS 501
FORMAT
(I2)
IF (KS.NE.0)
GO TO 1
STOP END C
FORMATION OF H20 SUBROUTINE WATER
(H,OH,H20,KH20,KH,PH)
KWM=200 DO 106 KW=I,KWM
910
IF (PH.LT.7)
DTW=5.E-12*I0**PH
IF (PH.GE.7)
DTW=5.E-2*10**(-PH)
IF (KW.EQ.O)
WRITE(NFIL,910)
FORMAT
(IH ,4EI2.3,F6.2)
D=+KH*H*OH*DTW-KI420*H20*DTW H20=H20+D H=H-D OH=OH-D TW=TW+DTW 106
CONTINUE RETURN END
TW,H20,H,OH,PH
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
539
Appendix B: Program with Curve Fitter
****************************************************************** C C A C C E P T S D A T A IN P E R R E L L A S F O R M A T KINETICS OF HEMOGLOBIN-CO, ALPHA AND BETA CHAINS C A C C E P T S C O N C E N T R A T I O N OF UP T O 9 I N T E R M E D I A T E S
(B) A T
(QUENCHING)
C TIME C A C C E P T S S T A R T I N G V A L U E S OF i0 C H O S E N O N R A T E C O N S T A N T S C C A L C U L A T E S B E S T FIT F O R T H E R A T E C O N S T A N T S BY L E A S T S Q U A R E S C PROCEDURE C T H E D A T A A R E READ F R O M A S P E C I A L I N P U T FILE: C O . D A T C T H E S T A R T I N G V A L U E S A R E R E A D F R O M A N I N P U T FILE: H E X P I 0 . D A T C A L L C O N C S IN M O L S / L I T E R C
1 M/L
=
4 HEMES
C R A T E S IN L I T E R S / M O L / S E C C F I L E #4 IS T H E O U T P U T FILE
T H E F U N C T I O N A L V A L U E S F O R C U R V E - F I T A R E C R E A T E D BY S I M U L A T I N G T H E K I N E T I C S I T S E L F F R O M T=0 TO T=TC, N O T A M A T H FUNCTION. IMPLICIT DOUBLE PRECISION
(A-H,K,O-Z)
CHARACTER*40 INFILE,OFILE C H A R A C T E R * I IX DIMENSION PARM(30),IRK(30,3),RK(30) D I M E N S I O N X(15,15) COMMON NDATA,NPARM,YDAT(60),F(60),ERR(60),HB(20),CO0(20) COMMON PARTL(30,60),DT,TC(20),N,B(0:5,0:5,60),NC,N~SW,SF,ISWI COMMON GRAD(30),EI(60),A(20,20),C(20,20),P(30),ALAM,IT,S0,SI COMMON KSW,IPARM(30,3),IYDAT(30,3,20) COMMON RON(12),ROFF(12),NEXP,NBDAT N B D A T = C O U N T E R F O R T H E # O F B'S E N T E R E D N P A R M = # P A R M S TO BE F I T T E D NK
= C O U N T E R FOR T H E # OF R K ' S E N T E R E D
540
MATHEMATICAL ANALYSIS AND MODELING
C
N
= INDEX # OF EXPERIMENT
C
NEXP = # OF EXPERIMENTS
C
I
= # OF LIGANDED SITES
C
J
= # OF VALUES
C
N
= # OF EXPERIMENTS
C 250 C O N T I N U E C S E T V A R I A B L E S T O ZERO NPARM=O NBDAT=O NK=0 C
IN B(I,L,N)
I is t h e # of l i g a n d e d sites, N is t h e # of t h e
experiment. D O 2 I=l, 12
RON(I)=O. R O F F (I) =0. 2 CONTINUE D O 28 I=0,4 D O 28 L = 0 , 4 D O 28 N = 1 , 2 0
B(I,L,N)=0. 28 C O N T I N U E D O Ii J = l , 3 0 ERR(J)=O. F(J)=0. PARM(J) =0. Y D A T (J) =0. 11 C O N T I N U E D O 6 J=l, 16 DO 6 I=i,16 A(J, I)=O. C(J, I)=O.
6
CONTINUE
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
C C ********************************************************** 51
FORMAT
(/' I N P U T F I L E N A M E :
2001
FORMAT
(IA40)
2003
FORMAT
(' D A T A FILE:
502
FORMAT
(312,2X,EI5.5)
902
FORMAT
(' ',3i2)
900
FORMAT
(312)
61
FORMAT
(/' # O F E X P E R I M E N T S ,
504
FORMAT
(' N B D A T
503
FORMAT
(212,1X,EIO.4)
901
FORMAT
(2FI0.2)
903
FORMAT
(10(F6.4,1AI))
904
FORMAT
('
505
FORMAT
(9X,I2)
83
FORMAT
(F5.3)
1589
FORMAT
(E15.5)
C COMPUTE C
',IA40)
NEXP =
')
',I2)
',10(F6.4,1A1))
MODE
READ-IN
DATA FILE
5 PRINT READ
=
')
51 (*,2001)
INFILE
OPEN(UNIT=3,FILE=INFILE,STATUS='OLD') READ(3,2003) WRITE(5,2003)
ASSUMPTION: NK=
INFILE INFILE
FAST CONSTS
FIXED(NOT
PARMS)
12
NPARM=
i0
D O 168 K = I , N K R E A D (3,502) IRK(K, i), IRK(K, 2), IRK(K, 3) ,RON(K) WRITE(5,502) 168 C
CONTINUE ENTER
# OF B'S
IRK(K,I),IRK(K,2),IRK(K,3),RON(K)
541
542
MATHEMATICAL ANALYSIS AND MODELING READ
(3,900)
NBDAT
C A L L CLOSE(3)
E N T E R NO OF E X P E R I M E N T S PRINT
61
READ(*,900)
NEXP
C O N V E R T R A T E CONSTS TO PARAMS
C
PARM(1)=RON(1) PARM(2)=RON(4) PARM(3)=RON(8) PARM(4)=RON(2) PARM(5)=RON(5) PARM(6)=RON(9) PARM(7)=RON(3) PARM(8)=RON(6) PARM(9)=RON(7) PARM(10)=RON(10)
P R I N T 818 818
FORMAT READ IF
(' O U T P U T FILE,
(*,900)
YES(I)
OR NO(0)?
ISW1
(ISWI.EQ.0)
GO TO 122
P R I N T 814 814
FORMAT READ
2002
122
(*,2002)
FORMAT OPEN
(' O U T P U T FILENAME:
')
OFILE
(IA40)
(UNIT=4,FILE=OFILE,STATUS='NEW')
CONTINUE ISW=I
A C C E P T D A T A IN P E R R E L L A S FORMAT P R I N T 51
')
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS READ
(*,2001)
OPEN
(UNIT=3,FILE=INFILE,STATUS='OLD
543
INFILE ')
DO 1622 J=l,12 READ
(3,903)
(X(I,J),IX,I=I,10)
WRITE(5,904) 1622
(X(I,J),IX, I=I,10)
CONTINUE DO 120 N=I,10 Hb(N)=X(N,I) CO0(N)=X(N,2) TC(N)=X(N,3) WRITE(5,56)
Hb(N),CO0(N),TC(N)
IF(ISWI.GT.0) 56
WRITE(4,56)
FORMAT(/lX,'HB
C
ASSUMPTIONS
C
NO LIGANDS
= ',1PE15.5,/'
B(0,I,N)=X(N,4) C
O N E LIGAND B(I,I,N)=X(N,5) B(I,2,N)=X(N,6)
C
TWO LIGANDS B(2,I,N)=X(N,7) B(2,2,N)=0. B(2,3,N)=X(N,8) B(2,4,N)=X(N,9)
C
THREE LIGANDS B(3, I, N) =X(N, i0) B(3,2,N)=X(N,II)
C
HB(N),CO0(N),TC(N)
FOUR LIGANDS B(4,I,N)=X(N,12) I=0 L=I WRITE(5,158)I,L,N,B(I,L,N)
COO = ',E15.5/'
TC = ',E15.5)
544
MATHEMATICAL ANALYSIS AND MODELING IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
I=l DO 160 L=I,2 W R I T E (5,158) I,L,N, B (I,L,N) IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
160 CONTINUE I=2 DO 161 L=I,4 W R I T E (5,158) I, L, N, B (I, L, N) IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
161 CONTINUE I=3 DO 162 L=I,2 W R I T E (5,158) I, L,N, B (I, L,N) IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
162 CONTINUE I=4 L=I W R I T E (5,158) I, L,N, B (I, L, N) IF(ISWI.GT.0) 120
WRITE(4,158)
I,L,N,B(I,L,N)
CONTINUE CALL CLOSE (3)
DO 59 I=i,12 IF (RON(I).LT.IoE-15) WRITE(5,58) WRITE
GOTO 59
I,RON(I)
(4,58) I,RON(I)
59 CONTINUE 58 FORMAT('
KON(',I2,')
I' KOFF(',I2,')
= ',IPE15.5)
= ',1PE15.5)
55 CONTINUE 158 FORMAT('
B(',Ii,',',Ii,',',i2,')
= ',1PE15.5)
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
1588
FORMAT(2E15.5)
C C Here the B's for all the experiments C
and s t o r e d
are lumped together
in the array YDAT ready for curve-fit.
57 C O N T I N U E 95
format
(5i2)
N0=I DO 68 N=I,NEXP YDAT(0+N0)=B(0,I,N) YDAT(I+N0)=B(I,I,N) YDAT(2+N0)=B(I,2,N) YDAT(3+NO)=B(2,I,N) YDAT(4+N0)=B(2,3,N) YDAT(5+N0)=B(2,4,N) YDAT(6+N0)=B(3,I,N) YDAT(7+N0)=B(3,2,N) YDAT(8+N0)=B(4,I,N) WRITE(5,904) IF(ISWI.GT.0)
(YDAT(I+N0-1),IX,I=I,9) WRITE(4,904)
(YDAT(I+N0-1),IX,I=I,9)
N0=N0+9 68 C O N T I N U E
SCALE F A C T O R FOR MILLIMOLS
(KINETICS)
SF=I°E3 NDATA=NEXP*NBDAT C C 240
PROG R E S E T STARTS HERE CONTINUE PRINT
80
80
FORMAT READ IF
(' ENTER TIME STEP DT
(*,901)
DT
(ISWI. EQ. i) W R I T E (4,81) DT
(SEC):
')
545
546
MATHEMATICAL ANALYSIS AND MODELING
81
FORMAT
C
('/ TIME STEP DT
(SEC) =
[25]
',E15.2)
T H I S S W I T C H M A Y BE SET TO 1 TO PRINT THE K I N E T I C S NSW=O
C D E F I N E I N I T I A L V A L U E S OF P A R A M E T E R S A N D PASS TO C U R V E - F I T C A L L CFIT(PARM) C
CLOSE O U T P U T FILE IF(ISWI.EQ.I)
C A L L CLOSE(4)
P R I N T 35 35 F O R M A T READ
(i T Y P E 0 TO TERMINATE,
(*,900)
T Y P E 1 TO RUN N E W D A T A SET:
')
ISW3
IF(ISW3.EQ.I)
GO TO 250
STOP END C ************************************************************* SUBROUTINE
KINET(PARM)
IMPLICIT DOUBLE PRECISION
(A-H,K,O-Z)
D I M E N S I O N PARM(30) COMMON NDATA,NPARM,YDAT(60),F(60),ERR(60),HB(20),CO0(20) COMMON PARTL(30,60),DT,TC(20),N,B(0:5,0:5,60),NC,NSW,
COMMON G R A D ( 3 0 ) , E I ( 6 0 ) , A ( 2 0 ; 2 0 ) , C ( 2 0 , 2 0 ) , P ( 3 0 ) , A L A M , C O M M O N KSW, I P A R M ( 3 0 , 3 ) , I Y D A T ( 3 0 , 3 , 2 0 ) COMMON RON(12),ROFF(12),NEXP,NBDAT 904 C
FORMAT('
',10(F6.4,'
RESET TIMEAND
'))
ALL QUANTITIES
T=0. DO 2 I=0,4 DO 2 L=0,4 DO 2 N = I , 2 0 B(I,L,N)=0. 2 C
CONTINUE C O N V E R T P A R A M E T E R S B A C K TO RATE C O N S T A N T S RON(1)=PARM(1)
SF,ISW1
IT,S0,S1
[25]
FINITE ELEMENT NUMERICAL METHODS RON(4)=PARM(2) RON(8)=PARM(3) RON(2)=PARM(4) RON(5)=PARM(5) RON(9)=PARM(6) RON(3)=PARM(7) RON(6)=PARM(8) RON(7)=PARM(9) RON(10)=PARM(10)
36
FORMAT('
',I2,EIS.5,312,EIS.5)
C DO 500 N = I , N E X P C
SET INITIAL V A L U E S FOR U N L I G A N D E D H E M O G L O B I N AND CO B(0,I,N)=HB(N)
CO=CO0(N) T=0. C
KINETICS
SIMULATION
FROM T=0 TO T=TC
120 C O N T I N U E C C COMPUTE C ON R E A C T I O N S C FIRST L I G A T I O N
STEP
(from 0 to 1 ligand)
C B0 + CO -->
B(1,1)
ALPHA
C B0 + CO -->
B(I,2)
BETA
C
K(1)
AND K(2)
D=RON(1)*B(0,1,N)*CO*DT CO=CO-D
B(0,I,N)=B(0,I,N)-D B(I,I,N)=B(I,I,N)+D D=RON(2)*B(0,1,N)*CO*DT CO=CO-D B(0,I,N)=B(0,I,N)-D
SITE LIGANDED '
547
548
MATHEMATICAL ANALYSIS AND MODELING
B(I,2,N)=B(I,2,N)+D iii CONTINUE C C SECOND L I G A T I O N
STEP
C B(I,J)
+ CO --> B(2,J)
C
K(4),
K(3),
(from 1 to 2 ligands)
A N D K(5)
ONLY
D=RON(4)*B(1,I,N)*CO*DT
B(I,I,N)=B(I,I,N)-D B(2, l,N) =B (2, I,N) +D CO=CO-D D=RON(5)*B(1,2,N)*CO*DT B(I,2,N)=B(I,2,N)-D
B(2, I,N)=B(2, I,N) +D CO=CO-D D=RON(3)*B(I,1,N)*CO*DT CO=CO-D
B(I,I,N)=B(I,I,N)-D B(2,3,N)=B(2,3,N)+D D=RON(6)*B(1,2,N)*CO*DT CO=CO-D B(I,2,N)=B(I,2,N)-D B(2,4,N)=B(2,4,N) +D 113 C O N T I N U E C C THIRD LIGATION C B(2,J) C
K(7),
STEP
(from 2 to 3 ligands)
+ CO --> B(3,J) K(8),
A N D K(9)
ONLY
D=RON(8)*B(2,1,N)*CO*DT CO=CO-D B(2,I,N)=B(2,I,N)-D B(3,I,N)=B(3,I,N)+D D=RON(9)*B(2,1,N)*CO*DT
[05]
[25]
FINITE ELEMENT NUMERICAL METHODS CO=CO-D B(2,1,N) =B(2, I,N)-D B(3,2,N)=B(3,2,N)+D D=RON (7) *B(2,3 ,N) *CO*DT CO=CO-D B(2,3,N)=B(2,3,N)-D B(3, I,N)=B(3,1,N) +D D=RON (10) *B(2,4,N) *CO*DT CO=CO-D B(2,4,N)=B(2,4,N)-D B(3,2,N)=B(3,2,N) +D
114 CONTINUE
C FOURTH LIGATION C B(3,J)
STEP
(from 3 to 4 ligands)
+ CO --> B(4,J)
c K(11) ~ D
K(12)
D=RON (ii) *B (3,1, N) *CO*DT S(3,1,N)=B(3, I,N)-D CO=CO-D
B(4, I,N)=B(4, I,N) +D D=RON (12 ) *B (3,2, N) *CO*DT CO=CO-D B(3,2,N)=B(3,2,N)-D B(4, I,N)=B(4, I,N) +D 115 CONTINUE C C INCREMENT
T I M E COUNTER
T=T+DT IF(T.GE.TC(N) )GOTO 500 GO TO 120 500 CONTINUE
549
550
MATHEMATICAL ANALYSIS AND MODELING
501 C O N T I N U E C A L C U L A T E VALUES P R E D I C T E D BY K I N E T N0=I DO 68 N = I , N E X P F (0+N0) =B (0, I,N)
F(I+N0)=B(I,I,N) F (2+N0) =B (i, 2,N) F(3+N0) =B (2, l,N)
F(4+NO)=B(2,3,N) F(S+.0)=B(2,4,N) F (6+N0) =B(3, I,N) F(7+N0) =B(3,2,N)
F (8+NO) =B (4, I,N) N0=N0+9 901
FORMAT
(flO.4)
68 C O N T I N U E C C C
Calculate
errors b e t w e e n m e a s u r e d values
and the p r e d i c t e d values) DO 116 I = I , N D A T A
116
ERR(I) =F(I) -YDAT (I)
ii
FORMAT
(2F8.5)
RETURN END ******************************************* SUBROUTINE CURVE-FITS
CFIT(PARM) D A T A BY S I M U L A T I N G K I N E T I C S
MARQUARDT-LEVENBERGALGORITHM PARTIALS CALCULATED NUMERICALLY N D A T A = # OF O B S E R V A T I O N S PARM=
PARAMETERS
(DATA POINTS)
IN M O D E L
(the B's)
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
C
NPARM= # OF PARAMETERS
C
IT
C
NITER= MAX # OF ITERATIONS
C
ALAM = ITERATION PARAMETER
C
E, E R R = E R R O R
551
= ITERATION COUNTER
IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) DIMENSION PARM(30),DPARM(30),PNEW(30) COMMON NDATA,NPARM,YDAT(60),F(60),ERR(60),HB(20),CO0(20) COMMON PARTL(30,60),DT,TC(20),N,B(0:5,0:5,60),NC,NSW,SF,ISW1 COMMON GRAD(30),EI(60),A(20,20),C(20,20),P(30),ALAM,IT,S0,SI COMMON KSW,IPARM(30,3),IYDAT(30,3,20) COMMON RON(12),ROFF(12),NEXP,NBDAT REAL LAM,LAM0,ZI,Z2,S0,E,RSQ,SI,TOL P R I N T 810 810 F O R M A T READ
(' M A X # OF I T E R A T I O N S
')
(*,900) N I T E R
900
FORMAT
811
FORMAT
(312)
P R I N T 811
READ 901
(o I N I T I A L V A L U E O F LAMDA:
(*,901) A L A M
FORMAT(3F) P R I N T 816
816
FORMAT READ
(o L A M D A INCREMENT:
')
(*,901) D L A M
C INITIALIZE ITERATION COUNTER IT=0 ° C C
CALC SUM OF SQUARES OF ERRORS C A L L KINET(PARM) IF
(NITER.EQ.O) R E T U R N
S0=0. D O 1057 J = I , N D A T A
')
552
MATHEMATICAL ANALYSIS AND MODELING S 0 = S 0 + E R R (J) **2 El (J) = E R R (J)
1057
CONTINUE W R I T E (5,399) IF
(ISWI. EQ. i) W R I T E (4,399)
399
FORMAT('
ITER
SSE
PARAMS
C
S W I T C H FOR S P E C I A L D I S P L A Y S
')
(NOT USED)
NI=0
IF
(NI.EQ.0)
GO TO 230
P R I N T 813 813
FORMAT READ
(' NO. OF THE P A R A M TO BE SHOWN:
(*,900)
')
NP
C C ********************************************* 230
CONTINUE IT=IT+I
C
NEXT ITERATION
STARTS H E R E
DO 103 I = I , N P A R M GRAD (I) =0 103
CONTINUE DO 104 J = I , N P A R M DO 104 I = I , N P A R M A(I,J)=0. C(I,J)=0.
104
CONTINUE
C C C
CALCULATE NUMERICAL PARTIALS HERE DE/DPARM DON'T SHOW K I N E T I C S H E R E KSW=NSW NSW=0 DO 106 I = I , N P A R M D P A R M (I) =. 0 0 1 * P A R M (I)
[25]
[')5]
FINITE ELEMENT NUMERICAL METHODS PNEW (I) =PARM (I) +DPARM (I) PARM (I) =PARM (I) *i. 001 CALL KINET (PARM) P A R M (I) =PARM (I)/1. 001 DO 1058 J = I , N D A T A PARTL (I, J)=(El (J) -ERR(J) )/DPARM (I)
1058
CONTINUE
106
CONTINUE
ii
CONTINUE
C C
CALCULATE
GRADIENT
V E C T O R E*DE/DP
DO 1056 J = I , N D A T A CALC I-TH R O W OF JACOBIAN GRADIENT
VECTOR
DO 107 I=I,NPARM GRAD (I) =GRAD (I) +El (J) *PARTL (I, J) M A T R I X A IS COMPOSED
HERE
DO 108 II=I,NPARM A (I, Ii) =A (I, Ii) +PARTL (I, J) *PARTL (Ii, J) 108
CONTINUE
107
CONTINUE
1056
CONTINUE
C 365
CONTINUE DO 395 I=I,NPARM DO 390 J=I,NPARM
385
C(I,J)=A(I,J)
390
CONTINUE
395
CONTINUE ALAMI=I. +ALAM
OUTPUTS HERE
553
554 415
MATHEMATICAL
MODELING
[25]
CONTINUE IF
416
ANALYSIS AND
(NI.NE.0)
GO TO 417
CONTINUE WRITE(5,400) IF
IT,S0,(PARM(I),I=I,NPARM)
(ISWI.EQ.I)
400
FORMAT
417
CONTINUE
WRITE(4,400)
IT,S0,(PARM(I),I=I,NPARM)
(i ',I2,EII.4,12F8.1)
C C
IS P R E S E T NO OF I T E R A T I O N S IF (IT.EQ.NITER)
401
REACHED?
GO TO 687
CONTINUE
C C
PERTURB C-MATRIX DO iii I=I,NPARM
********************************************************************** DO 110 J = I , N P A R M ii0
C(I,J)=A(I,J)
**********************************************************************
C ( I , I ) = C ( I , I ) * ( i . +ALAM) 111
CONTINUE
C C
MATRIX
INVERSION HERE
C A L L MATIN4(C,NPARM)
CALCULATE
CHANGES IN PARAMS
DO 440 I = I , N P A R M P(I)=0. C
PRODUCT MATRIX
(DPARM)---(C)*(GRAD)
DO 115 J = I , N P A R M P (I) =P (I) +C (I, J) *GRAD (J) 115
CONTINUE
440
CONTINUE
[25]
FINITE ELEMENT NUMERICAL METHODS
IF
555
(NI.EQ.0) GO TO 418
WRITE(5,909)IT,PARM(NP),S0,GRAD(NP),C(NP,NP),ALAMI,P(NP) IF
(ISWI.EQ.I) W R I T E ( 4 , 9 0 9 ) I T , P A R M ( N P ) , S 0 , G R A D ( N P ) , C ( N P , N P ) ,
1ALAM1,P(NP) 909
FORMAT
(' ' , I 3 , F 7 . 0 , 3 E 1 5 . 5 , F 8 . 3 , F 8 . 0 )
418
CONTINUE
C C
UPDATE PARAMS
500
DO 501 I = I , N P A R M PARM(I)=PARM(I)+P(I)
501
CONTINUE
C C
CALCULATE NEW ERR
C
S H O W K I N E T I C S IF S P E C I F I E D NSW=KSW C A L L KINET(PARM) SNEW=0. DO 5211 J = I # N D A T A SNEW=SNEW+ERR(J)**2 EI(J)=ERR(J)
5211
CONTINUE
C C C C
T E S T IF S U M OF S Q U A R E S G E T T I N G L A R G E R IF
(SNEW.GT.S0)GOTO
NO, U P D A T E SSQ, GO A R O U N D I T E R A T I O N L O O P A G A I N ALAM=ALAM/DLAM S0=SNEW G O T O 230
C C 595
595
S N E W GT S O ALAM=ALAM*DLAM
556 C
MATHEMATICAL ANALYSIS AND MODELING
RESTORE
PARAMS TO PREV VALUES,
[25]
TRY A G A I N
CONTINUE DO 116 I=I,NPARM P A R M (I) =PARS(1) -P(I) 116
CONTINUE IT=IT+I GO TO 365
C **************************************************** C
ITERATIONS
635
CONTINUE
C
VARIANCE
DONE;
COMPUTE
FINAL Q U A N T I T I E S
OF THE FIT
V=S0/(NDATA-NPARM) VI=SQRT(V) WRITE(5,655)Vl IF 655
(ISWI.EQ.I)
WRITE(4,655)VI
FORMAT( 1' SD = ',E13.3) WRITE(5,675) IF
675
(ISWI.EQ.I)
FORMAT
WRITE(4,675)
(' FINAL PARAMS T E T R A M E R
DO 676 I=I,NPARM DO 676 J=I,NPARM 676
C(I,J)=A(I,J) CALL MATIN4(C,NPARM)
DO 123 I = I , N P A R M C
STAND.
ERR PARAM
D=VI*SQRT(C(I,I)) COEF=D/PARM(I) TET=PARM(I)/4. WRITE(5,685)PARM(I),TET,D,COEF
STD ERR PARAM
COEFF VAR')
[25]
557
FINITE ELEMENT NUMERICAL METHODS IF
(ISWI.EQ.I)
685
FORMAT
123
CONTINUE
WRITE(4,685)PARM(I),TET,D,COEF
(' ',F9.3,4X,F8.3,4X,2EI3.3)
C C A L L KINET(PARM) WRITE(5,780) IF(ISWI.EQ.I) 780
FORMAT
WRITE(4,780)
(' OBS NO
X
EST Y
ERR')
DO 125 J = I , N D A T A WRITE(5,815) IF 815
J,YDAT(J),F(J),ERR(J)
(ISWI.EQ.I)
FORMAT
WRITE(4,815)
J,YDAT(J),F(J),ERR(J)
(' ',I3,3FI0.5)
125 C O N T I N U E C C 700
SUBROUTINE
ENDS HERE FOR N O R M A L EXIT
RETURN
C C 687
P R E S E T ITERATIONS R E A C H E D CONTINUE PRINT 919
919
FORMAT READ IF
(' TYPE N O . F U R T H E R
ITERATIONS
DESIRED,OR
0 TO TERMINATE'
(*,900) MORE
(MORE.EQ.0)
GO TO 635
NITER=NITER+MORE GO TO 401 END C ********************************** S U B R O U T I N E MATIN4(ARRAY,NPARM) 3 INVERSION M A T R I X INVERSION BY G A U S S - J O R D A N NO P I V O T I N G IS USED ARRAY=MATRIX
TO BE INVERTED
N P A R M = SIZE OF M A T R I X TO BE INVERTED
ELIMINATION
558
MATHEMATICAL ANALYSIS AND MODELING ORIGINAL
IMPLICIT DIMENSION
MATRIX
[25]
IS DESTROYED
DOUBLE P R E C I S I O N
(A-H,K,O-Z)
ARRAY(20,20)
DO 604 I=I,NPARM STORE=ARRA¥(I,I) ARRAY(I,I)=I.0 DO 601 J = I , N P A R M IF 601
(STORE.EQ.0.0)
RETURN
~Y(I,J)=ARRAY(I,J)/STORE
C DO 604 IK=I,NPARM IF(IK-I) 602
602,604,602
STORE=ARRAY(IK,I) ARRAY(IK,I)=0.0 DO 603 J = I , N P A R M
603 A R R A Y (IK, J) =ARRAY (IK, J) - S T O R E * A R R A Y (I, J) C 604
CONTINUE RETURN END
Acknowledgments Acknowledgment is due to Stuart J. Davids for assistance in doing the computations in this chapter.
FINITE ELEMENT NUMERICAL METHODS
517
[25] S i m u l a t i o n o f H e m o g l o b i n K i n e t i c s U s i n g F i n i t e Element Numerical Methods
By
ROBERT L. BERGER, NORMAN DAVIDS,
and
M I C H E L E PERRELLA
Introduction Ordinary differential equations (ODEs) have been used for centuries as mathematical models of physical and chemical processes. The laws of mass action, for example, can be used to derive an ODE or system of ODEs describing a kinetic process. The advantage of using ODEs to express complex models arising from the laws of mass action is that they are easy to derive, even if they are difficult or impossible to solve analytically. Numerical solutions to such equations have been one of the great advantages of the computer age. The numerical analysis literature is full of solution methods, and a survey of that literature is beyond the scope of this chapter. The interested reader should see such topics as ordinary differential equations, stiff differential equations, initial value problems, and boundary value problems. In this chapter, we discuss one elementary method of analyzing biochemical kinetic problems to show what some of the issues are, including a discussion of how such methods appear to a user, in terms of ease of use. The general solution of a complex reaction, in this case hemoglobin (Hb) binding to carbon monoxide (CO), is then discussed in detail, including a global fit using a curve fitter. Finite Element Approach In simulating a rate process, we start from the law of mass action, which states that the rate of the forward reaction is proportional to the concentrations of the reactants, say A and B. If these are written as [A] and [B], Of = kf[A][B] and in a similar manner, vb = kb[C][D] if C and D are the products) In the simplest case we might have reactant A going to B, with the back reaction going so slowly that it can be neglected. Then we have
dA/dt = -kf[A]
(1)
where kf is the rate constant for the reaction and [A] is the concentration of the particular component. I W. J. Moore, "Physical Chemistry," 3rd ed., p. 169. Prentice-Hall, Englewood Cliffs, NJ, 1962.
METHODS IN ENZYMOLOGY,VOL. 232
Copyright © 1994by AcademicPress, inc. All rights of reproduction in any formreserved.
518
MATHEMATICAL ANALYSIS AND MODELING
[25]
Equation (1) may be interpreted in two ways: (1) that it is an ODE, and (2) that it is a statement about finite steps dA and dt. In the first case the ODE must be solved, which in this elementary case is easy, the solution in closed form being A = A0 e x p ( - k f t ) . We can, however, ignore this solution completely by refraining from having dA, dt pass to the limit, which is zero. In this way we bypass having to solve the differential equation either in closed form or numerically, because it is not generated in the first place. Thus Eq. (1) can be rewritten as dA = -kf[A] dt and used repetitively to provide an accounting of the increments of reaction products being created during a succession of time intervals dt. Following this step, we adjust or update the value of A by the relation [A]' = [A] + d[A]
(2)
These operations develop into a computer program when these steps are carried out for each reaction and each set of reactants, and include the p r o p e r i n p u t - o u t p u t statements and control statements. The particular scheme for doing this was developed by Davids and Berger 2,3 and is called the finite element or F E method. The above steps are inherent in the use of any digital procedure for simulating continuous processes: (I) The time must be divided into a discrete set of intervals so that the state of the system, stored as numerical values, can be updated repetitively, as described above; (2) the value of the time step dt and other numerical parameters must be chosen with care so that the computations are not only accurate but, even more fundamentally, remain stable; that is, do not cause a rapidly increasing error to build up. This point is discussed in more detail in the next section.
Stability Problem-Multiple Reaction Systems An inherent problem with any procedure that describes a physical system with finite sequential steps is that instability occurs when initially small errors grow, tending to build up large fluctuations in a system. Instability problems may arise when multiple reactions exist in a system having widely varying orders of magnitude in rate. As an illustration consider the following system: A + B C + D
kl k2
) C
(3)
, A + E
(4)
2 N. Davids and R. L. Berger, Commun. ACM 7, 547 (1964). 3 R. L. Berger, Discuss. Biophys. Soc., Biophys. J. 24, 2 (1978).
[25]
FINITE ELEMENT NUMERICAL METHODS
519
TABLE I COMPARISON OF EXACT SOLUTION VERSUS FINITE ELEMENT SOLUTION a
[A]I[Ao] kit
Analytical
Simulation
[B]/[B0]
0 0.2 0.4 0.6 0.8 1.0 2.0 3.0 4.0 5.0
1.000 0.819 0.670 0.549 0.449 0.368 0.135 0.050 0.018 0.007
1.000 0.818 0.669 0.547 0.448 0.366 0.134 0.049 0.018 0.007
0.000 0.164 0.268 0.331 0.363 0.374 0.302 0.204 0.144 0.114
a Comparison of the solutions of the kinetics of Eq. (5) according to the exact mathematical solution of Eqs. (6)-(8) and the finite element method simulation. Conditions: [A]0 = 1, kl = 1, k 3 = 1, k 2 = 0.1, and dt = 0.01 sec.
in which kl is several orders of magnitude faster than k2. In such a case the stability of the numerical model is maintained by choosing dtd for the fast reaction, say dt d = (1/100) dt, and then taking 100 steps for this reaction every time a single step of the slow reaction, which uses the time step dt, is calculated. We can think of this procedure as "synchronizing the clocks" associated with each of the reactions. Several kinetic cases have been taken from the literature to validate the procedure. As an example, consider the following case of two consecutive reactions: A
kl
> B
.
k3
" C
(5)
k2 A first, irreversible unimolecular reaction is followed by a reversible one. A general treatment of such reactions as given by Benson 4 and Johnson has given a two-exponential solution5 (extensive tables for this reaction 4 S. W. B e n s o n , " T h e Foundations of Chemical Kinetics." McGraw-Hill, New York, 1960 (see esp. Chapter III). 5 K. A. Johnson, this series, Vol. 134, p. 677.
520
MATHEMATICAL ANALYSISAND MODELING
[25]
based on exact m a t h e m a t i c a l formulas are available6). The rate equations for this reaction scheme are dA/dt = -kl[A] dB/dt = k~[A] + k2[C] - k3[B] dC/dt = -k2[C] + k3[B]
(6) (7) (8)
Table I gives a c o m p a r i s o n with the result c o m p u t e d f r o m the closed f o r m solution, taken f r o m K u b a and K o n o w a l o w 6 and from the F E method, assuming [A]0 = 1, kl = 1, k 3 = 1, and k2 = 0.1. The following is the F E p r o g r a m used, written in BASIC. The reaction of bicarbonate and hydrochloric acid is a well-known e x a m p l e o f this type of reaction scheme.7 On mixing, the reactants combine quickly (108 M -I sec -t) to f o r m carbonic acid. The b r e a k d o w n to carbon dioxide and w a t e r is also fast, 25 per second at 25 ° , and the reverse reaction is slow, e v e n at p H 7.4. A p p e n d i x A contains the Fortran program for solving this case b y the finite element method, but without the curve fitter. This requires making a first estimate o f the values of the rate constants. The calculation is then m a d e and a visual appraisal is made, comparing the F E output with the experimental data. F r o m the rates calculated, a new estimate can be m a d e and the process repeated until satisfactory a g r e e m e n t is reached. Appendix B shows a m o r e complicated case with the automatic c u r v e fitter added. This point is discussed in greater detail later on page 523.
Association Reaction between Hemoglobin and Carbon Monoxide Finite element analysis can be used to analyze the interaction between hemoglobin and CO in the general case, in which both the CO association and dissociation reactions occur, along with structural transition equilibria, during the course of the reaction. We consider here a special case in which substoichiometric a m o u n t s of CO are quantitatively bound by hemoglobin in a period of time during which the dissociation of the ligand is negligible. At the end of this preset reaction time, the liganded intermediates are q u e n c h e d and separated by cryogenic techniques. 8,9 The experimental data are in the f o r m of concentrations of the intermediates as a 6 D. W. Kuba and D. D. Konowalow, "Report ARL 63-44 (ASTIA)." University of Wisconsin, Madison, 1963. 7 L. Rossi-Bernardi and R. L. Berger, J. Biol. Chem. 243, 1297 (1968). 8 M. Perrella, L. Benazzi, L. Cremonesi, S. Vesely, G. Viggiano, and R. L. Berger, J. Biochem. Biophys. Methods 7, 187 (1983). 9 M. Perrella and L. Rossi-Bernardi, this volume [21].
[25]
FINITE ELEMENT NUMERICAL METHODS
521
30
'INITIALIZE:
40
A=I
'Initial
50
B=0
'
"
"
B
60
C=0
'
"
"
C
70
KI=I
'Rate const
for A - - - > B
80
K3=I
'Rate const
for B - - - > C
90
K2=.I
'Rate const
for b a c k reac B < - - - C
i00
T=0
ii0
DT=.01
'Set t i m e
120
KM=I01
'Set n u m b e r
130
'
conc
of A
step of
intervals
140
'START
REACTION
150
FOR K=I
TO KM
160
LPRINT
170
D=KI*A*DT
'Law o f M a s s
180
A=A-D
'Update
A
190
B=B+D
'Update
B
200
D=K3*B*DT-K2*C*DT
'Law o f M a s s
210
B=B-D
'Update
of B again
220
C=C+D
'Update
of C
230
T=T+DT
'Update
of time
240
NEXT
'End l o o p
250
STOP
K
USING
"#.##
###.###
#.###
#.###";;T,A,B,C Action
for A--->B
Action
for B<--->C
SCHEME I
function of fractional CO bound. Time information is not provided because at the high values of the hemoglobin concentration required by the cryogenic technique, the CO-binding reaction is so fast that sophisticated quench-flow methods are required to study the reaction on a time-dependent basis. Program Figure la is a general schematic showing the 10 ligation states and the pathways (16 in all) that can be taken from deoxyhemoglobin (Hb) to
522
MATHEMATICAL ANALYSIS AND MODELING
[25]
24
b
r"'-'-I k3
- 1231
k7
FIG. 1. (a) General schematic showing the 10 ligation states of hemoglobin and the 16 forward pathways leading from deoxyhemoglobin to fully liganded hemoglobin. Each ligation state is assigned indexes U, w h e r e j denotes the code of the isomer in ligation state i. Each rate constant k is assigned subscripts ijk: i, matrix from which transition originates;j, matrix element where transition terminates; k, number of sites liganded. (b) Simplified 12-rate constant path. Species 21 and 22 are grouped together because they cannot be resolved by the cryogenic technique for the isolation of the intermediates.
[25]
FINITE ELEMENT NUMERICAL METHODS
523
carbonmonoxyhemoglobin (HbCO). Each ligation state is located by the number of liganded sites (horizontal position). In a previous notation each state was also located by its position from the top of the diagram (columnwide). 3 We have changed this notation for the double-liganded species to conform to another such scheme now commonly used. l° A variable defining the rate constants is now introduced as (KON)iik, where i identifies the matrix element from which transition originates; j identifies the matrix element where transition terminates; and k is the number of sites liganded. Care must be taken in assigning statistical factors, 1 or 2, for the type of heme available for any given path. A value of zero for the variable corresponding to a particular transition in Fig. la prevents that transition from taking place. By putting in enough zeros, one can narrow the process to a single pathway or group of pathways. The general scheme for doing the simulation and optimizing the rate constants by global curve fitting is presented below, with applications to a simple fourstep model of CO hemoglobin binding and then using a 12-rate constant model (Fig. lb). Outline of Kinetics Curve-Fit Program Standard methods exist for determining the best fit by least squares to a set of data points, using a particular type of function, such as a polynomial or a series of exponentials, l~ The method of finite elements starts with the differential equations for the system in finite form. One then supplies the kinetics program with a set of trial parameters, that is, starting values or values determined from a previous iteration. We then go through the kinetics program and determine the error between the predicted curve and the data. These errors are then used to calculate the next mean-square error, which is compared with the previous one. A decision is then made as to whether to repeat the process with a disturbed set of parameters, until an optimum is reached. In this manner, it is possible to curve-fit a set of data by using kinetic models in place of numerical functions altogether. Figure 2 shows the scheme diagrammatically. 1. Read in and store the data, as identified by file name. 2. Enter the initial concentrations of the solutions, [Hb] and [CO]. 3. Enter the assumed transition sequence case. 4. Enter the assumed starting values for the rate constants. 5. Enter the number of iterations desired. l0 G. K. Ackers and F. R. Smith, Annu. Rev. Biophys. Biophys. Chem. 16, 583 (1987). 11 D. W. Marquardt, S I A M ] . 2, 431 (1963).
524
[25]
MATHEMATICAL ANALYSIS AND MODELING INITIAL VALUES of Rate Constants Bol, etc. Hb, CO for
~ 1
[
CURVE - FIT
I
[ [
~
Are these the best parameters?
Nfo-gl
I KINETICSl
FIG. 2. Outline of the finite element solution for the h e m o g l o b i n - C O reaction.
6. The c o m p u t e r then begins by passing the values of the given rate constants to the subroutine CFIT, where they are stored as the vector: parm;, i = 1, nparm. 7. These values are then passed to subroutine K I N E T , which works out a prediction function y(t) according to the transition sequence described there. The error is then calculated from e(t) = ydat(t) - y(t), where ydat(t) equals the observed data. 8. We then return to C F I T to obtain the sum of squares: s = X[e(t)] 2. 9. Each parameter is then perturbed in turn by the relation parm; = parm;(1 + 0.001). 10. With the above steps, the gradient vector can be set up normal to the surface o f a constant sum of squares in parameter space. The Jacobean matrix of the system can then be constructed. The inversion and diagonalization of this matrix, according to the G a u s s - N e w t o n procedure, provides an updated set of parameters and a new sum o f squares o f error, hopefully lower in magnitude. 11. The cycle of iterations terminates when the mean-square error converges. This can be determined by specifying a preset tolerance or the n u m b e r of iterations to be made. The user should review the display o f the successive sum of squares and then, by typing an integer, select additional iterations as desired.
Global Fit In general, any curve-fit scheme will yield more precise values for the parameters when more data points are available. In fact, for any kind of precision we must have the condition ndata >> nparams
(9)
[25]
FINITEELEMENTNUMERICALMETHODS
525
In the simple four-step scheme of the reaction, in which the subunit heterogeneity is ignored, there are five values of the concentration of the intermediates, B(i), available for determining four parameters, if time information is provided. If time information is not provided, there are three parameters relative to a fourth fixed parameter for making the determination. Because the sum of the concentrations of the intermediates at a certain value of the fractional CO bound equals the total hemoglobin concentration, only four values of the concentrations of the intermediates are independent. Thus for the simple four-step scheme it is possible to curve-fit the parameters by using the data from a single experiment at some fractional value of CO bound to hemoglobin. If data from experiments carried out at many different values of the fractional CO are available, a more reliable estimate of the error in the curve-fitted parameters is obtained by combining all of the data into what is called a "global" fit. This is done by simply "stacking the data" until we run out of data, as shown in Table II. Thus, if we have n = 10 experiments to work with in this fashion, we obtain 90 independent data points for the adjustable parameters relative to the reference parameter and condition (9) is fulfilled. If subunit functional heterogeneity is considered, I0 different ligation states are possible according to Fig. la. Because the cryogenic technique
TABLEII STACKINGMETHOD [BH] [B21] [B°~]t experiment1 [B31] [B41] [BI2] t [Bo2] [B22] experiment2
[B32] [B42]
[B.] I experimentn
526
MATHEMATICAL ANALYSIS AND MODELING
[25]
cannot determine the concentrations of intermediates B2~ and B22 separately, 9 values of the concentrations of the intermediates, B(i), are determined, of which only 8 values are independent, against a total of 12 parameters to curve-fit for, as shown in Fig. lb. Thus, even if k~l and k12 are assumed to be equal reference parameters, the number of adjustable parameters (10 values of the k) is greater than the number of data points, that is, 8 concentration values, and condition (9) is not fulfilled. In such a case the global fit using the data from several experiments at different values of fractional CO bound is mandatory. The Fortran program used to do this is shown in Appendix B.
Application to Simple Four-Step Scheme Because the functional heterogeneity of the subunits in the association reaction with CO is slight, 12the general process in Fig. la can be simplified by reducing it to a simple four-step process by putting kll I = kl21 k132 = k242
(10) (11)
kll 2 = k212 = k222 = k122
(12)
k313 = k423
(13) (14)
kll 3 -- k123 = k213 = k223 kll 4 = k214
(15)
Furthermore, when no time information is provided by the data, either the rate constants are calculated relative to a reference rate in the process or such a reference rate is assigned a value taken from the literature. Table III shows the results of the analyses of the experimental data, obtained by quench-flow and cryogenic isoelectric focusing methods s,9 to give the intermediates, to obtain the rate constants in the association reaction 12 between Hb and CO at 20°, in 0.1 M KC1, pH 7, assuming k4 = 6.0/zM -1 sec -l, as obtained by Gibson using optical stopped-flow techniques ~3 under similar conditions. The optimum time step dt was determined to be 0.01 msec. The analyses, starting with the data at the lowest fractional value of CO bound, were carried out by stepwise addition of the data at higher CO saturations to obtain the final global fit to all the
12 M. Perrella, N. Davids, and L. Rossi-Bernardi, J. Biol. Chem. 267, 8744 (1992). 13 Q. H. Gibson, J. Biol. Chem. 248, 1281 (1973).
[25]
FINITE ELEMENT NUMERICAL METHODS
527
TABLE III RATE CONSTANTS FOR REACTION OF HEMOGLOBIN PLUS
CO: FIVE-SPECIES MODELa
Fitted parameters (p.M -I sec -l)
Global fit up to experiment no.
Sco (%)
kl
k2
k3
SSE b (× 104)
I 2 3 4 5 6 7 8 9 10c
17.6 18.5 21.8 30.1 39.9 45.3 54.9 59.8 72.4 82.7
0.0777 0.0952 0.0944 0.0975 0.103 0.108 0.110 0.112 0.113 0.113
0.199 0.227 0.226 0.231 0.241 0.253 0.256 0.261 0.264 0.265
1.081 1.302 1.333 1.336 1.412 1.472 1.474 1.465 1.474 1.472
11.403 6.231 6.618 5.974 5.638 5.445 5.436 5.376 5.371 5.368
0.1024
0.242
1.382
6.286
Average valuesd:
a Rate constants for the association reaction between Hb and CO at 20°, pH 7, from data on the reaction intermediates as studied by cryogenic techniques. ~2Parameters (k 1, k2, k3) (intrinsic rate constants, i.e., apparent rate constants divided by the statistical factors 4, 3, and 2, respectively) were fitted by the finite element method as a function of the cumulative number of experiments at different CO saturations Sco, (1, 1 + 2 . . . . . 1 + 2 + -.- + 10), using k4 = 6.0/xM -1 sec-l. 12 b SSE, Sum of square errors (see text). c The error in the rate constants is - 10%. d Rate constants obtained by averaging the values of the parameters fitted to the data at each CO saturation value.
d a t a . T h e c a l c u l a t e d r a t e c o n s t a n t s at e a c h s t e p w e r e u s e d to fit t h e d a t a at all C O s a t u r a t i o n s ; t h e s u m o f s q u a r e e r r o r s f o r e a c h fit ( S S E ) a r e s h o w n in T a b l e I I I . Clearly, limiting values for the rate constants are reached inthe stepw i s e p r o c e s s a n d t h e m i n i m u m v a l u e f o r S S E is a t t a i n e d w i t h t h e p a r a m e t e r s f r o m t h e g l o b a l fit. F i g u r e 3 s h o w s t h e r e l a t i v e c o n c e n t r a t i o n s o f t h e s p e c i e s in d i f f e r e n t s t a t e s o f l i g a t i o n v e r s u s C O s a t u r a t i o n as c a l c u l a t e d b y t h e finite e l e m e n t method, assuming negligible functional heterogeneity of the subunits and u s i n g t h e v a l u e s o f t h e r a t e c o n s t a n t s o b t a i n e d b y t h e g l o b a l fit. If we now expand the possible pathways by utilizing the scheme of F i g . l b , f i x i n g t h e l a s t 2 r a t e c o n s t a n t s , k u = kl2 = 3.0 ~ M -1 s e c -1, a n d a g a i n d o t h e g l o b a l fit u s i n g all 10 e x p e r i m e n t s , w e o b t a i n t h e r e s u l t s s h o w n in T a b l e I V .
528
M A T H E M A T I C ANALYSIS AL AND MODELING
[25]
0.8
0
4
0.6
fi 0.4
0.2
0
0.2
0.4
0.6
0.8
Sco
FIc. 3. Fraction f/of liganded intermediate i (i = 0-4) versus CO saturation Sco in the association reaction between Hb and CO at 20° in 0.1 M KCI, pH 7) 2 The curves were calculated using the rate constants obtained by the global fit to the experimental data on the concentrations of the intermediates (see Table III), assuming negligiblefunctional heterogeneity of the subunits. T h e results o f this global fit to the CO-binding data, according to the 12-rate constants s c h e m e in Fig. Ib, confirm the finding of Perrella et a l ) 2 that the fl subunit reacts faster than the a subunit in the first CO-binding step (k2/k~ = 1.91 - 0.28). A similar subunit heterogeneity is o b s e r v e d b y comparing the rate constants for the second CO-binding step along the a a - ~ f l p a t h w a y (upper p a t h w a y s in Fig. Ib) and along the flfl-aa p a t h w a y (lower p a t h w a y s in Fig. lb), which yield the ratio k6/k 3 -- 2.98 - 1.76. Because species 21 and 22, both having one ct- and one fl-liganded subunit, are not equivalent and the cryogenic technique allows only the determination o f the sum of their concentrations, the rate constants k4, ks, k8, and k 9 are c o m p o s i t e constants. H o w e v e r , within experimental error, k 4 = k 5 and k 8 = k9. Thus the global fit confirms a previous conclusion that the CO-binding properties o f species 21 and 22 are similar ~2and that all binding steps are c o o p e r a t i v e with an acceleration after the binding of the second CO molecule. The simple model s h o w n above, and the m o r e complex 12-constant model, can be solved by analytical methods and the results are in close agreement. H o w e v e r , w h e n one goes to the reactions with o x y g e n and hemoglobin, analytical models of the 12 rate constants, including the b a c k reactions, which cannot be neglected as in the CO-binding case, b e c o m e
[25l
FINITE ELEMENT NUMERICAL METHODS
529
T A B L E IV RATE CONSTANTS FOR REACTION OF CO PLUS HEMOGLOBIN: NINE-SPECIES MODELa Ligation of ct subunits
Ligation of/3 subunits
kI k3 k5 k8 kj0
k2 = k0 = k6 = k7 = k9 =
= = = = =
0.0410 -+ 0.0042 0.0438 - 0.0215 0.153 -+ 0.054 0.454-+ 0.117 0.729 -+ 0.363
0.0782 -+ 0.0080 0.162 --- 0.03I 0.131 -+ 0.0428 0.133 +- 0.220 0.713 -+ 0.159
a Rate constants, k. (n = 1-12) b, for the binding of CO to deoxyhemoglobin at 20 ° in 0.1 M KCI, pH 7, calculated from the concentrations of the intermediates according to the reaction scheme of Fig. lb, assuming kll = k12 = 3.0 tzM -~ sec -~, that is, one-half the value measured by optical stopped-flow kinetics for the binding of the last CO moleculeJ 3 b Intrinsic, which is apparent constant k', divided by the statistical factor (in ~ M -j sec-1). The concentration is the heine concentration, kj = k l / 2 ; k 2 = kr/2; k 4 = k ' J 2 ; k 5 = k'5/2; ks = k~/2; k 9 = k ¥ 2 ; in all other cases k, = k,~.
intractable in general. We thus believe that the finite element method offers an important technique for solving or simulating large kinet~ problems. General
Conclusions
The FE method has the advantage that its steps are form,u'lated easily, as they deal directly with the physics or chemistry o f the problem. It simulates the problem in the sense that if the computation is stopped for any reason at an earlier time, it shows the state of the system at that time. Furthermore, it can easily deal with an additional boundary condition or the presence of additional effects, such as electric or magnetic fields, without invalidating the work up to that point. On the other hand, instabilities are difficult to deal with, just as with other numerical methods. In a more general sense, solving kinetic equations requires a number of methods. We have found that the program MATLAB 14 serves very well. We have presented a method of solving kinetic equations using any computer and a language such as Fortran, but we have found that the 14 MATLAB, PC or 386 version available from Mathworks, Inc., 24 Prime Park Way, Natick, MA 01760.
530
MATHEMATICAL ANALYSIS AND MODELING
[25]
utilization of MATLAB allows any of the numerical methods to be written down without concern of portability, as these programs will run on any version of MATLAB. These programs are still under development. It should be pointed out, however, that approaching a problem by numerical analysis of ODEs gives mathematical analysts the advantage of having at their disposal many established numerical and error analysis methods. They also have at their disposal such computer procedures as MATLAB, and so on, that contain within them the broad experience of many other analysts and that are universal in scope as well as convenience. Initial efforts to put the FE method into MATLAB and to utilize stiff differential equation solvers as well as curve fitters seem to be working and it is hoped these efforts will be reported in detail in the near future. 15 Appendix A: Program without Curve Fitter C
PROG:
C02.F4,
CARBON
DIOXIDE
KINETIC
MODELS
C C
REACTIONS
IN ACID
SOLUTIONS
C
COMPONENTS:
C
C02
= CARBON
C
H2CO3
= CARBONIC
C
HCO3-
= BICARBONATE
C
C03--
= CARBONATE
C
TOT
= TOTAL
C
PC02
= GASEOUS
C
KC02
= HYDRATION
C
KH2CO3=DEHYDRATION
C
HEATS:
DIOXIDE
(DISSOLVED)
ACID ION(-) ION
CARB.
(--),
DIOX.
HCO3=H+C03
CONTENT
F O R M O F C02, VELOCITY
OF SYSTEM
MEASURED
CONSTANT
VELOCITY
AS PARTIAL
PRES.
(THERMAL)
CONST.
C
DELHI
= HEAT
OF DE-IONIZATION
CAL/MOL
C
DELH2
= HEAT
OF DEHYDRATION
CAL/MOL
C
PROG.
CONTROL
C
KS=0,
TERMINATE
C
KS=I,
RECYCLE
DIMENSION
(=CO2+H2CO3+HC03-)
SWITCHES EXECUTION
PROGRAM
Y(500),TAU(500),Q(500)
15 R. I. Shrager, Division of Computer Research and Technology, National Institutes of
Health, Bethesda, MD (personal communication).
[25]
FINITE ELEMENT NUMERICAL METHODS
REAL*4 N A , N A H C O 3 , K C O 2 , K H 2 C O 3 , K C L , K H , K H 2 0 , K H C O 3 , K L 2 C 0 3 REAL*4 I S T R , M U C O E F F IST D I S S O C H2CO3
(KI)
A=3404.71 D=14.8435 C=0.032786 CONC OF H20
(MOLS/L)
H20=55.5 P R I N T 804 804
FORMAT
(' T E M P E R A T U R E O F E X P E R
(DEG-C):
R E A D 3, T E M P TABS=TEMP+273.16 H2CO3 DISSOC CONST PKI=A/TABS-D+C*TABS EI=I0**(-PKI) WRITE(21,911)TEMP,PKI,EI IF
(TEMP.EQ.37.) GO TO i0
ASSUMED RATE CONSTANTS
(/SEC):
D I S S O C D A T A F R O M LONG, BIOCHEM.HANDB. P43 D E H Y D R A T I O N C A R B O N I C ACID: H 2 C O 3 = C O 2 + H 2 0 VELOCITY CONST
(AT 25 DEG)
KH2C03=25.5 BACK REACTION RATE
(HYDRATION)
KCO2=.0257/H20 2ND D I S S O C C O N S T HCO3 PK=I0.329 A=2902.39 D=6.4980 C=0.02379 GO T O 20 C i0
CONTINUE
(K2), 25 D E G
'$)
531
532
MATHEMATICAL ANALYSIS AND MODELING VELOCITY
CONST
AT
37 D E G
KH2C03=72.00 KCO2=.0620/H20 20 C C
CONTINUE DISSOC
CONSTS
IONIZATION
C
EQUILIBR
C
PKWHEN
VS TEMP & FORMATION
CONST
H2CO3
CORRECTED
AT
FOR
OF CARBONIC 25 D E G
SALT
ACID:
H+HCO3=H2C03
IS 3 . 7 6 5
EFFECT
PK=3.5 PK AT
37 DEG:
PK=3.6 EH2CO3=I0**(-PK) ASSUMED
FORMATION
RATE
OF H2CO3
FROM
H+HCO3-
KHCO3=SEI0 DISSOCIATION
RATE
(IONIZATION)
OF H2CO3
KL2CO3=EH2CO3*KHCO3
FORMATION
OF H20:H+OH=H20
KH=I. 4Eli KH20=2.5E-5 DELHI=I205. DELH2=I455. SOLUBILITY
CONSTANTS
A=0. 0301 COMPUTER NFIL=23 KM=I00 K2=0 KPR = 1 VOL=I. 0 T=Oo 0 TI=0.0
INITIALIZATION
& CONTROL
CONSTS
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS TW=0.0 K=0 TW=0 QI=0. Q2=0. QMAX=0.
C
P R E P A R A T I O N OF S O L U T I O N
i: HCL
CALL A N M O D E WRITE(NFIL,912) 912
FORMAT
(' S O L U T I O N
1 (HCL)')
PRINT 801 801
FORMAT
(' M O L A R CONC SOLUTION 1 (HCL):
)
R E A D 3,HCL 3
FORMAT
(2F)
PRINT 802 802
FORMAT
(' CONC SALT ADDED:
)
R E A D 3,KCL C
H C L COMES TO E Q U I L I B R I U M
C
E Q U I L I B R I U M CONST H C L
(IONIZES):
HCL=H+CL
HI=HCL CLI=HCL PH0=-ALOGI0(H1) H=H1
PH=PH0 C
CHARGE B A L A N C E CHG=H-CL WRITE(NFIL,901)
C
PREP.
SOLUTION
K,K2,T,PH,HUM,OHUM,HCO3M,H2CO3M,CO2M,CHG
2:NAHCO3
PRINT 803 803
FORMAT
(' M O L A R CONC S O L U T I O N 2 ( N A H C O 3 ) : )
READ 3,NAHCO3 C
NAHCO3
COMES TO EQUILIBRIUM:
NAHCO3=NA+HCO3
533
534
MATHEMATICAL ANALYSIS AND MODELING NA=NAHCO3 HCO3=NAHC03 HCO3M=HCO3*IE3 WRITE(NFIL, 9 0 1 ) K , K 2 , T , P H , H U M
C 1 C
R E S T A R T OF C O M P U T A T I O N
HERE
CONTINUE M I X I N G A F T E R 1 MILLISEC,
NAHCO3+HCL=H2CO3+NACL
WRITE(NFIL, 903) 903
FORMAT
(IH ,'SOLUTIONS
1 & 2 B E C O M E MIXED')
CL=CL1/2.0 NA=HCO3/2.0 HCO3=HC03/2.0 H=HI/2.0 H2CO3=0.0 H2CO3M=H2CO3*IE3 HUM=H*I.E6*VOL OHUM=OH*I.E6*VOL C02M=CO2*IE3 CALL MOVABS(IX0,IY0) CO2=0.0 C
C A R B O N I C A C I D FORMED:
N+HCO3=H2C03
DT=5E-10 DO i01 K 2 = 1 , 4 0 D=KHCO3*H*HCO3*DT-KL2CO3*H2CO3*DT H2CO3=H2CO3+D H=H-D HCO3=HCO3-D C
EXIT AT E Q U I L I B R I U M TOL=0.1E-5 IF
(ABS(D).LT.TOL)
T=T+DT i01
CONTINUE
GO TO 13
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
13
CONTINUE
C
CHG BALANCE CHG=H+NA-CL-HCO3 HI=H
C
CONSTANTS
IN R O S S I - B E R N A R D I - B E R G E R
PAPER
A=.07 B=.I
CI=2.303*(I.+(EH2CO3/(B-A))) C2=2.303*(I.-(EH2CO3/(B-A))) C
IONIC STRENGTH AT START OF REAC
(M/L)
ISTR=0.1 C
DEBYE HUCKEL FORMULA MU=ISTR DB=0.505*SQRT(MU)/(I.+I.6*SQRT(MU))
C
MEAN IONIC ACTIVITY
COEFF
F=I0**(-DB) CI=2.303*(I.+(EH2CO3/F/(B-A))) C2=2.303*(I.-(EH2CO3/F/(B-A))) WRITE(21,911)CHG,DB,F,EH2CO3 911 C
FORMAT
(IH ,4E13.3)
KINETICS
LOOP HERE
DTI=.001 TMAX=KM*DTI DO 102 K=I,KM PH=-ALOGI0(N) CHG=H-HCO3-OH-CL+NA HCO3M=HCO3*IE3 N2CO3M=H2CO3*IE3 CO2M=CO2*IE3 IF((K-I)/EPR*KPR-(E-I))II,12,11 12
CONTINUE
535
536
MATHEMATICAL ANALYSIS AND MODELING P A R T I A L PRES.
CO2
(MMHG)
PCO2=(CO2M+H2CO3M) IA HUM=H*I.E6*VOL OHUM=OH*I.E6*VOL WRITE
(NFIL, 901) K , K 2 , T I , P H , H U M , O H U M , H C O 3 M , H 2 C O 3 M , C O 2 M , C H G
901
FORMAT(IH
ii
CONTINUE
,213,E10.2,F6.2,2X,2EI0.3,2X,3E10.3,2X,2EI0.3)
C
KINETICS
PLOT
C
PLOT OF PH VS T I M E IX=T1/TMAX*IXINT+IXO IY=(PH-2)/5.*IYINT+IY0 CALL DRWABS(IX,IY)
D E H Y D R A T I O N OF C A R B O N I C ACID:
H2CO3=H20+CO2
D=KH2CO3*H2CO3*DTI-KCO2*H20*CO2*DTI H2CO3=H2CO3-D C02=CO2+D H20=H20+D I N C R E M E N T A L HEAT P R O D U C E D BY D E H Y D R A T I O N DQ2=D*DELH2 Q2=Q2+DQ2 C
I O N I Z A T I O N OF W A T E R
C
H+OH=H20
C
CALL WATER(H,OH,H20,KH20,KH,PH)
C C 15
I O N I Z A T I O N OF C A R B O N I C A C I D
(H+)+(HCO3-)=H2CO3
DO 103 K 2 = I , 4 0 HUM=H*I.E6*VOL
C 902
WRITE
(NFIL,902)
FORMAT(IH
K2,H,HCO3,H2CO3,D
,I3,4E10.3)
D=KHCO3*H*HCO3*DT~KL2CO3*H2CO3*DT H2CO3=H2CO3+D
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
537
H=H-D HCO3=HCO3-D K2M=K2
EXIT TEST IF E Q U I L I B R I U M R E A C H E D IF
(ABS(D).LT.TOL)
GO TO 18
INCR HEAT PRODUCED BY D E - I O N I Z A T I O N DQI=D*DELHI QI=QI+DQI T=T+DT 103
CONTINUE
18
CONTINUE
C
TOTAL HEAT Q(K)=QI+Q2 IF
(Q(K).GE.QMAX)
IF
(T.GE.0.2)
QMAX=Q(K)
GO TO 2
TI=TI+DTI 102
CONTINUE
2
CONTINUE
C C
EQUATION ROSSI-BERGER
PAPER J. Biol.
Chem.
243:1297-1302,1968.
P H = - A L O G I 0 (HI) DO 104 K = I , K M H = 1 0 * * (-PH) Y (K)=CI*ALOGI0 (HI/H)+C2*ALOGI0 ((HI+B-A) / (H+B-A)) T A U (K) =Y (K)/KH2CO3 I X = T A U (K)/TMAX*IXINT+IX0 IY= (PH-2)/5. *IYINT+IY0 PH=PH+. i C A L L DRWABS (IX, IY) IF 104
(TAU(K).GT.TMAX)
CONTINUE
GO TO 17
538 17
M A T H E M A T I C A L ANALYSIS AND MODELING
CONTINUE WRITE(5,921)TMAX
921
FORMAT('
',F5.3)
C TI=0. DO 121 K=I,KM IX=T1/TMAX*IXINT+IX0 IY=Q(K)/100.*IYINT+IY0 TI=TI+DT1 121
CONTINUE WRITE
922
(5,922)QMAX
FORMAT
(1H ,E13.3)
READ 501,KS 501
FORMAT
(I2)
IF (KS.NE.0)
GO TO 1
STOP END C
FORMATION OF H20 SUBROUTINE WATER
(H,OH,H20,KH20,KH,PH)
KWM=200 DO 106 KW=I,KWM
910
IF (PH.LT.7)
DTW=5.E-12*I0**PH
IF (PH.GE.7)
DTW=5.E-2*10**(-PH)
IF (KW.EQ.O)
WRITE(NFIL,910)
FORMAT
(IH ,4EI2.3,F6.2)
D=+KH*H*OH*DTW-KI420*H20*DTW H20=H20+D H=H-D OH=OH-D TW=TW+DTW 106
CONTINUE RETURN END
TW,H20,H,OH,PH
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
539
Appendix B: Program with Curve Fitter
****************************************************************** C C A C C E P T S D A T A IN P E R R E L L A S F O R M A T KINETICS OF HEMOGLOBIN-CO, ALPHA AND BETA CHAINS C A C C E P T S C O N C E N T R A T I O N OF UP T O 9 I N T E R M E D I A T E S
(B) A T
(QUENCHING)
C TIME C A C C E P T S S T A R T I N G V A L U E S OF i0 C H O S E N O N R A T E C O N S T A N T S C C A L C U L A T E S B E S T FIT F O R T H E R A T E C O N S T A N T S BY L E A S T S Q U A R E S C PROCEDURE C T H E D A T A A R E READ F R O M A S P E C I A L I N P U T FILE: C O . D A T C T H E S T A R T I N G V A L U E S A R E R E A D F R O M A N I N P U T FILE: H E X P I 0 . D A T C A L L C O N C S IN M O L S / L I T E R C
1 M/L
=
4 HEMES
C R A T E S IN L I T E R S / M O L / S E C C F I L E #4 IS T H E O U T P U T FILE
T H E F U N C T I O N A L V A L U E S F O R C U R V E - F I T A R E C R E A T E D BY S I M U L A T I N G T H E K I N E T I C S I T S E L F F R O M T=0 TO T=TC, N O T A M A T H FUNCTION. IMPLICIT DOUBLE PRECISION
(A-H,K,O-Z)
CHARACTER*40 INFILE,OFILE C H A R A C T E R * I IX DIMENSION PARM(30),IRK(30,3),RK(30) D I M E N S I O N X(15,15) COMMON NDATA,NPARM,YDAT(60),F(60),ERR(60),HB(20),CO0(20) COMMON PARTL(30,60),DT,TC(20),N,B(0:5,0:5,60),NC,N~SW,SF,ISWI COMMON GRAD(30),EI(60),A(20,20),C(20,20),P(30),ALAM,IT,S0,SI COMMON KSW,IPARM(30,3),IYDAT(30,3,20) COMMON RON(12),ROFF(12),NEXP,NBDAT N B D A T = C O U N T E R F O R T H E # O F B'S E N T E R E D N P A R M = # P A R M S TO BE F I T T E D NK
= C O U N T E R FOR T H E # OF R K ' S E N T E R E D
540
MATHEMATICAL ANALYSIS AND MODELING
C
N
= INDEX # OF EXPERIMENT
C
NEXP = # OF EXPERIMENTS
C
I
= # OF LIGANDED SITES
C
J
= # OF VALUES
C
N
= # OF EXPERIMENTS
C 250 C O N T I N U E C S E T V A R I A B L E S T O ZERO NPARM=O NBDAT=O NK=0 C
IN B(I,L,N)
I is t h e # of l i g a n d e d sites, N is t h e # of t h e
experiment. D O 2 I=l, 12
RON(I)=O. R O F F (I) =0. 2 CONTINUE D O 28 I=0,4 D O 28 L = 0 , 4 D O 28 N = 1 , 2 0
B(I,L,N)=0. 28 C O N T I N U E D O Ii J = l , 3 0 ERR(J)=O. F(J)=0. PARM(J) =0. Y D A T (J) =0. 11 C O N T I N U E D O 6 J=l, 16 DO 6 I=i,16 A(J, I)=O. C(J, I)=O.
6
CONTINUE
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
C C ********************************************************** 51
FORMAT
(/' I N P U T F I L E N A M E :
2001
FORMAT
(IA40)
2003
FORMAT
(' D A T A FILE:
502
FORMAT
(312,2X,EI5.5)
902
FORMAT
(' ',3i2)
900
FORMAT
(312)
61
FORMAT
(/' # O F E X P E R I M E N T S ,
504
FORMAT
(' N B D A T
503
FORMAT
(212,1X,EIO.4)
901
FORMAT
(2FI0.2)
903
FORMAT
(10(F6.4,1AI))
904
FORMAT
('
505
FORMAT
(9X,I2)
83
FORMAT
(F5.3)
1589
FORMAT
(E15.5)
C COMPUTE C
',IA40)
NEXP =
')
',I2)
',10(F6.4,1A1))
MODE
READ-IN
DATA FILE
5 PRINT READ
=
')
51 (*,2001)
INFILE
OPEN(UNIT=3,FILE=INFILE,STATUS='OLD') READ(3,2003) WRITE(5,2003)
ASSUMPTION: NK=
INFILE INFILE
FAST CONSTS
FIXED(NOT
PARMS)
12
NPARM=
i0
D O 168 K = I , N K R E A D (3,502) IRK(K, i), IRK(K, 2), IRK(K, 3) ,RON(K) WRITE(5,502) 168 C
CONTINUE ENTER
# OF B'S
IRK(K,I),IRK(K,2),IRK(K,3),RON(K)
541
542
MATHEMATICAL ANALYSIS AND MODELING READ
(3,900)
NBDAT
C A L L CLOSE(3)
E N T E R NO OF E X P E R I M E N T S PRINT
61
READ(*,900)
NEXP
C O N V E R T R A T E CONSTS TO PARAMS
C
PARM(1)=RON(1) PARM(2)=RON(4) PARM(3)=RON(8) PARM(4)=RON(2) PARM(5)=RON(5) PARM(6)=RON(9) PARM(7)=RON(3) PARM(8)=RON(6) PARM(9)=RON(7) PARM(10)=RON(10)
P R I N T 818 818
FORMAT READ IF
(' O U T P U T FILE,
(*,900)
YES(I)
OR NO(0)?
ISW1
(ISWI.EQ.0)
GO TO 122
P R I N T 814 814
FORMAT READ
2002
122
(*,2002)
FORMAT OPEN
(' O U T P U T FILENAME:
')
OFILE
(IA40)
(UNIT=4,FILE=OFILE,STATUS='NEW')
CONTINUE ISW=I
A C C E P T D A T A IN P E R R E L L A S FORMAT P R I N T 51
')
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS READ
(*,2001)
OPEN
(UNIT=3,FILE=INFILE,STATUS='OLD
543
INFILE ')
DO 1622 J=l,12 READ
(3,903)
(X(I,J),IX,I=I,10)
WRITE(5,904) 1622
(X(I,J),IX, I=I,10)
CONTINUE DO 120 N=I,10 Hb(N)=X(N,I) CO0(N)=X(N,2) TC(N)=X(N,3) WRITE(5,56)
Hb(N),CO0(N),TC(N)
IF(ISWI.GT.0) 56
WRITE(4,56)
FORMAT(/lX,'HB
C
ASSUMPTIONS
C
NO LIGANDS
= ',1PE15.5,/'
B(0,I,N)=X(N,4) C
O N E LIGAND B(I,I,N)=X(N,5) B(I,2,N)=X(N,6)
C
TWO LIGANDS B(2,I,N)=X(N,7) B(2,2,N)=0. B(2,3,N)=X(N,8) B(2,4,N)=X(N,9)
C
THREE LIGANDS B(3, I, N) =X(N, i0) B(3,2,N)=X(N,II)
C
HB(N),CO0(N),TC(N)
FOUR LIGANDS B(4,I,N)=X(N,12) I=0 L=I WRITE(5,158)I,L,N,B(I,L,N)
COO = ',E15.5/'
TC = ',E15.5)
544
MATHEMATICAL ANALYSIS AND MODELING IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
I=l DO 160 L=I,2 W R I T E (5,158) I,L,N, B (I,L,N) IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
160 CONTINUE I=2 DO 161 L=I,4 W R I T E (5,158) I, L, N, B (I, L, N) IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
161 CONTINUE I=3 DO 162 L=I,2 W R I T E (5,158) I, L,N, B (I, L,N) IF(ISWI.GT.0)
WRITE(4,158)
I,L,N,B(I,L,N)
162 CONTINUE I=4 L=I W R I T E (5,158) I, L,N, B (I, L, N) IF(ISWI.GT.0) 120
WRITE(4,158)
I,L,N,B(I,L,N)
CONTINUE CALL CLOSE (3)
DO 59 I=i,12 IF (RON(I).LT.IoE-15) WRITE(5,58) WRITE
GOTO 59
I,RON(I)
(4,58) I,RON(I)
59 CONTINUE 58 FORMAT('
KON(',I2,')
I' KOFF(',I2,')
= ',IPE15.5)
= ',1PE15.5)
55 CONTINUE 158 FORMAT('
B(',Ii,',',Ii,',',i2,')
= ',1PE15.5)
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
1588
FORMAT(2E15.5)
C C Here the B's for all the experiments C
and s t o r e d
are lumped together
in the array YDAT ready for curve-fit.
57 C O N T I N U E 95
format
(5i2)
N0=I DO 68 N=I,NEXP YDAT(0+N0)=B(0,I,N) YDAT(I+N0)=B(I,I,N) YDAT(2+N0)=B(I,2,N) YDAT(3+NO)=B(2,I,N) YDAT(4+N0)=B(2,3,N) YDAT(5+N0)=B(2,4,N) YDAT(6+N0)=B(3,I,N) YDAT(7+N0)=B(3,2,N) YDAT(8+N0)=B(4,I,N) WRITE(5,904) IF(ISWI.GT.0)
(YDAT(I+N0-1),IX,I=I,9) WRITE(4,904)
(YDAT(I+N0-1),IX,I=I,9)
N0=N0+9 68 C O N T I N U E
SCALE F A C T O R FOR MILLIMOLS
(KINETICS)
SF=I°E3 NDATA=NEXP*NBDAT C C 240
PROG R E S E T STARTS HERE CONTINUE PRINT
80
80
FORMAT READ IF
(' ENTER TIME STEP DT
(*,901)
DT
(ISWI. EQ. i) W R I T E (4,81) DT
(SEC):
')
545
546
MATHEMATICAL ANALYSIS AND MODELING
81
FORMAT
C
('/ TIME STEP DT
(SEC) =
[25]
',E15.2)
T H I S S W I T C H M A Y BE SET TO 1 TO PRINT THE K I N E T I C S NSW=O
C D E F I N E I N I T I A L V A L U E S OF P A R A M E T E R S A N D PASS TO C U R V E - F I T C A L L CFIT(PARM) C
CLOSE O U T P U T FILE IF(ISWI.EQ.I)
C A L L CLOSE(4)
P R I N T 35 35 F O R M A T READ
(i T Y P E 0 TO TERMINATE,
(*,900)
T Y P E 1 TO RUN N E W D A T A SET:
')
ISW3
IF(ISW3.EQ.I)
GO TO 250
STOP END C ************************************************************* SUBROUTINE
KINET(PARM)
IMPLICIT DOUBLE PRECISION
(A-H,K,O-Z)
D I M E N S I O N PARM(30) COMMON NDATA,NPARM,YDAT(60),F(60),ERR(60),HB(20),CO0(20) COMMON PARTL(30,60),DT,TC(20),N,B(0:5,0:5,60),NC,NSW,
COMMON G R A D ( 3 0 ) , E I ( 6 0 ) , A ( 2 0 ; 2 0 ) , C ( 2 0 , 2 0 ) , P ( 3 0 ) , A L A M , C O M M O N KSW, I P A R M ( 3 0 , 3 ) , I Y D A T ( 3 0 , 3 , 2 0 ) COMMON RON(12),ROFF(12),NEXP,NBDAT 904 C
FORMAT('
',10(F6.4,'
RESET TIMEAND
'))
ALL QUANTITIES
T=0. DO 2 I=0,4 DO 2 L=0,4 DO 2 N = I , 2 0 B(I,L,N)=0. 2 C
CONTINUE C O N V E R T P A R A M E T E R S B A C K TO RATE C O N S T A N T S RON(1)=PARM(1)
SF,ISW1
IT,S0,S1
[25]
FINITE ELEMENT NUMERICAL METHODS RON(4)=PARM(2) RON(8)=PARM(3) RON(2)=PARM(4) RON(5)=PARM(5) RON(9)=PARM(6) RON(3)=PARM(7) RON(6)=PARM(8) RON(7)=PARM(9) RON(10)=PARM(10)
36
FORMAT('
',I2,EIS.5,312,EIS.5)
C DO 500 N = I , N E X P C
SET INITIAL V A L U E S FOR U N L I G A N D E D H E M O G L O B I N AND CO B(0,I,N)=HB(N)
CO=CO0(N) T=0. C
KINETICS
SIMULATION
FROM T=0 TO T=TC
120 C O N T I N U E C C COMPUTE C ON R E A C T I O N S C FIRST L I G A T I O N
STEP
(from 0 to 1 ligand)
C B0 + CO -->
B(1,1)
ALPHA
C B0 + CO -->
B(I,2)
BETA
C
K(1)
AND K(2)
D=RON(1)*B(0,1,N)*CO*DT CO=CO-D
B(0,I,N)=B(0,I,N)-D B(I,I,N)=B(I,I,N)+D D=RON(2)*B(0,1,N)*CO*DT CO=CO-D B(0,I,N)=B(0,I,N)-D
SITE LIGANDED '
547
548
MATHEMATICAL ANALYSIS AND MODELING
B(I,2,N)=B(I,2,N)+D iii CONTINUE C C SECOND L I G A T I O N
STEP
C B(I,J)
+ CO --> B(2,J)
C
K(4),
K(3),
(from 1 to 2 ligands)
A N D K(5)
ONLY
D=RON(4)*B(1,I,N)*CO*DT
B(I,I,N)=B(I,I,N)-D B(2, l,N) =B (2, I,N) +D CO=CO-D D=RON(5)*B(1,2,N)*CO*DT B(I,2,N)=B(I,2,N)-D
B(2, I,N)=B(2, I,N) +D CO=CO-D D=RON(3)*B(I,1,N)*CO*DT CO=CO-D
B(I,I,N)=B(I,I,N)-D B(2,3,N)=B(2,3,N)+D D=RON(6)*B(1,2,N)*CO*DT CO=CO-D B(I,2,N)=B(I,2,N)-D B(2,4,N)=B(2,4,N) +D 113 C O N T I N U E C C THIRD LIGATION C B(2,J) C
K(7),
STEP
(from 2 to 3 ligands)
+ CO --> B(3,J) K(8),
A N D K(9)
ONLY
D=RON(8)*B(2,1,N)*CO*DT CO=CO-D B(2,I,N)=B(2,I,N)-D B(3,I,N)=B(3,I,N)+D D=RON(9)*B(2,1,N)*CO*DT
[05]
[25]
FINITE ELEMENT NUMERICAL METHODS CO=CO-D B(2,1,N) =B(2, I,N)-D B(3,2,N)=B(3,2,N)+D D=RON (7) *B(2,3 ,N) *CO*DT CO=CO-D B(2,3,N)=B(2,3,N)-D B(3, I,N)=B(3,1,N) +D D=RON (10) *B(2,4,N) *CO*DT CO=CO-D B(2,4,N)=B(2,4,N)-D B(3,2,N)=B(3,2,N) +D
114 CONTINUE
C FOURTH LIGATION C B(3,J)
STEP
(from 3 to 4 ligands)
+ CO --> B(4,J)
c K(11) ~ D
K(12)
D=RON (ii) *B (3,1, N) *CO*DT S(3,1,N)=B(3, I,N)-D CO=CO-D
B(4, I,N)=B(4, I,N) +D D=RON (12 ) *B (3,2, N) *CO*DT CO=CO-D B(3,2,N)=B(3,2,N)-D B(4, I,N)=B(4, I,N) +D 115 CONTINUE C C INCREMENT
T I M E COUNTER
T=T+DT IF(T.GE.TC(N) )GOTO 500 GO TO 120 500 CONTINUE
549
550
MATHEMATICAL ANALYSIS AND MODELING
501 C O N T I N U E C A L C U L A T E VALUES P R E D I C T E D BY K I N E T N0=I DO 68 N = I , N E X P F (0+N0) =B (0, I,N)
F(I+N0)=B(I,I,N) F (2+N0) =B (i, 2,N) F(3+N0) =B (2, l,N)
F(4+NO)=B(2,3,N) F(S+.0)=B(2,4,N) F (6+N0) =B(3, I,N) F(7+N0) =B(3,2,N)
F (8+NO) =B (4, I,N) N0=N0+9 901
FORMAT
(flO.4)
68 C O N T I N U E C C C
Calculate
errors b e t w e e n m e a s u r e d values
and the p r e d i c t e d values) DO 116 I = I , N D A T A
116
ERR(I) =F(I) -YDAT (I)
ii
FORMAT
(2F8.5)
RETURN END ******************************************* SUBROUTINE CURVE-FITS
CFIT(PARM) D A T A BY S I M U L A T I N G K I N E T I C S
MARQUARDT-LEVENBERGALGORITHM PARTIALS CALCULATED NUMERICALLY N D A T A = # OF O B S E R V A T I O N S PARM=
PARAMETERS
(DATA POINTS)
IN M O D E L
(the B's)
[25]
[25]
FINITE ELEMENT NUMERICAL METHODS
C
NPARM= # OF PARAMETERS
C
IT
C
NITER= MAX # OF ITERATIONS
C
ALAM = ITERATION PARAMETER
C
E, E R R = E R R O R
551
= ITERATION COUNTER
IMPLICIT DOUBLE PRECISION(A-H,K,O-Z) DIMENSION PARM(30),DPARM(30),PNEW(30) COMMON NDATA,NPARM,YDAT(60),F(60),ERR(60),HB(20),CO0(20) COMMON PARTL(30,60),DT,TC(20),N,B(0:5,0:5,60),NC,NSW,SF,ISW1 COMMON GRAD(30),EI(60),A(20,20),C(20,20),P(30),ALAM,IT,S0,SI COMMON KSW,IPARM(30,3),IYDAT(30,3,20) COMMON RON(12),ROFF(12),NEXP,NBDAT REAL LAM,LAM0,ZI,Z2,S0,E,RSQ,SI,TOL P R I N T 810 810 F O R M A T READ
(' M A X # OF I T E R A T I O N S
')
(*,900) N I T E R
900
FORMAT
811
FORMAT
(312)
P R I N T 811
READ 901
(o I N I T I A L V A L U E O F LAMDA:
(*,901) A L A M
FORMAT(3F) P R I N T 816
816
FORMAT READ
(o L A M D A INCREMENT:
')
(*,901) D L A M
C INITIALIZE ITERATION COUNTER IT=0 ° C C
CALC SUM OF SQUARES OF ERRORS C A L L KINET(PARM) IF
(NITER.EQ.O) R E T U R N
S0=0. D O 1057 J = I , N D A T A
')
552
MATHEMATICAL ANALYSIS AND MODELING S 0 = S 0 + E R R (J) **2 El (J) = E R R (J)
1057
CONTINUE W R I T E (5,399) IF
(ISWI. EQ. i) W R I T E (4,399)
399
FORMAT('
ITER
SSE
PARAMS
C
S W I T C H FOR S P E C I A L D I S P L A Y S
')
(NOT USED)
NI=0
IF
(NI.EQ.0)
GO TO 230
P R I N T 813 813
FORMAT READ
(' NO. OF THE P A R A M TO BE SHOWN:
(*,900)
')
NP
C C ********************************************* 230
CONTINUE IT=IT+I
C
NEXT ITERATION
STARTS H E R E
DO 103 I = I , N P A R M GRAD (I) =0 103
CONTINUE DO 104 J = I , N P A R M DO 104 I = I , N P A R M A(I,J)=0. C(I,J)=0.
104
CONTINUE
C C C
CALCULATE NUMERICAL PARTIALS HERE DE/DPARM DON'T SHOW K I N E T I C S H E R E KSW=NSW NSW=0 DO 106 I = I , N P A R M D P A R M (I) =. 0 0 1 * P A R M (I)
[25]
[')5]
FINITE ELEMENT NUMERICAL METHODS PNEW (I) =PARM (I) +DPARM (I) PARM (I) =PARM (I) *i. 001 CALL KINET (PARM) P A R M (I) =PARM (I)/1. 001 DO 1058 J = I , N D A T A PARTL (I, J)=(El (J) -ERR(J) )/DPARM (I)
1058
CONTINUE
106
CONTINUE
ii
CONTINUE
C C
CALCULATE
GRADIENT
V E C T O R E*DE/DP
DO 1056 J = I , N D A T A CALC I-TH R O W OF JACOBIAN GRADIENT
VECTOR
DO 107 I=I,NPARM GRAD (I) =GRAD (I) +El (J) *PARTL (I, J) M A T R I X A IS COMPOSED
HERE
DO 108 II=I,NPARM A (I, Ii) =A (I, Ii) +PARTL (I, J) *PARTL (Ii, J) 108
CONTINUE
107
CONTINUE
1056
CONTINUE
C 365
CONTINUE DO 395 I=I,NPARM DO 390 J=I,NPARM
385
C(I,J)=A(I,J)
390
CONTINUE
395
CONTINUE ALAMI=I. +ALAM
OUTPUTS HERE
553
554 415
MATHEMATICAL
MODELING
[25]
CONTINUE IF
416
ANALYSIS AND
(NI.NE.0)
GO TO 417
CONTINUE WRITE(5,400) IF
IT,S0,(PARM(I),I=I,NPARM)
(ISWI.EQ.I)
400
FORMAT
417
CONTINUE
WRITE(4,400)
IT,S0,(PARM(I),I=I,NPARM)
(i ',I2,EII.4,12F8.1)
C C
IS P R E S E T NO OF I T E R A T I O N S IF (IT.EQ.NITER)
401
REACHED?
GO TO 687
CONTINUE
C C
PERTURB C-MATRIX DO iii I=I,NPARM
********************************************************************** DO 110 J = I , N P A R M ii0
C(I,J)=A(I,J)
**********************************************************************
C ( I , I ) = C ( I , I ) * ( i . +ALAM) 111
CONTINUE
C C
MATRIX
INVERSION HERE
C A L L MATIN4(C,NPARM)
CALCULATE
CHANGES IN PARAMS
DO 440 I = I , N P A R M P(I)=0. C
PRODUCT MATRIX
(DPARM)---(C)*(GRAD)
DO 115 J = I , N P A R M P (I) =P (I) +C (I, J) *GRAD (J) 115
CONTINUE
440
CONTINUE
[25]
FINITE ELEMENT NUMERICAL METHODS
IF
555
(NI.EQ.0) GO TO 418
WRITE(5,909)IT,PARM(NP),S0,GRAD(NP),C(NP,NP),ALAMI,P(NP) IF
(ISWI.EQ.I) W R I T E ( 4 , 9 0 9 ) I T , P A R M ( N P ) , S 0 , G R A D ( N P ) , C ( N P , N P ) ,
1ALAM1,P(NP) 909
FORMAT
(' ' , I 3 , F 7 . 0 , 3 E 1 5 . 5 , F 8 . 3 , F 8 . 0 )
418
CONTINUE
C C
UPDATE PARAMS
500
DO 501 I = I , N P A R M PARM(I)=PARM(I)+P(I)
501
CONTINUE
C C
CALCULATE NEW ERR
C
S H O W K I N E T I C S IF S P E C I F I E D NSW=KSW C A L L KINET(PARM) SNEW=0. DO 5211 J = I # N D A T A SNEW=SNEW+ERR(J)**2 EI(J)=ERR(J)
5211
CONTINUE
C C C C
T E S T IF S U M OF S Q U A R E S G E T T I N G L A R G E R IF
(SNEW.GT.S0)GOTO
NO, U P D A T E SSQ, GO A R O U N D I T E R A T I O N L O O P A G A I N ALAM=ALAM/DLAM S0=SNEW G O T O 230
C C 595
595
S N E W GT S O ALAM=ALAM*DLAM
556 C
MATHEMATICAL ANALYSIS AND MODELING
RESTORE
PARAMS TO PREV VALUES,
[25]
TRY A G A I N
CONTINUE DO 116 I=I,NPARM P A R M (I) =PARS(1) -P(I) 116
CONTINUE IT=IT+I GO TO 365
C **************************************************** C
ITERATIONS
635
CONTINUE
C
VARIANCE
DONE;
COMPUTE
FINAL Q U A N T I T I E S
OF THE FIT
V=S0/(NDATA-NPARM) VI=SQRT(V) WRITE(5,655)Vl IF 655
(ISWI.EQ.I)
WRITE(4,655)VI
FORMAT( 1' SD = ',E13.3) WRITE(5,675) IF
675
(ISWI.EQ.I)
FORMAT
WRITE(4,675)
(' FINAL PARAMS T E T R A M E R
DO 676 I=I,NPARM DO 676 J=I,NPARM 676
C(I,J)=A(I,J) CALL MATIN4(C,NPARM)
DO 123 I = I , N P A R M C
STAND.
ERR PARAM
D=VI*SQRT(C(I,I)) COEF=D/PARM(I) TET=PARM(I)/4. WRITE(5,685)PARM(I),TET,D,COEF
STD ERR PARAM
COEFF VAR')
[25]
557
FINITE ELEMENT NUMERICAL METHODS IF
(ISWI.EQ.I)
685
FORMAT
123
CONTINUE
WRITE(4,685)PARM(I),TET,D,COEF
(' ',F9.3,4X,F8.3,4X,2EI3.3)
C C A L L KINET(PARM) WRITE(5,780) IF(ISWI.EQ.I) 780
FORMAT
WRITE(4,780)
(' OBS NO
X
EST Y
ERR')
DO 125 J = I , N D A T A WRITE(5,815) IF 815
J,YDAT(J),F(J),ERR(J)
(ISWI.EQ.I)
FORMAT
WRITE(4,815)
J,YDAT(J),F(J),ERR(J)
(' ',I3,3FI0.5)
125 C O N T I N U E C C 700
SUBROUTINE
ENDS HERE FOR N O R M A L EXIT
RETURN
C C 687
P R E S E T ITERATIONS R E A C H E D CONTINUE PRINT 919
919
FORMAT READ IF
(' TYPE N O . F U R T H E R
ITERATIONS
DESIRED,OR
0 TO TERMINATE'
(*,900) MORE
(MORE.EQ.0)
GO TO 635
NITER=NITER+MORE GO TO 401 END C ********************************** S U B R O U T I N E MATIN4(ARRAY,NPARM) 3 INVERSION M A T R I X INVERSION BY G A U S S - J O R D A N NO P I V O T I N G IS USED ARRAY=MATRIX
TO BE INVERTED
N P A R M = SIZE OF M A T R I X TO BE INVERTED
ELIMINATION
558
MATHEMATICAL ANALYSIS AND MODELING ORIGINAL
IMPLICIT DIMENSION
MATRIX
[25]
IS DESTROYED
DOUBLE P R E C I S I O N
(A-H,K,O-Z)
ARRAY(20,20)
DO 604 I=I,NPARM STORE=ARRA¥(I,I) ARRAY(I,I)=I.0 DO 601 J = I , N P A R M IF 601
(STORE.EQ.0.0)
RETURN
~Y(I,J)=ARRAY(I,J)/STORE
C DO 604 IK=I,NPARM IF(IK-I) 602
602,604,602
STORE=ARRAY(IK,I) ARRAY(IK,I)=0.0 DO 603 J = I , N P A R M
603 A R R A Y (IK, J) =ARRAY (IK, J) - S T O R E * A R R A Y (I, J) C 604
CONTINUE RETURN END
Acknowledgments Acknowledgment is due to Stuart J. Davids for assistance in doing the computations in this chapter.