Numerical analysis of coupled systems of ODEs and applications to enzymatic competitive inhibition by product

Enzymatic inhibition is one of the key regulatory mechanisms in cellular metabolism, especially the enzymatic competitive inhibition by product. This inhibition process helps the cell regulate enzymatic activities. In this paper, we derive a mathematical model describing the enzymatic competitive inhibition by product. The model consists of a coupled system of nonlinear ordinary di erential equations for the species of interest. Using nondimensionalization analysis, a formula for product formation rate for this mechanism is obtained in a transparent manner. Further analysis for this formula yields qualitative insights into the maximal reaction velocity and apparent Michaelis-Menten constant. Integrating the model numerically, the e ects of the model parameters on the model output are also investigated. Finally, a potential application of the model to realistic enzymes is brie y discussed. Mathematics Subject Classi cation (2010): 45E25, 35F20, 34-XX

in cellular metabolism to produce the metabolites needed for the cells . Cells use many regulatory mechanisms to regulate the concentrations of metabolites at physiological levels. Enzymatic inhibition processes are very common mechanisms in cellular metabolism [3,4].
There are three main inhibition processes of enzymes and they are competitive inhibition, allosteric (noncompetitive) inhibition, and uncompetitive inhibition processes. In the competitive inhibition process, the enzyme molecule has one binding site for the substrate and inhibitor. Inhibitor molecules compete with substrate molecules for the binding sites of enzyme molecules. The binding of a inhibitor molecule to an enzyme molecule prevents the substrate molecules from binding to the enzyme molecule and form an enzymeinhibitor complex, so this enzyme molecule is not able to catalyse reactions to form product. Competitive product inhibition of an enzyme is an competitive inhibition process in which the product plays the role of inhibitor [3,4,5]. Mini-hexokinase I is one of the enzymes that are inhibited by their products [6,7]. Figure 1 graphically depicts the enzymatic reactions for the kinetic mechanism of competitive product inhibition of an enzyme. In the model, the enzyme molecule has one binding site that can accommodate both the substrate and the product. When the product of the active site binds to the binding site, it prevents the substrate molecules from binding to the binding site. A minimal set of chemical reactions representing competitive product inhibition is given by where E, S, and P denote an enzyme molecule, a substrate molecule, and a product molecule, respectively. The complexes ES and EP have the obvious interpretation; see Figure 1.
It is well known that the Michaelis-Menten formula for the product formation rate of enzymatic reactions is derived from a simple model for the kinetic mechanism of an enzyme, and it is a widely used tool in studying the kinetic mechanism of enzymes [5].
The remainder of this paper is organised as follows. In Section 2, we describe the formulation of the mathematical model. The non-dimensionalization analysis of the model is described in Section 3, and the eects of the model parameters on the product formation rate are given in Section 4. In Section 5, an application of the model to the phosphorylation of glucose by mini-hexokinase I is briey introduced. Finally, we nish with several conclusions in Section 6.

The mathematical model
In this section, we develop a minimal model that describes the kinetic mechanism of competitive product inhibition of an enzyme. The model is based on the mechanism described in Section 1. Although the formula for the product formation rate is already available in the literature [8,9,10], we shall derive it in a transparent manner by non-dimensionalising the equations and making rational approximations. We begin by listing our modelling assumptions.

Modelling assumptions
(i) It is assumed throughout that the mixture of substrate and enzyme is well-stirred. This implies that diusive eects in the degradation process can be neglected, and that the concentrations of the various species in the mixture can be described by functions of time only. This further implies that the evolution of the system can be modelled by a coupled system of nonlinear ordinary dierential equations, and that a partial dierential equations model is not required [11].
(ii) We assume mass action kinetics throughout; this implies that the rate of a reaction is taken to be proportional to the product of the concentrations of the reactants. We emphasize here that more complex formulas, such as the MichaelisâMenten formula for the rate of product production in an enzyme-catalysed reaction, are derivable from more fundamental mass action considerations under simplifying assumptions [5,11].

The governing ordinary dierential equations
Applying the law of mass action in the usual way, the corresponding governing ordinary dierential equations are given by where [X] = [X](t) denotes the concentration of species X at time t. It is not necessary to discuss all of these equations here. However, we do briey discuss two of them to illustrate how the governing equations are constructed. The chemical reactions for the model are displayed in Figure 1. We begin by considering the equation for S. This is given by  The enzyme here has one binding site that can accommodate both the substrate and the product. In (a) a substrate binds to a free enzyme molecule to form an enzyme-substrate complex. The enzyme then catalyses the substrate to form a product. In (b) binding of a product molecule to a free enzyme molecule to form an enzyme-product complex prevents the enzyme molecule from binding with a substrate. a -this term accounts for the reduction in concentration of S due to enzyme binding.
b -the increase in concentration of S due to enzyme unbinding from the complex ES.
Next consider the equation for the enzyme E, given by where 1 -this term accounts for the reduction in concentration of E due to enzyme binding to the substrate S.
2 -the reduction in concentration of E due to enzyme binding to the product P .
3 -the increase in concentration of E due to enzyme catalysing the complex ES and releasing the product then.
4 -the increase in concentration of E due to enzyme unbinding from the complex ES.
5 -the increase in concentration of E due to enzyme unbinding from the complex EP .
The remaining equations (2b), (2c), and (2e) are interpreted similarly. These equations are to be solved subject to the initial conditions where e 0 , s 0 are positive constants corresponding to the initial concentrations of enzyme and substrate, respectively. Summing the rst three equations, (2a) + (2b) + (2c), in (2) and integrating yields which is an expression of conservation of enzyme.

Nondimensionalization analysis
We non-dimensionalize the equations by introducing the dimensionless variables The governing equations may be written in the equivalent dimensionless form are dimensionless parameters. We have omitted the equation for e here since this can be determined from the dimensionless form for (4), given by e + c 1 + c 2 = 1.
These equations are solved subject to the initial conditions e(t = 0) = 1, In applications, the amount of substrate initially present typically greatly exceeds the enzyme present, so that e 0 s 0 . Hence, it is of value to consider the behavior of (4),(6) in the limit ε → 0. There is an initial transient at τ = O(ε) as ε → 0, but this behavior is of limited practical interest, and its discussion is omitted here. For τ = O(1), we have at leading order as ε → 0 that (see (4)) and these expressions may be manipulated to give Substituting these expressions into (4) gives Reverting to dimensional variables, the rate of formation of product is now given by and using (7), this leads to where Here V max = k 0 e 0 is the maximum production rate for the enzyme, K m is the Michaelis constant for the enzyme in the absence of product inhibition, and K D is the dissociation constant for product binding to the enzyme. Notice that we can write (9) as where is the apparent Michaelis-Menten constant that takes account of competitive product inhibition. It is noteworthy here that the maximal production rate for the enzyme V max is unaected by product inhibition. However, the apparent Michaelis-Menten constant increases linearly with product concentration; see Figure  2. Figure 3 shows the Lineweaver-Burk plots of the product formation rate formula (10) for [P ] = 0 and one value with [P ] > 0. It can be seen that the maximal product formation rate does not depend on the product concentration, while the slope of the line corresponding to [P ] > 0 is greater than that of the line corresponding to [P ] = 0, as would be expected since the presence of product slows the speed of product formation. Figure 2: Plots of the rate of product formation for the case of competitive product inhibition. The relevant formula is given in equation (9), and the plots shown illustrate the eect of the product concentration on the product formation rate. The parameter values used to generate these plots are given by Vmax

The eects of the model parameters on the model output
In this section, we numerically integrate the model to study the eects of the model parameters on the output by using the odeint solver [12] of the module integrate of the Scipy library [13] of Python [14]. The product concentration is considered as the model output. The initial conditions used here are given by: To investigate the eects of the catalytic rate constant k 0 on the output, we numerically solved the model for the three values k 0 , 1.3k 0 , and 0.7k 0 , while the values of other parameters are xed. The same process was applied to the rest of parameters. Figure 4 shows the output for the dierent values of k 0 . The blue line corresponds to the value k 0 = 1000.0 s −1 , the red dashed line corresponds to the value k 0 = 1300.0 s −1 , and the green dashed-dotted line corresponds to the value k 0 = 700.0 s −1 . It can be seen that the output is proportional to the value of k 0 , that is the higher the value of k 0 is, the higher the output is.   In Figure 5, we plot concentrations of product for the dierent values of k 1 . The blue line corresponds to the value k 1 = 10000.0 mM −1 s −1 , the red dashed line corresponds to the value k 1 = 13000.0 mM −1 s −1 , and the green dashed-dotted line corresponds to the value k 1 = 7000.0 mM −1 s −1 . It can be seen that the output is proportional to the value of k 1 , that is the higher the value of k 1 is, the higher the output is.  It can be seen that the output is inversely proportional to the value of k 2 , that is the higher the value of k 2 is, the lower the output is. This implies that reducing the binding rate constant of the product to the enzyme may be a potential approach to accelerate the rate of the enzymatic reaction.
In Figure 7, we plot concentrations of product for the dierent values of k −2 . The blue line corresponds to the value k −2 = 1000.0 s −1 , the red dashed line corresponds to the value k −2 = 1300.0 s −1 , and the green dashed-dotted line corresponds to the value k −2 = 700.0 s −1 . It can be seen that the output is proportional to the value of k −2 , that is the higher the value of k −2 is, the lower the output is. This implies that increasing the unbinding rate constant of the product from the enzyme-product complex may be a potential approach to accelerate the rate of the enzymatic reaction.
In the next section, we provide a brief introduction of a potential application of the model to a real enzyme.

Model and glucose phosphorylation by mini-hexokinase I
Glucose glycolysis is a key pathway for the production of energy in a cell, and glycolytic intermediates form precursors for the biosynthesis of other key cellular constituents, such as glycogen, nucleotide sugars, and hyaluronan. The rst step of glycolysis is the transformation of glucose into glucose-6-phosphate. This is achieved via a phosphorylation that is catalysed by an enzyme called hexokinase. There are four isozymes of hexokinase found in mammalian tissue [15,16], and these are usually referred to as hexokinase I, II, III, and IV (glucokinase). The product of glucose phosphorylation, glucose-6-phosphate (G6P ), inhibits the activity of hexokinase I, II, and III [17]. Only the C domain of hexokinase I contains the catalytic site, whereas the N domain does not [6,17,18]. Hexokinase I has binding sites for adenosine triphosphate (AT P ), glucose, and G6P in both N and C domains [6,7]. Furthermore, the binding sites for AT P and G6P have a common part, so they compete together for binding sites [6]. Figure 8 represents a single hexokinase I molecule, with blue being used for the N terminal domain and green for the C terminal domain. Each hexokinase I molecule possesses binding sites for glucose, AT P , and G6P in the both C and N domains, even though the C domain only is catalytically active [17,18,19]. In Figure 8, the binding sites for glucose on the C and N domains are depicted by the shape, and the binding sites for AT P , glucose-6-phosphate, and inorganic phosphate are represented by a ∨ cleft. In 1969, J. Ning, D. L. Purich, and H. J. Fromm [20] proposed a random Bi Bi kinetic mechanism for hexokinase I. This mechanism can be represented by the following set of chemical equations where here E, A, B, C, D represent hexokinase enzyme, AT P , glucose, adenosine diphosphate (ADP ), and G6P , respectively. Moreover, K X , K X with X = A, B, C, D are the dissociation constants for the four species. Finally, k 1 and k 2 are forward and backward rate constants, respectively, for the catalytic reaction. Many investigators have found experimental evidence in support of the Bi Bi mechanism for hexokinase I [17,21]. The model just described may be applicable to a mini-hexokinase I system. A mini-hexokinase molecule consists of the C terminal half only of the hexokinase I enzyme [19,22], and this corresponds to the enzyme species E in our model above; see Figure 9. The active site here is the binding site for AT P , so that AT P corresponds to the substrate S in the model. The product P here is G6P since G6P competes with AT P for the AT P binding site. However, the correspondence between the model and the mini-hexokinase I system falls down here since G6P is not a direct product of AT P binding -recall that a glucose molecule must also be bound to its site in the C terminal domain in order for G6P to be formed. However, G6P would be the eective product of AT P binding if AT P binding is the rate-limiting step for product formation. This would be the case for suciently high concentrations of glucose, for example. The model also does not take account of phosphate binding, and so would only apply to the mini-hexokinase system if phosphate concentrations are suciently low. However, in these circumstances, the rate of production of G6P may be approximated by (see (9)

Conclusions
Enzymatic competitive inhibition is one of the key regulatory mechanisms that cells use to regulate the concentrations of metabolites in physiological levels, especially the competitive product inhibition of an enzyme. Hence, from the point of view of applications, the development of reliable mathematical models for the enzymatic competitive inhibition by product is clearly desirable. Although the formula for the product formation rate is already available in the literature [8,9], we derive it in a transparent manner by nondimensionalizing the equations and making rational approximations in this paper. The eects of the model parameters on the model output were investigated using a Python software to numerically integrate the model. The model was seen to model, under appropriate conditions, the phosphorylation of glucose by minihexokinase I. Also, the model may be applied to describe the reaction of β-galactosidase [24]. The system under consideration may be further studied in the form of fractional dierential equations, such as Caputo fractional dierential equations [25,26].