- Email: [email protected]
Theoretical and experimental studies of the electron transfer reactions in the cytochrome chain
ARCHIVES
OF BIOCHEMISTRY
Theoretical
and
AND
BIOPHYSICS
147, 675-682 (1971)
Experimental
Reactions
Studies
of the Electron
in the Cytochrome
II. Fitting of the Experimental MICHAEL Johnson
WAGNER,
MARIA
ERECI@XA
Transfer
Chain’ Data AND
MARTIN
PRING
University of Pennsylvania, Philadelphia, Pennsylvania 19104 Received August 30, 1971;accepted September 17, 1971
Research Fomdation,
The experimental curves of cytochromes a, c, and cl were fitted to the kinetic model presented in the first part of this study. Quantitative agreement between the experiments and the model has been found. Estimation of the initial conditions for computer treatment of the cytochrome kinetics, the applicability of the nonlinear optimization program to the biological sciences and the determination of the errors with respect to the constants are presented.
In contrast to modern physics, there are still no explicit means in the mathematical language of chemical kinetics to relate the observer and the measuring process to the analysis of the results. Thus the only way to establish the validity of a mode, is to compare the experimental curves with those derived theoretically. But even then it is diacult to say how far the agreement between theory and experiments confirms the model. However, a qualitative treatment based on accurately estimated parameters can lead to certain predictions in the experimental field. This is a positive manner in which mathematical tools can be utilized in modern biology. The need for continuous critical comparison of the simulated models and experiments in the process of searching for a reaction mechanism is also suggested by the absence of adequate mathematical language for the evaluation of biochemical 1 Supported by USPH GM Computer time was supported
12202 and 16501. by NIH Grant
#RR-15. All the programs used in this work are available upon request. The Davidon optimization program was modified for biochemical purposes, and the plotting programs are written in FORTRAN IV; the integrating program is in PDP-B-10 MACRO assembly language.
kinetics. There is no commonly valid criterion leading to the acceptance of one model as “true” or “convenient”2 and the rejection of another as “false” or “inconThe criteria of mathematical venient”. statistics and testing of hypotheses might sometimes be applied (l), but the most important criterion is the “predictive of the model for other indeproperty” pendent experiments, independent at least in the mathematical sense if not in practical organization of the experiments. In the first part of this study (3) we have shown qualitative agreement between the effect of inhibitors of the cytochrome system and a kinetic model able to simulate most of the important responses. We will present quantitative results obtained by using a nonlinear optimization technique to determine velocity constants and compare more closely the model presented in the first part of this study with the experimental behavior of the cytochrome system. This work is the result of fitting experimental curves of cytochrome oxidation in the presence of two different concentra2 The term convenient is normally used in more abstract sciences (a), but is equally justified in biochemistry.
675
WAGNER, TABLE DETERMINATION PARSMETERS EXPERIMENT BLOCK'"*~
ERECIP;TSKA,
I
OF THE ACCURACY OF FOR THE FINAL FIT FROM WITH HIGH ANTIMYCIN
THE THE
A
Velocity constant Nbv-
Value
1 2 3 4 5 6 7
5.8 221 65 28 20 0.18 0.02
E
0.67 1.51 0.99 1.54 1.14 2.53 4.97b
A(+40%)
0.11 0.0088 0.0097 0.016
-4
0.023 x 10-b
A(-40%)
0.49 0.026 0.016 0.077 0.017 0.015 -4 x 10-s
a In column 3 the accuracy of the parameters is estimated from the values of the diagonal members of the Hessian matrices. In columns 4 and 5 the accuracy of the parameters is estimated from the changes of S&R if the constants are changed by 40’%. A =
(S&R&
+
40%)-SQR(kj))
The dimensionless numbers correspond to the E = Hi,i/ki and are interpreted as relative error. b Some nondiagonal members of matrix still negative. ’ Reaction scheme and S&R identical with those in table II.
tions of antimycin A. It represents the situation in which the numerical values of all constants except that of the reaction No. 5 of Table I are predicted by one experiment (in the presence of 0.12 pg antimycin A/mg prot. A) and the accuracy of quantitative of prediction for the second experimental (in the presence of 0.06 kg of antimycin A/mg prot.) is tested. The results of previous fittings of the experimental curves of cytochrome kinetics (H), lead to the detection of several inaccuracies in the assumptions and experimental approach to the study of oxidation in the cytochrome system. Therefore, an attempt has been made to obtain a fit of the experimental data with high accuracy with respect to the sum of squares of deviations as well as to the parameters.a The results presented here are intended to a The term parameters means unknown variables, experimentally unavailable. In our case. these are the values of the velocity constants (which are never determined by direct experiment) and the concentration of the substance “POOL”.
AND
PRING
increase the accuracy of previous fits. Fitting the experimental data by means of our optimization programs led us to an understanding of the limits of the present model. In addition, optimization produced a better estimate of the velocity constants to be used for simulation studies with extended models. METHODS Experimental. The experimental data were obtained in a stopped flow apparatus using pigeon heart mitochondria, as described in the first paper of this series. A low concentration of antimycin A was added to the system inhibited with rotenone, TTFB and CO. The amount of antimycin A was doubled and the traces were retaken. Since the fitting process is much more sensitive to experimental inaccuracies than a simulation study, precautions were taken to ensure the necessary high reproducibility of the experimental data. Oxygen was delivered to the anaerobic mitochondrial mixture at nearly identical intervals (5 min) and the experimental traces were taken in duplicate. The choice of model. The most difficult and vulnerable part of any study which is neither purely experimental nor purely theoretical is the relation between the experiments and the model used to interpret the results, i.e., t.he choice of some model and the translation of one model to another. For this reason, and because the individual steps in our reaction scheme might seem arbitrary without a more extensive discussion, we describe here in some detail the design of a model used in this work. The reactions considered are those concerning the experimentally measured cytochromes a, c, cl . An attempt was made to approximate the part of the reaction scheme concerning the cytochromes not discussed in this work by as few reaction steps and variables aa possible. The system of reactions transferring electrons from that part of the cytochrome system studied here to oxygen (via CO uninhibited cytochromes a-a3 can be approximated under our experimental conditions as one first order reaction. This is possible because in differential equations for all reactions between a3 (red) + oxygen = a3 (ox) + c (red) given by the equations:
d&l dt = - klZlZz
(1)
dZe/ dt = klZiZa - kaZaZ,
(2)
oxygen Zl Zz + Zs = oxidized + reduced constant reduced acceptor. Z4
(uaJ
complex
=
ELECTRON
TRANSFER
The uninhibited aa cytochrome will become almost immediately oxidized and further changes of reaction rates will be negligible with respect to changes of the cytochromes discussed in our paper. The choice of acceptor Zn in equation 2 wss made between the cytochrome a (belonging to the cytochrome a-us complex where aa is inhibited by CO) and cytochrome c. All three possible choices: (1) that cytochrome a is the only acceptor; (2) that cytochromes a and c are acceptors; and (3) that cytochrome e is the only acceptor; were fitted for identical sets of experimental data. Alternative 3 consistently gives the best fit-in accordance with the experimental results of Wohlrab on CO inhibited, cytochrome c depleted mitochondria (7). The reason that the complicated kinetics of the reactions on the substrate side of cytochrome cl could be approximated by one chemical and two reactions is that the ubiquinone shows a monophasic kinetics under the experimental conditions used in our work. This means (8) that ubiquinone is not rapidly equilibrating with cytochromes a, c and cl . This kinetic conclusion is in agreement with thermodynamic data which show that in uncoupled mitochondria the equilibrium between cytochromes br and cr strongly favors reduction of cytochrome ci . Hence any terms describing the reverse reaction between ubiquinone and cytochrome cl is negligible and we may approximate the carriers on the substrate side as one chemical species as long as only the kinetics of oytochromes a, c, cl are concerned. The final scheme is given in Fig. 1 of the previous paper (3) and the detailed reactions used are presented in Table I of this work. The computer studies (simulation as well as fitting studies) were carried out using a computer program which directly translates the listed reaction to an integration program based on Euler’s numeric integration method (9, 4). For the sake of completeness, the list of differential equations stemming from the scheme in Column I of Table II of this work is given in Footnote’. 4
AR = 5,.
2’2 = k2x4Xl - k3x3xa
A0 = x1
x’p = klX3 - k*x421 + kaxaxs - ~GW~ + ksx~xa (4)
CR = x1
xf8 = ktxgb
-
ksxsxa
(3)
(5)
- kcxsx~ co = x4 CiR = x!i Cl0
=
X8
POOL R = x1 POOL 0 = Zg
~‘8 = k~x~r - k7xs
03
With conservation conditions: x1 + x2 = c0nst.r ; x3 + x4 = const.~; x6 + 26 = c0nst.r; x7 + ~8 = const. 4, expressing the fact that any cytochrome
REACTIONS
677
The estimation of the initial conditions. The difference between the most reduced level (obtainable in the absence of oxygen and the highest oxidized level (obtainable in the presence of oxygen and two noncompetitive inhibitors, rotenone and TTFB, and one competitive inhibitor, malonate) was taken as the amount of biologically active cytochromes. As long ss some electron flow from the substrates occurs in the respiratory chain, the carriers are never 196% oxidized. If the amount of inhibitor used in the experiment is too small, the differences between the highest oxidized level and loO’% oxidation are excessively large. On the other hand if the inhibition of the reaction between the substrate and the cytochrome chain is complete, mitochondria never return to the reduced, anaerobic form. In the present experimental series a concentration of inhibitors has been used which gives 93-97yo oxidation of cytochromes a, c, cl and consumption of 16 PM 02 in about l-2 min. The choice of the proper initial condition for any model is complicated by the fact that the components of the respiratory chain are membrane bound and their active concentration within the membrane cannot be changed or calculated as if they were free in solution. For this reason we introduce the use of the relative concentration of the components as their dimensionless ratio with respect to cytochrome c (defined as equal to 1). This makes the concentration of the individual cytochromes independent of the dilution factor employed in model studies as is the case experimentally. Since the relative concentrations are dimensionless numbers, the second order velocity constants presented here have a formal dimension of set-1. For direct physical interpretation in absolute units it would be necessary to take account of the quantity of cytochrome c in the mitochondrial membranes. Since, however, the pseudo first-order rate constants are dilution independent as well as concentrations of the cytochromes inside the mitochondria, the validity of second-order rate conctant derived in this way would be doubtful. The optimization procedure and accuracy. For optimization of the parameters, a rapidly convergent program based on the Davidon (9), Fletcher and Powell (10) method was employed
might exist only in the reduced or oxidized state, and the sum of both states should be constant. The conservation conditions permit us to eliminate the differential equations for 2’1 , 5’3 ) 2’3 , 2’8 . The constants (l-4) correspond to the total analytical concentrations of the cytochromes.
675
WAGNER,
ERECIRSKA, TABLE
AND PRING II
LINT OF REACTIONS AND VALUES OF RATE CONSTANTS FOR THE FITS WITH Low AND Hrc+n AMOUNTS OF ANTIMYCIN Aa Only constants 6, 7 Rate All constantsfitted High and concentration of All constants fitted LOW Reaction ,E”,; antimycin A experiment pool fitted Low anti- antimycin A experiment my& A experiment 1
Cred
+
=
cox
2 Co, + Ared = Aox + 3 Cred + A,, = Co, + 4 co, + Clred = clox + 5 cred + clox = Clred + 6 &ax + Prsd = POX + 7 Pox + = Pred
1 2 3 Cred 4 cOx 5 Chd 6 7 Pool concentration Cred Ared
5.8 220 65 28 20 0.17 0.020 6.3 S&R 0.08142 points SQR/l point 0.0022
2.0 0.19 6.9 0.20181 points 0.0625
5.1 237 65 31 18 1.7 0.20 6.2 0.12/81 points 0.0015
a Relative concentration of cytochromes: a = 0.44; c = 1.0; c, = 0.52. with the initial estimation of the diagonal members of the Hessian matrix equal to the square of the corresponding constants. For the variable “POOL” the corresponding diagonal member wa8 taken aa 1. The normalized sum of the squares of deviations (SQR) between the experimental and computer evaluated concentrations was used as the function of interest. The Hessian matrix (11-12) which is evaluated during the fitting process in Davidon’s type of program, gives a picture of the accuracy with which the parameters are estimated with respect to the given set of experimental data and within the given model. This is possible, since at the point where the values of constants give the minimum value of the function of interest, the Hessian matrix should correspond to the variance-covariance matrix (8). Due to the rounding errors and the fact that S&R is not a quadratic form of the function of parameters, the Hessian matrix is giving, in our case, only first approximation of this error. However, for biochemical experiment this approximation is usually sufficiently sophisticated. The inaccuracy with respect to the parameter usually becomes greater with an increase in the number of parameters, whereas, at the same time, S&R usually decreases. The accuracy of the determination of parameters at the point of best fit may also be estimated by the following criteria (A and B). (A). S&R is a function of the parameters and the optimization process finds, by a fitting process, the values of the parameters for which S&R is the smallest. If we consider S&R to be a one-dimensional function of one of the parameters, with other parameters as constants, we get the functions shown in Fig. 1 for constants kr and kt of our model. Under these circumstances, it is clear that a better determined
3 FIG. 1. The function of interest aa a one-dimensional function of a chosen velocity constant. Ordinate: 10% change of the value of constant; Abscissa: value of S&R X 10. S&R was evaluated for experimental points of cytochrome c, ci and a from an experiment with high antimycin A block. Curve 1 is for constant 1 and curve 2 for constant 2 of Table III.
parameter is that which has the sharper “valley” and hence for which the corresponding diagonal member in the Hessian matrix is smaller. This is the independent control of the accuracy of estimation of the constant from Hessian matrices. (B) The other criterion of accuracy comes from the Lipschitz condition (11) which can be presented without complicated mathematics aa follows : S&R might be already very low and the diagonal members of the Hessian matrix small, but this does not necessarily tell us whether we have already found the minimum with respect to all the constants, i.e., that we are on the bottom of the
ELECTRON
TRANSFER
“valley” for all parameters. Further decreases in some direction are still possible as long aa any of the nondiagonal members of the second derivative matrix (and the corresponding member of the Hessian matrix) are negative. There exists a possibility that the minimum found is not the only local minimum existing over the range of all the values of the constants. There are two different cases for which more than one minimum might be found. The first case is when the amount of experimentally available information is not sufficient for the model used. This situation is mentioned together with numerical examples in the discussion. The second situation occurs when the initial estimation of the constants or of the initial conditions is too inaccurate. In such a case a false minimum is sometimes found by the computer. The S&R in this kind of minimum is usually much higher than in the “true” one and the graphical representation of a solution usually indicates that some important features (the presence of the inflex point, or tbe maxima, or the number of steady states, etc.) on the theoretical and experimental curves do not correspond to each other. RESULTS
The experimental curves for the kinetics of oxidation of cytochromes a, c and cl obtained in the presence of two antimycin A concentrations (see methods) are shown in Fig. 2. The curves exhibit characteristics typical of antimycin A and CO blocked systems described in paper I of this series (3). The initial fit is first obtained for the
experimental curves with a high amount of antimycin A (Fig. 3). Estimations of the accuracy with which the parameters are determined are shown in Table I. The concentration variable (“POOL”) and constant k7 can be chosen only approximately and are not estimated properly with respect to these sets of experimental data. This is shown by the fact that some of the non-diagonal members of the Hessian matrix for these two parameters are still negative. Except for constant k, and the “POOL” concentration, for which the Lipschitz condition is not completely obeyed for the experiment with high antimycin A level, all the other constants are determined with sticient accuracy to be relevant to the model employed. After the initial fit is completed and the constants obtained, the values for all these constants except B are transferred as the initial approximation to the experiment with a low antimycin A block. Velocity constants which have been estimated with relative error determined from Hessian matrices larger than 200% (,+,, size of pool) and constant i& are fitted first; when no further improvement of the fit can be obtained, all constants are fitted Cytochrome 0
Flow Velocity Trace
Starts stops Ir 4iz-L
679
REACTIONS
c
100%
ox.
Qn cyio. c 1
s/afitops
I f
D 0
set
0% ox.
3
550-540nm Absorbance Increase 1 4
klsec
f 17fiM 0,
A
f 17pM 0,
+
tdsec
LB
FIG. 2. The effect of low and high antimycin
A block on the oxidation traces of cytochrome c in CO blocked pigeon heart mitoohondria. Medium: 0.3 M mannitol-sucrose, 40 mu MOPS pH 6.9, 6.0 mM succinate, 1.5 mM glutamate, 5.0 pCcM rotenone, 2.0 mM malonate, 30 /JM CO, 0.06 pg antimycin A/mg prot (A) and 0.12 pg antimycin A/mg prot (B) 2.7 mg prot/ml.
set
0%
ox.
3
FIG. 3. Fit of calculated oxidation levels (solid lines) to experimental data ( q I) for the high antimycin A block. The values of constants are given in column 2 of Table II.
680
WAGNER,
ERECIRSKA,
together (Fig. 4). The values of velocity constants for the fits in both antimycin A concentrations are summarized in Table II. With the basis for initial estimations taken from the high antimycin A block experiment the constants obtained for the low inhibitor block are approximately x to x better estimated than those from experiments using a complete antimycin A block. Cytochrome
c
6
0
0
6
set
0% ox
0% ox.
FIG. 4. Fit of calculated oxidation levels (solid lines) to experimental data (0) for the low antimycin A block. The values of constants are given in column 3 of Table II.
COMPUTER OBTAINED FITS
FOR
Reactions
1 3 4 5 6 7 8
SQR
c,, + Cred + C 10x + C bed +
DISCUSSION
The numerical results of the fitting procedure for cytochrome kinetics show that the model presented in the first part of this work is in quantitative agreement with the experimental data and can be used BS the basis for future studies of the kinetics of the cytochrome system. In various fields of modern science where the possibility exists of comparing the experimental results with mathematical models curve fitting computer studies become of increasing importance (12-14). The method however has never enjoyed confidence among biochemists or biologists since there exists a certain scientific justification for caution towards curve fitting. In order to illustrate the dangers of the method we would like to present an example from which it can be easily deduced how important is the question of the relevance of the experimental data with respect to the model used. The example is taken from an older fitting attempt where the primary input of electrons to unblocked cytochrome a3 was assigned to cytochrome a instead of c. The experiment was carried out at 17”, otherwise the conditions were identical to those presented in this paper. Two independent fittings were carried out from the same experimental data, initial conditions and initial approximation of the constants. The only difference was in the use of different initial estimation of the Hessian matrix.
TABLE III ELECTRON TRANSFER REACTIONS IN ANTIMYCIN PIGEON HEART
Number
AND PRING
Ctred --f Cred + cl,, &XC -+ Clred + CO, Ctred Gor
A AND CO BLOCKED
MITOCHONDRIA~ Fit
lb
0.75 0.22 1.03 2.30 10-Z 4.2 10.2 0.2 0.010/20 points
Fit
2’
0.72 0.20 0.74 1.60 0.27 2.10 1.50 0.003 0.011/20 points
o Used for two different initial estimations of the Hessian matrix. Cytochromes c and a were obtained experimentally. Temperature 17’. 60 PM CO, 0.2 fig antimycin A/mg mitochondrial protein. b Fit 1 with initial estimation of diagonal members of the Hessian matrix equal 1 (identity matrix). c Fit 2 with initial estimation of diagonal members of the Hessian matrix equal to squares of the corresponding constants.
ELECTRON
TRANSFER
The results for both fits are summarized in Table III. The sum of squares is, in both cases, very low. However the constants of the reactions of cytochrome cl, which was not measured in this experiment, show great discrepancies for both fits. The explanation is that the number of experimentally unavailable substances involved in the model must be limited to those which make mutually distinguishable effects upon the curves of experimentally available substances. In the example given here the total effect of all substances on the substrate side of cytochrome c are indistinguishable by the fitting algorithm. The set of experimental data for cytochrome a and c alone did not sufficiently determine the value for all the constants in the model used: hence, the experimental data were not completely adequate with respect to the chosen model (or vice versa). In the fitting case presented in this work the experimental data are relevant to the model. The region of the conditions for which quantitative agreement between the data and the model can be achieved is necessarily n.arrower than that for a qualitative one. Cytochrome us has to be suflicient’ly blocked by CO (95-99 %) to enable the approximation of the unblocked part of the chain by one reaction and to provide a suitable density of points on the experimental traces to ensure sufficient accuracy for the fitting procedure. The amount of information required to obtain the velocity constants in the quantitative study is seemingly higher than in the qualitative one. To achieve this goal the methods of linked comparison of parameters were used: i.e., carrying the parameters from a less complicated system (representing a svstem kinetically isolated by means of inhibitors from a ‘kinetic neighborhood”) to a more complicated model. For an adequate model only the parameters which were changed in the experimental system are allowed to change on fitting of the data from experiments without inhibitors. The results of the fitting can be briefly summarized as follows: in a high antimycin A block the level of the steadystate is determined by the values of constants k1 and lcs and the differences in slopes of the curves of the individual cyto-
REACTIONS
681
chromes by the “inner constants” kz - kS . In the low antimycin A block the jump to the first steady-state is given by the combined effects of constants kl and lea and the size of the “POOL”. The relations between constants kl and k, determine the second steady state. The ‘Smer constants” are determined from the differences in the slopes of the curves of the individual cytochromes. It is very interesting that the ratio of the inner constants-the equilibrium constant is given with greater accuracy than the individual constants themselves. This is also confirmed from the value of corresponding nondiagonal terms in the Hessian matrix which can in this case be considered as a variancecovariance matrix. The possibility of determining the equilibrium constants by means of this study is of considerable importance for the confrontation of the kinetic and thermodynamic data, for which only the ratio of the constants is used. The results of the fitting show no observable discrepancy between the experimental data and the model nor between the values of constants determined from curves obtained under different antimycin A levels. This confirms the suggestion that the binding of cytochromes c, cl and a in the mitochondrial membrane is weak enough to permit different kinds of mutual interaction which could be described by models of our type. The fitting procedure finds as well a suitable mutual stoichiometric relation between the three cytochromes studied here. No specific effect of rigid structural organization is necessary to explain the kinetic behavior of those cytochromes in the membranes treated with different inhibitors. ACKNOWLEDGMENTS The authors wish to express their gratitude to Dr. Joseph Higgins and Dr. David Garfinkel for their helpful discussion. REFERENCES 1. EXNER, 0. Collect. Czech. Chem. Commun. 35, 187. (1970). 2. HEMPEL, C. G. “Minnesota
Studies in the Philosophy of Science”, Vol. 3, 98, University of Minnesota Press, 1966. 3. WAGNER, M. AND ERECINSKA, M. Arch. Biothem. Biophys. 147,666 (1971).
652
WAGNER,
ERECINSKA,
4. CHANCE, B. AND PRING, M. Coil. Gesell. Biolog. Chemie., MosbachlBaden, p. 102, (1969). Springer Verlag, Berlin/Heidelberg. 5. PRING, M. “Concepts and Models in Biomathematics,” (F. Heinmets, M. Dekker, Eds.), Chaper 2, pp. 75-103, New York, (1970). 6. CHANCE, B., ERECINSKA, M., AND WAGNER, M. Ann. N.Y. Acad. Sci. 174,193 (1970). 7. WOHLRAB, H. Biochemistry 9,474 (1970). 8. J. HIGGINS. Thesis, University of Pennsylvania. “A theoretical study of the kinetic properties of sequential enzyme reactions” 1959.
AND
PRING
9. DAVIDON, W. C. A.E.C. Res. Development Rep., ANL-5990 (rev), (1959). 10. FLETCHER, R., AND POXJELL, M. J. D. Computer J. 6, 163 (1963). 11. HENRICI, A. “Numeric Analysis,” Chapter 2, 3, Elsevier, New York, (1968). 12. KOI~ALIK, J., AND OSBORNE, M. R. “Methods for Unconstrained Optimization Problem,” Chapter 3.6-3.7, p. 4148, Elsevier, New York, (1968). 13. Swann, W. H. FEBSLetters 2,suppl. 539 (1969). 14. FLETCHER, R., AND POWELL, M. J. D. “Constrained Optimization,” Pergamon Press, (1970).
OF BIOCHEMISTRY
Theoretical
and
AND
BIOPHYSICS
147, 675-682 (1971)
Experimental
Reactions
Studies
of the Electron
in the Cytochrome
II. Fitting of the Experimental MICHAEL Johnson
WAGNER,
MARIA
ERECI@XA
Transfer
Chain’ Data AND
MARTIN
PRING
University of Pennsylvania, Philadelphia, Pennsylvania 19104 Received August 30, 1971;accepted September 17, 1971
Research Fomdation,
The experimental curves of cytochromes a, c, and cl were fitted to the kinetic model presented in the first part of this study. Quantitative agreement between the experiments and the model has been found. Estimation of the initial conditions for computer treatment of the cytochrome kinetics, the applicability of the nonlinear optimization program to the biological sciences and the determination of the errors with respect to the constants are presented.
In contrast to modern physics, there are still no explicit means in the mathematical language of chemical kinetics to relate the observer and the measuring process to the analysis of the results. Thus the only way to establish the validity of a mode, is to compare the experimental curves with those derived theoretically. But even then it is diacult to say how far the agreement between theory and experiments confirms the model. However, a qualitative treatment based on accurately estimated parameters can lead to certain predictions in the experimental field. This is a positive manner in which mathematical tools can be utilized in modern biology. The need for continuous critical comparison of the simulated models and experiments in the process of searching for a reaction mechanism is also suggested by the absence of adequate mathematical language for the evaluation of biochemical 1 Supported by USPH GM Computer time was supported
12202 and 16501. by NIH Grant
#RR-15. All the programs used in this work are available upon request. The Davidon optimization program was modified for biochemical purposes, and the plotting programs are written in FORTRAN IV; the integrating program is in PDP-B-10 MACRO assembly language.
kinetics. There is no commonly valid criterion leading to the acceptance of one model as “true” or “convenient”2 and the rejection of another as “false” or “inconThe criteria of mathematical venient”. statistics and testing of hypotheses might sometimes be applied (l), but the most important criterion is the “predictive of the model for other indeproperty” pendent experiments, independent at least in the mathematical sense if not in practical organization of the experiments. In the first part of this study (3) we have shown qualitative agreement between the effect of inhibitors of the cytochrome system and a kinetic model able to simulate most of the important responses. We will present quantitative results obtained by using a nonlinear optimization technique to determine velocity constants and compare more closely the model presented in the first part of this study with the experimental behavior of the cytochrome system. This work is the result of fitting experimental curves of cytochrome oxidation in the presence of two different concentra2 The term convenient is normally used in more abstract sciences (a), but is equally justified in biochemistry.
675
WAGNER, TABLE DETERMINATION PARSMETERS EXPERIMENT BLOCK'"*~
ERECIP;TSKA,
I
OF THE ACCURACY OF FOR THE FINAL FIT FROM WITH HIGH ANTIMYCIN
THE THE
A
Velocity constant Nbv-
Value
1 2 3 4 5 6 7
5.8 221 65 28 20 0.18 0.02
E
0.67 1.51 0.99 1.54 1.14 2.53 4.97b
A(+40%)
0.11 0.0088 0.0097 0.016
-4
0.023 x 10-b
A(-40%)
0.49 0.026 0.016 0.077 0.017 0.015 -4 x 10-s
a In column 3 the accuracy of the parameters is estimated from the values of the diagonal members of the Hessian matrices. In columns 4 and 5 the accuracy of the parameters is estimated from the changes of S&R if the constants are changed by 40’%. A =
(S&R&
+
40%)-SQR(kj))
The dimensionless numbers correspond to the E = Hi,i/ki and are interpreted as relative error. b Some nondiagonal members of matrix still negative. ’ Reaction scheme and S&R identical with those in table II.
tions of antimycin A. It represents the situation in which the numerical values of all constants except that of the reaction No. 5 of Table I are predicted by one experiment (in the presence of 0.12 pg antimycin A/mg prot. A) and the accuracy of quantitative of prediction for the second experimental (in the presence of 0.06 kg of antimycin A/mg prot.) is tested. The results of previous fittings of the experimental curves of cytochrome kinetics (H), lead to the detection of several inaccuracies in the assumptions and experimental approach to the study of oxidation in the cytochrome system. Therefore, an attempt has been made to obtain a fit of the experimental data with high accuracy with respect to the sum of squares of deviations as well as to the parameters.a The results presented here are intended to a The term parameters means unknown variables, experimentally unavailable. In our case. these are the values of the velocity constants (which are never determined by direct experiment) and the concentration of the substance “POOL”.
AND
PRING
increase the accuracy of previous fits. Fitting the experimental data by means of our optimization programs led us to an understanding of the limits of the present model. In addition, optimization produced a better estimate of the velocity constants to be used for simulation studies with extended models. METHODS Experimental. The experimental data were obtained in a stopped flow apparatus using pigeon heart mitochondria, as described in the first paper of this series. A low concentration of antimycin A was added to the system inhibited with rotenone, TTFB and CO. The amount of antimycin A was doubled and the traces were retaken. Since the fitting process is much more sensitive to experimental inaccuracies than a simulation study, precautions were taken to ensure the necessary high reproducibility of the experimental data. Oxygen was delivered to the anaerobic mitochondrial mixture at nearly identical intervals (5 min) and the experimental traces were taken in duplicate. The choice of model. The most difficult and vulnerable part of any study which is neither purely experimental nor purely theoretical is the relation between the experiments and the model used to interpret the results, i.e., t.he choice of some model and the translation of one model to another. For this reason, and because the individual steps in our reaction scheme might seem arbitrary without a more extensive discussion, we describe here in some detail the design of a model used in this work. The reactions considered are those concerning the experimentally measured cytochromes a, c, cl . An attempt was made to approximate the part of the reaction scheme concerning the cytochromes not discussed in this work by as few reaction steps and variables aa possible. The system of reactions transferring electrons from that part of the cytochrome system studied here to oxygen (via CO uninhibited cytochromes a-a3 can be approximated under our experimental conditions as one first order reaction. This is possible because in differential equations for all reactions between a3 (red) + oxygen = a3 (ox) + c (red) given by the equations:
d&l dt = - klZlZz
(1)
dZe/ dt = klZiZa - kaZaZ,
(2)
oxygen Zl Zz + Zs = oxidized + reduced constant reduced acceptor. Z4
(uaJ
complex
=
ELECTRON
TRANSFER
The uninhibited aa cytochrome will become almost immediately oxidized and further changes of reaction rates will be negligible with respect to changes of the cytochromes discussed in our paper. The choice of acceptor Zn in equation 2 wss made between the cytochrome a (belonging to the cytochrome a-us complex where aa is inhibited by CO) and cytochrome c. All three possible choices: (1) that cytochrome a is the only acceptor; (2) that cytochromes a and c are acceptors; and (3) that cytochrome e is the only acceptor; were fitted for identical sets of experimental data. Alternative 3 consistently gives the best fit-in accordance with the experimental results of Wohlrab on CO inhibited, cytochrome c depleted mitochondria (7). The reason that the complicated kinetics of the reactions on the substrate side of cytochrome cl could be approximated by one chemical and two reactions is that the ubiquinone shows a monophasic kinetics under the experimental conditions used in our work. This means (8) that ubiquinone is not rapidly equilibrating with cytochromes a, c and cl . This kinetic conclusion is in agreement with thermodynamic data which show that in uncoupled mitochondria the equilibrium between cytochromes br and cr strongly favors reduction of cytochrome ci . Hence any terms describing the reverse reaction between ubiquinone and cytochrome cl is negligible and we may approximate the carriers on the substrate side as one chemical species as long as only the kinetics of oytochromes a, c, cl are concerned. The final scheme is given in Fig. 1 of the previous paper (3) and the detailed reactions used are presented in Table I of this work. The computer studies (simulation as well as fitting studies) were carried out using a computer program which directly translates the listed reaction to an integration program based on Euler’s numeric integration method (9, 4). For the sake of completeness, the list of differential equations stemming from the scheme in Column I of Table II of this work is given in Footnote’. 4
AR = 5,.
2’2 = k2x4Xl - k3x3xa
A0 = x1
x’p = klX3 - k*x421 + kaxaxs - ~GW~ + ksx~xa (4)
CR = x1
xf8 = ktxgb
-
ksxsxa
(3)
(5)
- kcxsx~ co = x4 CiR = x!i Cl0
=
X8
POOL R = x1 POOL 0 = Zg
~‘8 = k~x~r - k7xs
03
With conservation conditions: x1 + x2 = c0nst.r ; x3 + x4 = const.~; x6 + 26 = c0nst.r; x7 + ~8 = const. 4, expressing the fact that any cytochrome
REACTIONS
677
The estimation of the initial conditions. The difference between the most reduced level (obtainable in the absence of oxygen and the highest oxidized level (obtainable in the presence of oxygen and two noncompetitive inhibitors, rotenone and TTFB, and one competitive inhibitor, malonate) was taken as the amount of biologically active cytochromes. As long ss some electron flow from the substrates occurs in the respiratory chain, the carriers are never 196% oxidized. If the amount of inhibitor used in the experiment is too small, the differences between the highest oxidized level and loO’% oxidation are excessively large. On the other hand if the inhibition of the reaction between the substrate and the cytochrome chain is complete, mitochondria never return to the reduced, anaerobic form. In the present experimental series a concentration of inhibitors has been used which gives 93-97yo oxidation of cytochromes a, c, cl and consumption of 16 PM 02 in about l-2 min. The choice of the proper initial condition for any model is complicated by the fact that the components of the respiratory chain are membrane bound and their active concentration within the membrane cannot be changed or calculated as if they were free in solution. For this reason we introduce the use of the relative concentration of the components as their dimensionless ratio with respect to cytochrome c (defined as equal to 1). This makes the concentration of the individual cytochromes independent of the dilution factor employed in model studies as is the case experimentally. Since the relative concentrations are dimensionless numbers, the second order velocity constants presented here have a formal dimension of set-1. For direct physical interpretation in absolute units it would be necessary to take account of the quantity of cytochrome c in the mitochondrial membranes. Since, however, the pseudo first-order rate constants are dilution independent as well as concentrations of the cytochromes inside the mitochondria, the validity of second-order rate conctant derived in this way would be doubtful. The optimization procedure and accuracy. For optimization of the parameters, a rapidly convergent program based on the Davidon (9), Fletcher and Powell (10) method was employed
might exist only in the reduced or oxidized state, and the sum of both states should be constant. The conservation conditions permit us to eliminate the differential equations for 2’1 , 5’3 ) 2’3 , 2’8 . The constants (l-4) correspond to the total analytical concentrations of the cytochromes.
675
WAGNER,
ERECIRSKA, TABLE
AND PRING II
LINT OF REACTIONS AND VALUES OF RATE CONSTANTS FOR THE FITS WITH Low AND Hrc+n AMOUNTS OF ANTIMYCIN Aa Only constants 6, 7 Rate All constantsfitted High and concentration of All constants fitted LOW Reaction ,E”,; antimycin A experiment pool fitted Low anti- antimycin A experiment my& A experiment 1
Cred
+
=
cox
2 Co, + Ared = Aox + 3 Cred + A,, = Co, + 4 co, + Clred = clox + 5 cred + clox = Clred + 6 &ax + Prsd = POX + 7 Pox + = Pred
1 2 3 Cred 4 cOx 5 Chd 6 7 Pool concentration Cred Ared
5.8 220 65 28 20 0.17 0.020 6.3 S&R 0.08142 points SQR/l point 0.0022
2.0 0.19 6.9 0.20181 points 0.0625
5.1 237 65 31 18 1.7 0.20 6.2 0.12/81 points 0.0015
a Relative concentration of cytochromes: a = 0.44; c = 1.0; c, = 0.52. with the initial estimation of the diagonal members of the Hessian matrix equal to the square of the corresponding constants. For the variable “POOL” the corresponding diagonal member wa8 taken aa 1. The normalized sum of the squares of deviations (SQR) between the experimental and computer evaluated concentrations was used as the function of interest. The Hessian matrix (11-12) which is evaluated during the fitting process in Davidon’s type of program, gives a picture of the accuracy with which the parameters are estimated with respect to the given set of experimental data and within the given model. This is possible, since at the point where the values of constants give the minimum value of the function of interest, the Hessian matrix should correspond to the variance-covariance matrix (8). Due to the rounding errors and the fact that S&R is not a quadratic form of the function of parameters, the Hessian matrix is giving, in our case, only first approximation of this error. However, for biochemical experiment this approximation is usually sufficiently sophisticated. The inaccuracy with respect to the parameter usually becomes greater with an increase in the number of parameters, whereas, at the same time, S&R usually decreases. The accuracy of the determination of parameters at the point of best fit may also be estimated by the following criteria (A and B). (A). S&R is a function of the parameters and the optimization process finds, by a fitting process, the values of the parameters for which S&R is the smallest. If we consider S&R to be a one-dimensional function of one of the parameters, with other parameters as constants, we get the functions shown in Fig. 1 for constants kr and kt of our model. Under these circumstances, it is clear that a better determined
3 FIG. 1. The function of interest aa a one-dimensional function of a chosen velocity constant. Ordinate: 10% change of the value of constant; Abscissa: value of S&R X 10. S&R was evaluated for experimental points of cytochrome c, ci and a from an experiment with high antimycin A block. Curve 1 is for constant 1 and curve 2 for constant 2 of Table III.
parameter is that which has the sharper “valley” and hence for which the corresponding diagonal member in the Hessian matrix is smaller. This is the independent control of the accuracy of estimation of the constant from Hessian matrices. (B) The other criterion of accuracy comes from the Lipschitz condition (11) which can be presented without complicated mathematics aa follows : S&R might be already very low and the diagonal members of the Hessian matrix small, but this does not necessarily tell us whether we have already found the minimum with respect to all the constants, i.e., that we are on the bottom of the
ELECTRON
TRANSFER
“valley” for all parameters. Further decreases in some direction are still possible as long aa any of the nondiagonal members of the second derivative matrix (and the corresponding member of the Hessian matrix) are negative. There exists a possibility that the minimum found is not the only local minimum existing over the range of all the values of the constants. There are two different cases for which more than one minimum might be found. The first case is when the amount of experimentally available information is not sufficient for the model used. This situation is mentioned together with numerical examples in the discussion. The second situation occurs when the initial estimation of the constants or of the initial conditions is too inaccurate. In such a case a false minimum is sometimes found by the computer. The S&R in this kind of minimum is usually much higher than in the “true” one and the graphical representation of a solution usually indicates that some important features (the presence of the inflex point, or tbe maxima, or the number of steady states, etc.) on the theoretical and experimental curves do not correspond to each other. RESULTS
The experimental curves for the kinetics of oxidation of cytochromes a, c and cl obtained in the presence of two antimycin A concentrations (see methods) are shown in Fig. 2. The curves exhibit characteristics typical of antimycin A and CO blocked systems described in paper I of this series (3). The initial fit is first obtained for the
experimental curves with a high amount of antimycin A (Fig. 3). Estimations of the accuracy with which the parameters are determined are shown in Table I. The concentration variable (“POOL”) and constant k7 can be chosen only approximately and are not estimated properly with respect to these sets of experimental data. This is shown by the fact that some of the non-diagonal members of the Hessian matrix for these two parameters are still negative. Except for constant k, and the “POOL” concentration, for which the Lipschitz condition is not completely obeyed for the experiment with high antimycin A level, all the other constants are determined with sticient accuracy to be relevant to the model employed. After the initial fit is completed and the constants obtained, the values for all these constants except B are transferred as the initial approximation to the experiment with a low antimycin A block. Velocity constants which have been estimated with relative error determined from Hessian matrices larger than 200% (,+,, size of pool) and constant i& are fitted first; when no further improvement of the fit can be obtained, all constants are fitted Cytochrome 0
Flow Velocity Trace
Starts stops Ir 4iz-L
679
REACTIONS
c
100%
ox.
Qn cyio. c 1
s/afitops
I f
D 0
set
0% ox.
3
550-540nm Absorbance Increase 1 4
klsec
f 17fiM 0,
A
f 17pM 0,
+
tdsec
LB
FIG. 2. The effect of low and high antimycin
A block on the oxidation traces of cytochrome c in CO blocked pigeon heart mitoohondria. Medium: 0.3 M mannitol-sucrose, 40 mu MOPS pH 6.9, 6.0 mM succinate, 1.5 mM glutamate, 5.0 pCcM rotenone, 2.0 mM malonate, 30 /JM CO, 0.06 pg antimycin A/mg prot (A) and 0.12 pg antimycin A/mg prot (B) 2.7 mg prot/ml.
set
0%
ox.
3
FIG. 3. Fit of calculated oxidation levels (solid lines) to experimental data ( q I) for the high antimycin A block. The values of constants are given in column 2 of Table II.
680
WAGNER,
ERECIRSKA,
together (Fig. 4). The values of velocity constants for the fits in both antimycin A concentrations are summarized in Table II. With the basis for initial estimations taken from the high antimycin A block experiment the constants obtained for the low inhibitor block are approximately x to x better estimated than those from experiments using a complete antimycin A block. Cytochrome
c
6
0
0
6
set
0% ox
0% ox.
FIG. 4. Fit of calculated oxidation levels (solid lines) to experimental data (0) for the low antimycin A block. The values of constants are given in column 3 of Table II.
COMPUTER OBTAINED FITS
FOR
Reactions
1 3 4 5 6 7 8
SQR
c,, + Cred + C 10x + C bed +
DISCUSSION
The numerical results of the fitting procedure for cytochrome kinetics show that the model presented in the first part of this work is in quantitative agreement with the experimental data and can be used BS the basis for future studies of the kinetics of the cytochrome system. In various fields of modern science where the possibility exists of comparing the experimental results with mathematical models curve fitting computer studies become of increasing importance (12-14). The method however has never enjoyed confidence among biochemists or biologists since there exists a certain scientific justification for caution towards curve fitting. In order to illustrate the dangers of the method we would like to present an example from which it can be easily deduced how important is the question of the relevance of the experimental data with respect to the model used. The example is taken from an older fitting attempt where the primary input of electrons to unblocked cytochrome a3 was assigned to cytochrome a instead of c. The experiment was carried out at 17”, otherwise the conditions were identical to those presented in this paper. Two independent fittings were carried out from the same experimental data, initial conditions and initial approximation of the constants. The only difference was in the use of different initial estimation of the Hessian matrix.
TABLE III ELECTRON TRANSFER REACTIONS IN ANTIMYCIN PIGEON HEART
Number
AND PRING
Ctred --f Cred + cl,, &XC -+ Clred + CO, Ctred Gor
A AND CO BLOCKED
MITOCHONDRIA~ Fit
lb
0.75 0.22 1.03 2.30 10-Z 4.2 10.2 0.2 0.010/20 points
Fit
2’
0.72 0.20 0.74 1.60 0.27 2.10 1.50 0.003 0.011/20 points
o Used for two different initial estimations of the Hessian matrix. Cytochromes c and a were obtained experimentally. Temperature 17’. 60 PM CO, 0.2 fig antimycin A/mg mitochondrial protein. b Fit 1 with initial estimation of diagonal members of the Hessian matrix equal 1 (identity matrix). c Fit 2 with initial estimation of diagonal members of the Hessian matrix equal to squares of the corresponding constants.
ELECTRON
TRANSFER
The results for both fits are summarized in Table III. The sum of squares is, in both cases, very low. However the constants of the reactions of cytochrome cl, which was not measured in this experiment, show great discrepancies for both fits. The explanation is that the number of experimentally unavailable substances involved in the model must be limited to those which make mutually distinguishable effects upon the curves of experimentally available substances. In the example given here the total effect of all substances on the substrate side of cytochrome c are indistinguishable by the fitting algorithm. The set of experimental data for cytochrome a and c alone did not sufficiently determine the value for all the constants in the model used: hence, the experimental data were not completely adequate with respect to the chosen model (or vice versa). In the fitting case presented in this work the experimental data are relevant to the model. The region of the conditions for which quantitative agreement between the data and the model can be achieved is necessarily n.arrower than that for a qualitative one. Cytochrome us has to be suflicient’ly blocked by CO (95-99 %) to enable the approximation of the unblocked part of the chain by one reaction and to provide a suitable density of points on the experimental traces to ensure sufficient accuracy for the fitting procedure. The amount of information required to obtain the velocity constants in the quantitative study is seemingly higher than in the qualitative one. To achieve this goal the methods of linked comparison of parameters were used: i.e., carrying the parameters from a less complicated system (representing a svstem kinetically isolated by means of inhibitors from a ‘kinetic neighborhood”) to a more complicated model. For an adequate model only the parameters which were changed in the experimental system are allowed to change on fitting of the data from experiments without inhibitors. The results of the fitting can be briefly summarized as follows: in a high antimycin A block the level of the steadystate is determined by the values of constants k1 and lcs and the differences in slopes of the curves of the individual cyto-
REACTIONS
681
chromes by the “inner constants” kz - kS . In the low antimycin A block the jump to the first steady-state is given by the combined effects of constants kl and lea and the size of the “POOL”. The relations between constants kl and k, determine the second steady state. The ‘Smer constants” are determined from the differences in the slopes of the curves of the individual cytochromes. It is very interesting that the ratio of the inner constants-the equilibrium constant is given with greater accuracy than the individual constants themselves. This is also confirmed from the value of corresponding nondiagonal terms in the Hessian matrix which can in this case be considered as a variancecovariance matrix. The possibility of determining the equilibrium constants by means of this study is of considerable importance for the confrontation of the kinetic and thermodynamic data, for which only the ratio of the constants is used. The results of the fitting show no observable discrepancy between the experimental data and the model nor between the values of constants determined from curves obtained under different antimycin A levels. This confirms the suggestion that the binding of cytochromes c, cl and a in the mitochondrial membrane is weak enough to permit different kinds of mutual interaction which could be described by models of our type. The fitting procedure finds as well a suitable mutual stoichiometric relation between the three cytochromes studied here. No specific effect of rigid structural organization is necessary to explain the kinetic behavior of those cytochromes in the membranes treated with different inhibitors. ACKNOWLEDGMENTS The authors wish to express their gratitude to Dr. Joseph Higgins and Dr. David Garfinkel for their helpful discussion. REFERENCES 1. EXNER, 0. Collect. Czech. Chem. Commun. 35, 187. (1970). 2. HEMPEL, C. G. “Minnesota
Studies in the Philosophy of Science”, Vol. 3, 98, University of Minnesota Press, 1966. 3. WAGNER, M. AND ERECINSKA, M. Arch. Biothem. Biophys. 147,666 (1971).
652
WAGNER,
ERECINSKA,
4. CHANCE, B. AND PRING, M. Coil. Gesell. Biolog. Chemie., MosbachlBaden, p. 102, (1969). Springer Verlag, Berlin/Heidelberg. 5. PRING, M. “Concepts and Models in Biomathematics,” (F. Heinmets, M. Dekker, Eds.), Chaper 2, pp. 75-103, New York, (1970). 6. CHANCE, B., ERECINSKA, M., AND WAGNER, M. Ann. N.Y. Acad. Sci. 174,193 (1970). 7. WOHLRAB, H. Biochemistry 9,474 (1970). 8. J. HIGGINS. Thesis, University of Pennsylvania. “A theoretical study of the kinetic properties of sequential enzyme reactions” 1959.
AND
PRING
9. DAVIDON, W. C. A.E.C. Res. Development Rep., ANL-5990 (rev), (1959). 10. FLETCHER, R., AND POXJELL, M. J. D. Computer J. 6, 163 (1963). 11. HENRICI, A. “Numeric Analysis,” Chapter 2, 3, Elsevier, New York, (1968). 12. KOI~ALIK, J., AND OSBORNE, M. R. “Methods for Unconstrained Optimization Problem,” Chapter 3.6-3.7, p. 4148, Elsevier, New York, (1968). 13. Swann, W. H. FEBSLetters 2,suppl. 539 (1969). 14. FLETCHER, R., AND POWELL, M. J. D. “Constrained Optimization,” Pergamon Press, (1970).