3842
DELOSF. DETARAND CARLETON E. DETAR
General Computer Techniques for Evaluating the Time-Concentration Relationships Predicted by Reaction Mechanisms, Including Complex Enzyme Mechanisms1.2
by DeLos F. DeTar and Carleton E. DeTar Department of Chemistry and the Institute of Molecular Biophysics, Florida Slate University, Tallahassee, Florida (Receaaed M a y 17, 1966)
General techniques are described for computation of reactant and product concentrations as a function of time for a proposed reaction mechanism of almost any complexity. The mechanisms may include equilibria and steady-state systems.
Introduction The study of reaction mechanisms from the standpoint of kinetics has proved very fruitful, but full use of the technique has been limited since mathematical expressions describing a given mechanism often become intractable with even quite simple mechanisms. Methods of chemical analysis have reached the point that experimental sophistication has outstripped the capabiIities of traditional theoretical approaches. A number of workers have therefore turned to computer techniques for the study of complex system^.^-^ I n most cases these have concerned solutions of the specific problems at hand. A computer program provides tables of numerical values of concentrations at given times. Even when analytical solutions are available, as for example in eq 7 and 12 of ref 3d, these may be so complex as to be useful primarily for preparing similar tabular output. The analytical solutions provide the necessary reference points, but the contrast between the complexity of these hard won expressions and the simplicity of the corresponding iterative computer programs is striking. The introduction of powerful and sophisticated computer techniques into the study of reaction mechanisms provides a wholly new dimension and may be expected to yield a wealth of new developments.
Method of Computation General. The basic computation involves numerical integration of systems of differential equations. The equation for the decrement D of a given reaction step The Journal of Physical Chemistry
is eq 1. The simplest form of integration in eq 1 con-
D
hlxl
. ..
XnAt
(1)
sists of computing reaction decrements for each reaction in turn using the reactant concentrations present a t the beginning of the interval. The concentration of a given compound is decremented by D for each equa(1) This work was supported by the Division of Biology and Medicine, U. S. Atomic Energy Commission, under Contract No. AT(40-1)-2690. (2) We wish to express appreciation to Dr. E. P. hliles and to other members of the Computing Center at Florida State University for generous help and to the computing Center for making available the many hours of computer time needed for developing the programs. (3) Among the many excellent references t o traditional treatments are the following: (a) A. A. Frost and R. G Pearson, “Kinetics and Mechanism,” 2nd ed, John Wiley and Sons, Inc., New York, N. Y., 1961; (b) S. Benson, “The Mathematical Foundation of Chemical Kinetics,” McGraw-Hi11 Book Co., Inc., New York, N. Y., 1960. The mathematical complexities of several relatively simple systems have been elegantly treated, for example, by (e) R. A. Alberty, J . Am. Chem. SOC.,80, 1777 (1958); (d) IT. G. Miller and R. A. Alberty, ibid., 8 0 , 5146 (1958). (4) (a) B. Chance, D. Garfinkel, J. Higgins, and B. Hess, J . Biol. Chem., 235, 2426 (1960), and earlier papers: (b) D. Garfinkel and B. Hess, ibid., 239, 971 (1964); (c) B. Chance, ibid., 235, 2440 (1960). ( 5 ) C. Perrin and F. H. Westheimer, J . Am. Chem. SOC.,85, 2773 (1963). (6) K. Wiberg, “Computer Programming for Chemists,” W. A. Benjamin, Inc., New York, N. Y., 1965, p 168. (7) C. F. Walter and 31. F. Morales, J . Biol. Chem., 239, 1277 (1964). (8) (a) C. J. Collins, J. E. Cristee, and V. F. Raaen, J . A m . Chem. SOC.,83, 4267 (1961); (b) C. J. Collins, Advan. Phys. Org. Chem., 2 , 44 (1964). (9) See also the RLICH series and RJICHSS, simple programs of moderate generality: D. F. DeTar, J . Chem. Educ., in press.
COMPUTER TECHNIQUES FOR EVALUATING THE TIME-CONCENTRATION RELATIONSHIPS
tion in which it appears as a reactant and incremented by D for each equation in which it is a product. Stoichiometry is preserved by applying the incrementing or decrementing to each appearance of a given reactant or product. Hence for the second-order reaction A A + B (with rate constant k , ) the amount of B produced is D, and the amount of A which disappears is 2 0 . The effect is to interpret the rate constant as equal to IC with respect to B and 2k with respect to A. Certain of the more sophisticated forms of integration are more efficient.’” A modification of the trapezoidal method (eq 2) appears to be the best compromise. The more elaborate Runge-Kutta method is lower.^
are defined in eq 1.) I n eq 4 the summation over p includes only reactions in which S, is a product and
+
A,
=
1 . 5 0 , - O.5Di-1
sn =
(5kiAiBi)/&kjAj)
(3)
spective coreactants. A and B are the concentrations of A and B, respectively. Inclusion of further terms is obvious if there are more than two reactants or products. A,, B,, and A, may also be steady-state intermediates (including s,). If the set of equations for a given steady-state intermediate does not contain any steady-state intermediates among the A,, etc., then the concentrations are obtained without iteration. I n the more general case the solutions are only a p proximate. One possible iterative procedure is to keep substituting the values of S, found in the previous iteration. However, this primitive approach is not very effective. The first step toward a workable method is to replace eq 3 by the equivalent equations (4 and 5). (Dvalues
(5)
&,(new) = R,S,(old)
that over r only reactions in which S, is a reactant. If S, has its correct value, then R, = 1, while R, > 1 means that the estimate of S, is too small. Equation 4 permits incorporation of a damping factor (DF) to reduce the tendency toward over compensation and hence oscillation. &(new) = [(DF)*R,
(2)
Steady-State Intermediates. Obviously, the simple incrementing and decrementing approach is not applicable to a reaction involving steady-state intermediates, for a definition of a steady-state intermediate appropriate to an interative calculation is a species whose concentration is less than, perhaps very much less than, the reaction flux, i.e., the decrement D for that reaction. Therefore, the Concentration of a steady-state intermediate is never incremented or decremented but is determined instead by direct calculation. The problem of calculating concentrations of a set of steady-state intermediates admits of no simple general solution Iterative procedures must therefore be used, and two independent approaches have been devised. I n the first, the regular method, the concentration of a given steady-state intermediate, S,, is given by eq 3; the i summation includes all those reactions in which S, is a product, and the j summation includes all those reactions in which S, is a reactant and Aj’s are the re-
3843
+ (for R,
> 1)
(64
1 - D F ] (for R,
< 1)
(6b)
1 - DF]S,(old) &(new) = S,(old)/[DF/R,
+
This method of calculation has worked very well with several different mechanisms. The first round of iterations may run to 20-30 to bring all S, to within 0.295, i.e., until R = 1.000 i 0.002 for each R. Once the “reaction” is underway, however, from three to five iterations usually suffice to touch up the concentrations of steady-state intermediates. I n some mechanisms it is necessary to keep the total concentrations of certain families of steady-state intermediates constant. This corresponds physically to a closed system such as an enzyme system in which one family of steady-state intermediates is represented by all the species containing this enzyme. I n this case the concentration of each species is normalized by a factor equal to the ratio of the apparent total enzyme concentration to the desired total concentration. I n the second method, the perturbation method of obtaining steady-state concentrations, the D, is written in generalized form as eq 7 ; S,I is the first steady-state intermediate in reaction n, SnZ is the second, X,1 is the first other type of reactant ( X = concentration of X), etc. Usually there will be only one or two reactants,
D,
=
ICnS,ISn2 . . . X,lXnz
,
. , At
(7)
but it is necessary to formulate general expressions. It is now assumed that for each steady-state intermediate in each equation there is available a preliminary estimate S’ which is reasonably good, and that it is desired to find the correction term 6S which gives the correct value S, through eq 8. Substitution of (8) (IO) (a) E. Whittaker and G. Robinson, “The Calculus of Observations,” 4th ed, Blackie & Son Ltd., London, 1944, p 132; (b) H. Margenau and G. M. Murphy, “The Mathematic of Physics and Chemistry,” D. Van Nostrand Co., Inc., New York, N. Y., 1943, p 450.
Volume 70, Number 18 December 1966
3844
DELOSF. DETARAND CARLETON E. DETAR
into (7) gives (9) where gn is the number of steadyEn1
+Snit
4- 6Sn1
(8)
(9) state intermediates in react,ion n. All terms involving products of two or more 6s have been discarded as representing high-order corrections. If the method converges, any error so introduced is corrected by iteration. The general steady-state relationship (cf. eq 4 with Rn = 1) is eq 10. All summations are taken over only those D terms in which radical S, is a product (summation over p) or a reactant (summation over r). SubstiP
r
EDi - EDj
= 0
tution of eq 9 into eq 10 gives eq 11. It may be noted that eq 11 is linear in the SS terms and that with one
(11) equation for each intermediate it is possible to solve for these corrections. To treat the case of closed systems, the first equation for each family which comprises a closed system is arbitrarily replaced by the constraining equation 12 in which the total concentration is specified as the sum of the concentration of the component species Vi,. The perturbation form of eq 12 is eq 13. i=l
Which version, the regular or the perturbation, is better must be ascertained by trial. I n some merhanisms one is clearly superior; occasionally both must be used in a cyclic order. The direct method is effective for many very complex sets of steady-state intermediates but does not converge well if both numerator and denominator contain a functionally related dominant term. Such an example is a polymerization mechanism with respect to the dominant polymerization step: RA!!. M + RRI.. The perturbation method requires solution of an N x N determinant where N is the number of steady-state species. This rapidly loses efficiency as N becomes large. Acid-Base Equilibria. l 1 , I 2 The equilibria referred to in this section are those that determine the concen-
+
The Journal of Physical Chemistry
tration of hydrogen ion or other lyonium ion or conversely of hydroxide or other lyate ion. For n acid-base systems with mi species in the i t h system, with qij representing the charge and eij the concentration of the j t h species of the ith system, and cit the family total for the ith system, then eq 14 is the conservation of mass expression and eq 15 defines a term Q which may be called the charge excess. If a buffer system is established by mixing certain volumes of j=1
0.1 M phosphoric acid and 0.1 M sodium hydroxide, then Q is equal to the concentration of sodium ion, for example. Acidity constants for each dissociation step are denoted Kij for dissociation of a proton from 1 species of the ith the j t h species to produce t h e j system (eq 16). The concentration of the j t h species of the ith system is given in terms of the initial species by eq 17, with the proviso that the product function is unity for j = 1. Finally, the concentration of the
+
CijilCH+
=
CzjKij
(16)
j t h species of the ith system is given in terms of eft, the total concentration of the system, by combining eq 17 and 14 to give eq 18.
The computation proceeds as follows. A preliminary very rough estimate for cHt is supplied. The individual species concentrations are computed through eq 18, and an estimate of Q ( = & I ) is evaluated by eq 15. If Q' is smaller than Q, then the cH+estimate was too small. By means of an iterative procedure involving logarithmic steps, successive values of C H t are tested until the change is less than some predetermined amount, e.g., c H + accurate to 0.01%. Equilibrium Reactions. These concern potentially complex systems such as those involving carboxylic acid monomers and dimers and their salts in nondissociating solvents. The general equations are given ~~
~
(11) Cf.L. G. Sill&, "Treatise on Analytical Chemistry," I. M. Kolthoff and P. J. Elving, Ed., The Interscience Encyclopedia Inc., New York, N. Y., 1959, p 277. (12) (a) A. R. Emery, J . Chem. Educ., 42, 131 (1965); (b) A. J. Bard and D. M. King, ibid., 42, 127 (1965); (c) S. L. Cooke, Jr., ibid., 42, 620 (1965).
COMPUTER TECHNIQUES FOR EVALUATING THE TIME-CONCENTRATION RELATIONSHIPS
in eq 19, in which the Ri,’s are reactants and the Pij’s are products of the ith reaction with equilibrium constant K , , and the equilibrium expression eq 20. Each of the Rij and Pij may be assigned to a specific family
Rz1 4-Riz
+ . . + Rtm,
-+ ki
A
Pi1 m,
PlZ
+ . . . + PI,,
K,IIR
