Computer-based derivation, reformulation, and simulation of enzyme reaction models

MATHEMATICAL

507

BIOSCIENCES

Computer-Based

Derivation,

Reformulation,

and Simulation

of

Enzyme Reaction Models DONALD

D. FISHER

Department

of Computing

Oklahoma

State

Stillwater,

Oklahoma

and Information

Sciences

University

AND

ARTHUR

R. SCHULZ

Department

of Biochemistry,

Indianapolis,

Indiana

University

Indiana

Communicated

by Richard Bellman

ABSTRACT A computer-based algorithm to derive the rate equation of an enzyme-catalyzed reaction under the assumption of steady state is the basis for additional algorithms that provide the enzyme kineticist with powerful tools to investigate models of the reaction. These algorithms include one to define kinetic constants associated with the mechanism and one to reformulate the rate equation from coefficient form to kinetic form. The reformulation algorithm is a nonnumeric, combinatorial, algebraic expander. A model of the reaction need not be unique. By using a rational approximation to fit experimental data, one can obtain uniqueness of a model to within an equivalence class of models defined in terms of the rational approximation. The ability to derive and reformulate the rate equation on a computer allows the kineticist to propose several models from which he may select those to study further through experimental means. After experimental data are obtained and the model is fitted by a rational approximation, further refinements may be made by iterating the foregoing procedure. The final model is unique to within the equivalence class defined by the rational approximation and order of accuracy prescribed.

INTRODUCTION

the

Earlier [1 , 21 we reported on a new computer-based method for derivation of the rate equation of an enzyme-catalyzed reaction Mathematical

Copyright

Biosciences

6 (1970), 507-517

@ 1970 by American Elsevier Publishing Company,

Inc.

508 mechanism

DONALD

D. FISHER

that is based on the description

AND

ARTHUR

of the mechanism

R. SCHULZ

in terms

of a connection matrix. This method allows the enzyme kineticist to investigate mechanisms that are unfeasible to consider by noncomputer techniques. Furthermore, it obtains all terms of the rate equation in coefficient form systematically, thus eliminating the possibility of mistakes inherent in manual methods [3]. It also has made obsolete the classification of mechanisms by rate equations. Finally, our method for the derivation of rate equations leads directly to other computer-based algorithms, which translate the rate equation from coefficient form to kinetic form and which simulate the reaction based on kinetic parameters that are determined by initial rate studies. In this article an algorithm is presented for the reformulation by a digital computer of the rate equation in kinetic form given that equation in coefficient form. At present the algorithm is set up for batch processing; it is expected, however, that a proposed interactive edition of the algorithm will be better suited for this type of combinatorial problem in which the enzyme kineticist guides the computer to a reformulation of the rate equation in terms of kinetic constants that may be determined experimentally. The ability to derive rapidly the rate equation in kinetic form allows the kineticist to propose and investigate several different models of a reaction. Based on the occurrence of key terms in the denominator of the rate equation, he then can devise experiments to check the models through initial rate, product inhibition, or isotopic exchange studies. Experimental data thus obtained are used as input to computer simulation programs to further verify, alter, and refine the models. It is this type of augmentation of human intellect that shows great promise in the manmachine interactive environment. The algorithm, even though explained in terms of the simple mechanism of Fig. I, is general and applies to a broad spectrum of mechanisms. It is a nonnumeric, combinatorial, algebraic expander. The essential steps in the procedure are 1. derivation of the rate equation in coefficient form; 2. generation of kinetic constant definitions: 3. designation of preferred inhibition constant definitions; 4. designation of undesirable inhibition constant definitions; 5. reformulation of the rate equation in terms of kinetic constants using preferred inhibition constant definitions; 6. interactive iteration of step 5 to arrive at an “optimum” kinetic formulation of the rate equation. Mathenzatical Biosciences 6 (1970), 507-517

COMPUTER

EQ

509

DERIVATIVE

EP

ER

FIG. 1. A model of monoamine

oxidase.

The algorithms we have developed for steps 1, 2, 5, and 6 apply to both manual and computer manipulations ; however, their effectiveness becomes increasingly apparent as one progresses from simple mechanisms to more complex ones. For example, one nine-node mechanism of interest contains over 4000 terms in the denominator and each of four inhibition constants may be defined in more than 90 different ways, a rather staggering combinatorial problem if one were to consider all possibilities. The rate equation reformulation usually is not unique; hence, the description of the reaction in terms of a model need not be unique. For example, the four models depicted in Fig. 2 give rise to the same stoichiometry and agree with experimental findings to the same order of precision for the reaction (thyroid monoamine oxidase) described by A+B+P+Q+R

(models

A+B&P+Q+R

(model

I, II, III) IV).

These models may be said to belong to the same equivalence class as defined by the stoichiometry and experimental findings. A desirable property is one that defines uniqueness of a model for a reaction with respect to an equivalence class. Such a possibility provides the enzyme kineticist with an additional powerful tool in his investigation of a reaction. He must be able to define the equivalence classes, however, if he is to make use of their implications. Enzyme kineticists have used the Lineweaver-Burk plot as a method for defining equivalence classes of reaction models [4]. The slope and intercept of this linear (first-degree) fit to the data provide the kineticist Mathematical 32

Biosciences

6 (1970), 507-517

DONALD

D. FISHER

AND

ARTHUR

R. SCHULZ

EP

1

(P)

FQ

FQB +

(B) MODEL I

EQ

FQ

FQ

EP

MODEL III

MODEL IV

FIG. 2. Models from the same equivalence class.

with two parameters to approximate in an attempt to obtain a reasonable model of the reaction. For the model depicted in Fig. 1 the linear fit may be written

in the form v = cg + c,(A).

All similar models of the reaction for which the linear approximation fits the data to a specified degree of accuracy are defined to belong to the same equivalence class. The approximation procedure may be generalized from a linear approximation to a rational approximation [5] (quotient of two polynomials), which expands the equivalence class possibilities manyfold. Mathematical

Biosciences

6 (1970), 507-517

COMPUTER

511

DERIVATIVE

For the model depicted

in Fig. 1 a rational v=

fit may be of the form

c0 + c,(A) b, + b,(A) *

Note that if 6, = 0, b, = 1, the rational fit becomes a linear fit; hence, the Lineweaver-Burk plot is a special case of the rational approximation. For the rational case, equivalence classes are defined in terms of the highest exponents that appear in the numerator and in the denominator of the approximation. Thus for a given degree of accuracy, the model of the mechanism is unique to within the equivalence class defined by the approximation. KINETIC

CONSTANT

DEFINITIONS

Variables representing the concentrations of the substrates and products present in the rate equation are arranged in an arbitrary sequence. In our example substrates are represented by A and B; products are represented by P, Q, and R. Output from our rate equation derivation program consists of the rate equation numerator followed by the rate equation denominator. The computer output format for the denominator term (A)(B) is illustrated below. A)B) + 2 0 4 5 1 1 1 A)B) + 2 3 4 0 1 1 1 A@) + 2 3 4 5 0 1 1 Each such term is called a (valid) path vector [ 11. To explain the notation, the rate constant for the reaction path from source node 1 to destination node 2 is written

as ki.

The position

in the path vector

represents

the

source node number, and value in that position represents the destination node number. We define kf, to be 1 for all j. Thus the term A)B) + 2 0 4 5 1 1 1 is the representation

for the more familiar

expression

+k;k;k;kfk;k$4)(B). The procedure for redefining the rate equation in terms of kinetic constants depends on the presence of a term such as (A)(B) but not the rate constants themselves. For this reason the three terms in (,4)(B) are lumped together and are referred to as (A)(B)

x coeff(A)(B). Mathematical

Biosciences 6 (1970), 507-517

512

DONALD

Next we define an exponent which consists

of the exponents

in the term.

For

our example,

D. FISHER

AND ARTHUR

R. SCHULZ

vector associated

with each lumped

of each substrate

and product

let the first position

term

appearing

of the exponent

vector correspond to the exponent of A, the second position to the exponent of B, and so on for P, Q, R. For example, the exponent vector for the term in (A)(B) is (1, l,O, (40). The terms that appear (as a sum) in the rate equation denominator and the corresponding exponent vectors of our example are given in Table I (a right parenthesis is used to separate the reactants). TABLE I Term

Exponent vector

4 B) PI 4B) AY-‘) A)Q) BY) B)Q) B)R)

pm P)Q) PW) Q)R) AP)Q) 4P)Q) B)Q)R) f’)P)Q) P)Q)Q) P)Q)R)

For our purposes the kinetic constants fall into three categories, namely, maximal velocities, Michaelis, and inhibition. The equilibrium constant, a thermodynamic constant, may be expressed in terms of the kinetic constants by means of the Haldane relationship. Inhibition constants may have multiple definitions, whereas the other constants have unique definitions. Mathematical Biosciences 6 (1970), 507-517

513

COMPUTER DERIVATIVE

Let P, denote

products,

Ai denote

substrates,

and let the equilibrium

and kinetic constants be denoted as Keq, equilibrium constant ; V,,

V,, maximal

velocities

in forward

and

reverse

directions,

respectively; K,,, Kai, Michaelis constants for product Pi and substrate Ai, respectively; Hpi, Hai, inhibition constants for product Pi and substrate Ai, respectively. Let tcr, . . . , tck be the exponents of the product terms Pi that appear in the numerator

of the rate equation,

that is, terms of the form

Pi’ . * * Kk. Let fir, . . . , p1 be the exponents of the product terms that appear denominator of the rate equation, that is, terms of the form

in the

Similarly, let yr, . . . , ym and 6r, . . . , S, be exponents of substrate terms Ai that appear in the numerator and denominator, respectively, of the rate equation. Let GCbe the maximum value for which GC~ = . - * = ak for some term in the numerator. Define /!?, y, 6 as analogous maximum values of terms with equal exponents. For example, if the numerator contained two product terms of the form PfPi and PtPt, the value of M would be 2. The equilibrium constant is defined as K erl The maximum

S coefficient

of (A: . . . AY,) in the numerator

coefficient

of (Pi f * * Pi) in the numerator

velocities

y g r

’

are defined as

coefficient coefficient

of (Pp . . . Pi) in the numerator of (Pf . . * Pf) in the denominator’

and I/ p f

The Michaelis

coefficient coefficient constant

of (A: * * . AY,) in the numerator of (At . . . At) in the denominator

for each product

is defined

Mathematical

from

Biosciences

* denominator 6 (1970), 507-517

514

DONALD

terms of the rate equation K,, 2

coefficient

Kat g

AND

ARTHUR

R. SCHULZ

as of ((~~._l,j,i

coefficient and for each substrate

D. FISHER

P$)Pa-l)

of (JJC1

in the denominator

P$) in the denominator

as

coefficient

of ((~~zl,j,i

coefficient

Aq)AP1) in the denominator

of (ny=r

A’$ in the denominator

product of the where JJ=, Pjs z P! . . . Pf stands for the arithmetic Inhibition constants, HDi for products and Hai for chemical products. substrates, are defined as a ratio of two terms from the denominator of the rate equation in a form not unlike that of the Michaelis constants except that both substrate and product terms may be present and the rules for An inhibition constant for product exponents are somewhat different. P, is defined as any ratio of the form H,( g

coefficient

of (lJY1

coefficient

A;] n:Sl,j+l

of (nj”=l

(Pfj)Pf”-l)

A;’ n)=l

in the den.

Pfj) in the den.

’

that is, corresponding exponents of substrates Aj are the same as are corresponding exponents of all products Pj except for Pi. The exponent of Pi that appears in the numerator of the ratio is one less than the exponent of Pi that appears in the denominator of the ratio. For our example (PI E P, Pz z Q), some definitions of HP, E H, are constant H, 2% coefficient coefficient S coefficient coefficient = coefficient

term of [P] of [Q] of [P . Q] of [A . Q]

etc.

of [A . P. Q] ’

The definition of inhibition constants for substrates is analogous to those for products It is an easy matter for a computer to determine these definitions from the rate equation, particularly if the terms in the numerator are placed in order and if the terms in the denominator are placed in order. Expressed as exponent vectors, some of the kinetic constants of our example are defined below (where a, b stand for products; p, q, r stand for substrates). Mathematical

Biosciences

6 (1970), 507-517

COMPUTER DERIVATIVE

Michaelis

constants

inhibition

constants

515

:

(preferred):

~640,1,1,0)

H

a

H

(0,0,1,0,0)

(l,O,l,l,O

g

‘(l,O,l,O,O)’

(090, LO, 0) (1 0, 0, 0, 0 5

(0, 12% 0, 0)

D (0,0,2,0,0)‘(1,0,1,0,0)‘(0,1,LW)’ H

* H

pw40,0,0) (O,l,O,O,O) (l,O,l,O,O) (l,O,O, W>‘(O,130, LWW, 1, LO)’ (0, LO, 0, 0) . p

T equilibrium

P,l,O,O,

1)

’

constant: K

~u,Lmo) eq

x (1,0,0,0,0

(1,0,0,0,0)

(l,O, 1,070)

(l,o,Lo,o)x

(1,0,1,1,0) x 64% 1, LO) x

(LO,Ll,O)

(0,0,1,1,0)

(0, 0, 1, 1, 1)

= KbH,;Hg3Hl;:K, (HD, stands for the second definition of the inhibition constant associated with product P). There are additional definitions for the inhibition constants that may have to be introduced in order to recast the equation in kinetic form. In our example, one reformulation requires three additional definitions of H, and two additional definitions of H,.. Mathematical

Biosciences 6 (1970), 507-517

DONALD

516

REFORMULATION

D. FISHER

OF RATE EQUATION

AND

IN KINETIC

The procedure to rewrite the rate equation 1. generate kinetic constant definitions;

ARTHUR

R. SCHULZ

FORM

in kinetic form is

2. introduce preferred inhibition constant definitions; 3. introduce undesired inhibition constant definitions; 4. multiply numerator and denominator by 1 coefficient

of A: . . . A: ’

5. expand each denominator term in such a way as to introduce constants or their reciprocals (this may require the definition of additional inhibition constants). To illustrate the expansion technique, consider the terms (0, 0, 0, 1, l), (1, 1, 0, 1, 0), and (0, 1, 0, 1, 1) from our example. These may be rewritten as follows.

kinetic

(0,0,0,1,1)=(0,0,0,1,1) (1, l,O,O,O)

x (O,l,O,l,l)

(O,l,O,

191)

(O,l,O,

(0, l,O, 130)

190)

(O,l,O,O,O)

x (0, 1, 0, 0,O) x (1, 1, 0, 090)

(1, 1, 0, 190) = (1, l,O,O,O)

H-l

I24

or (l,O, 0, 1,O) x (1,0,0,0,0)

(l,O,O,

090)

x (l,l,O,O,O)

(O,l,O, 131)JO, 13071,1)x (O,O.0, 1, 1) (1,1,0,0,0) (O,O,O, 171) [ (1, l,O,O,O) 1 22H,-:H,,H$H,-,'K II = HzlH,-,lK a' A new definition of the inhibition constant for R is introduced reformulation of the term in (0, 0, 0, 1, 1). The term in (1, 1, illustrates two different reformulations, both of which require inhibition constant definition. The third example illustrates the Mathematical

Biosciences

6 (1970), 507-517

in the 0, 1,O) a new use of

COMPUTER

previously

517

DERIVATIVE

reformulated

expressions

in a reformulation

and

the possi-

bility of cancellations. A general flow chart for the reformulation

follows.

It should be noted

that most of these steps require several pages of flow charts and hence are not included here.

GENERAL

FLOW

CHART

derive rate equation numerator and denominator place numerator, denominator terms in ascending order generate kinetic constant definition tables read preferred definitions for inhibition constants read undesirable definitions for inhibition constants select next term from denominator apply algebraic expander to this term cancel where appropriate express result as a product of kinetic constants repeat for all terms in denominator

ACKNOWLEDGMENTS This research was supported in part by PHS grant AM08944. Computations were performed at Indiana University Medical Center’s Research Computation Center, which is supported in part by PHS grant FR00162.

REFERENCES 1 D. D. Fisher and A. R. Schulz, Math. Biosci. 4(1969), 189-200. 2 A. Schulz and D. D. Fisher, Seventh Annual Symposium on Biomathematics and Computer Science in the Life Sciences, Houston (1969). 3 A. R. Schulz and D. D. Fisher, Canad. J. Biochem. 47(1967), 889-894. 4 W. W. Cleland, Biochim. Biophys. Actn 67(1963), 188. 5 Malcolm E. Turner, Robert J. Monroe, and Henry L. Lucas, Biometrics 17(1961), 120.

Mathematical

Bioscieaces

6 (1970), 507-5 17