Modelling the Phosphorylation of Glucose by Human hexokinase I

: In this paper, we develop a comprehensive mathematical model to describe the phosphorylation of glucose by the enzyme hexokinase I . Glucose phosphorylation is the ﬁrst step of the glycolytic pathway, and as such, it is carefully regulated in cells. Hexokinase I phosphorylates glucose to produce glucose-6-phosphate, and the cell regulates the phosphorylation rate by inhibiting the action of this enzyme. The cell uses three inhibitory processes to regulate the enzyme: an allosteric product inhibitory process, a competitive product inhibitory process, and a competitive inhibitory process. Surprisingly, the cellular regulation of hexokinase I is not yet fully resolved, and so, in this study, we developed a detailed mathematical model to help unpack the behaviour. Numerical simulations of the model produced results that were consistent with the experimentally determined behaviour of hexokinase I . In addition, the simulations provided biological insights into the abstruse enzymatic behaviour, such as the dependence of the phosphorylation rate on the concentration of inorganic phosphate or the concentration of the product glucose-6-phosphate. A global sensitivity analysis of the model was implemented to help identify the key mechanisms of hexokinase I regulation. The sensitivity analysis also enabled the development of a simpler model that produced an output that was very close to that of the full model. Finally, the potential utility of the model in assisting experimental studies is brieﬂy indicated.


Introduction
Glucose is a major source of energy for most living organisms. Glucose glycolysis is a key pathway for the production of energy in a cell [1], and glycolytic intermediates form precursors for the biosynthesis of other key cellular constituents, such as glycogen, nucleotide sugars, and hyaluronan. The first 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 [2,3], and these are usually referred to as hexokinase I, II, III, and IV (glucokinase). The molecular weights for hexokinase I, II, and III are all approximately 100 kDa. However, hexokinase IV is a smaller molecule, with a molecular weight of approximately 50 kDa [4].
Previous studies have identified some of the functions and expression levels for the various hexokinase isoforms. Hexokinase I is present in all tissues, where it regulates the rate-limiting step of glycolysis; the mechanism of this regulation forms the topic of the current study. It is the predominant form present in brain cells and red blood cells [5,6]. Hexokinase II is known to be highly expressed in skeletal muscle and adipose tissue [7,8].
Hexokinase III is typically present at low levels in most tissues, with the highest levels being found in the lung, the kidney, and the liver [9][10][11]. Finally, glucokinase is primarily expressed in hepatocytes and pancreatic β cells [12,13].
Abnormal levels of hexokinase in cells are known to be linked to a number of diseases. Hexokinase I deficiency is known to be associated with hemolytic anaemia; see [14,15] and the references therein. Hemolytic anaemia is a form of anaemia caused by an abnormally high level of red blood cell destruction in the body. On the other hand, overexpression of hexokinase I is associated with increased levels of basal insulin release [16]. For a therapeutic effect, the results of [17] suggest that overexpression of hexokinase II in endothelial cells can help reverse microvascular disease.
Hexokinase I and II can bind to the outer membrane of mitochondria, a process that has been associated with the prevention of cell death [18,19]. Hexokinase III does not bind to mitochondria and exists predominantly in the cytoplasmic fraction, although there is evidence for hexokinase III perinuclear binding [20]. Hexokinase III overexpression has been associated with a reduction in cell death [21]. As hexokinase III, hexokinase IV (glucokinase) cannot bind to mitochondria and is localised in the cytoplasm, where it plays a key role in the regulation of glucose homoeostasis [22].
The product of glucose phosphorylation, glucose-6-phosphate (G6P), inhibits the activity of hexokinase I, II, and III (but not glucokinase) at physiological levels. Inorganic phosphate (P i ), however, antagonises the inhibition of hexokinase I by glucose-6-phosphate at low concentrations (few millimolar) and becomes an inhibitor of hexokinase I at high concentrations. In addition, inorganic phosphate inhibits hexokinases II and III at all concentrations [11,21,23]. Only the C-terminal half of hexokinase I contains the catalytic sites, whereas the N-terminal half does not [23,24], but is involved in the P i -antagonism of the product inhibition [24,25]. In contrast, both the Cand N-terminal halves of hexokinase II are catalytically active and sensitive to G6P levels [26,27]. Furthermore, both hexokinase I and II have binding sites for ATP, glucose, G6P, and P i in both Nand C-terminal halves [4,24,28]. Similar to hexokinase I, only the C-terminal half of hexokinase III is catalytically active [29,30]. A detailed description of the kinetic mechanism for hexokinase I is given in the next section.
In this study, we for the first time synthesize the current state of knowledge of the cellular regulation of hexokinase I into a comprehensive mathematical model. The model incorporates numerous bound states for the enzyme, together with their associated activation status. Numerical simulations of the model produced results that were consistent with the experimentally determined behaviour of hexokinase I. Simulations of the model provided biological insights into the abstruse enzymatic behaviour, such as the dependence of the phosphorylation rate on the concentration of inorganic phosphate or the concentration of the product glucose-6-phosphate. A global sensitivity analysis of the model was implemented to help identify the key mechanisms of hexokinase I regulation. The sensitivity analysis enabled the development of a simpler model that produced an output that was very close to that of the full model. Finally, the potential utility of the model in assisting experimental studies is indicated.

Mathematical Model
Many cellular factors can influence the phosphorylation of glucose by hexokinase I. In the current paper, we constructed a mathematical model that describes the cellular regulation of glucose phosphorylation. One of the principal aims of the modelling was to gain insight into the roles of G6P and P i in regulating the phosphorylation process. The model consisted of a system of ordinary differential equations that tracks the evolution in time of the concentrations of various relevant species, including hexokinase enzyme, glucose, G6P, ATP, ADP, and P i . We give schematic representations for each of these species in Figure 1. Figure 1a represents a single hexokinase I molecule, with blue being used for the Nterminal domain and green for the C-terminal domain. Each hexokinase I molecule possesses binding sites for glucose, ATP, G6P, and inorganic phosphate in the both the Cand Ndomains, even though the C-domain is only catalytically active [4,23,24]. In Figure 1a, the binding sites for glucose on the Cand N-domains are depicted by the shape and the binding sites for ATP, glucose-6-phosphate, and inorganic phosphate are represented by a ∨ cleft. In 1969, Ning et al. [31] proposed a random Bi Bi kinetic mechanism for hexokinase I. This mechanism is represented by the chemical equations shown in Figure 2. There is considerable experimental evidence in support of the Bi Bi mechanism for hexokinase I; see [4,[32][33][34][35]. In the context of the current study, it forms a part of a larger kinetic mechanism we developed for hexokinase I. Here, E, A, B, C, and D represent the hexokinase enzyme, ATP, glucose, ADP, and G6P, respectively. K X and K X are the dissociation constants for X = A, B, C, D.
The mathematical model we developed for the phosphorylation of glucose is necessarily large and complex since it describes multiple binding sites, numerous species, almost 150 chemical reactions, and various inhibitory mechanisms. The phosphorylation mechanism of glucose by hexokinase I was already briefly introduced in the Introduction Section, and schematic representations of the relevant species arising are given in Figure 1. Recall that each hexokinase I molecule has two subunits, an Nand a C-terminal domain. Each subunit has its own binding site for glucose and another binding site for ATP, P i , and G6P [4]. In Figure 3, we depict the eight possible configurations of the molecule where only one of the binding sites is occupied.

The Kinetic Mechanism
The kinetic mechanism for the phosphorylation of glucose by hexokinase I may be summarized as follows:
Product production: The product here is G6P, and it is produced by the phosphorylation of glucose by hexokinase I. The phosphorylation is achieved via a Bi Bi mechanism [31], as represented in Figure 4a. However, in the current study, we considered the simplified process represented in Figure 4b. This is partially justified by noting that the simplified mechanism produces the same ultimate products as the Bi Bi mechanism, that is E, ADP, and G6P. This is because all of the reactions in the right branch of the Bi Bi mechanism are essentially irreversible; see Section 16.2 of [43]. Here, for a G6P molecule to be produced, an ATP molecule must be bound to its C-domain site and a glucose molecule must be bound to its C-domain site; see Figure 5;

3.
Product production is regulated: The phosphorylation of glucose is inhibited via the following three mechanisms: (a) Competitive inhibition: P i competes with ATP for the C-domain binding site; see Figure 3f.
(b) Allosteric product inhibition [44]. Binding of a molecule of G6P to the Nbinding site makes a conformational change to the C-domain binding site for ATP. This conformational change disables the binding of ATP to the Cdomain, resulting in the deactivation of the enzyme [37,38]; see Figure 3g. This inhibition is mitigated by the presence of P i and ATP, which compete with G6P for the N-domain binding site; (c) Competitive product inhibition: The product G6P competes with ATP for its C-domain binding site, inhibiting product production; see Figure 3h;

4.
Other details: The following information concerning glucose phosphorylation is also available in the literature: (a) Only one molecule of G6P can bind to an enzyme molecule at a time [37,38]; The binding of P i to the N-domain binding site weakens the binding of G6P to the C-domain binding site (that is, it increases the dissociation constant) [39,41]; (c) The ATP binding sites of free or complexed enzyme are open, except for the case where a G6P molecule is bound at the N-binding site; (d) The high-affinity binding site for G6P is in the C-domain, while the high affinity binding site for P i is in N-domain [39].

Modelling Assumptions
(a) It was assumed throughout that the cellular mixture of hexokinase I, glucose, ATP, and P i is well-stirred. This implies that diffusive effects in the phosphorylation 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 differential equations and that a partial differential equation model is not required; (b) We assumed mass action kinetics throughout [45]; this implies that the rate of a reaction is taken to be proportional to the product of the concentrations of the reactants. We emphasise here that more complex formulae, 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 [46]; (c) We focused our attention solely on the phosphorylation of glucose and made no attempt to model in detail the evolution of the intracellular glucose concentration. Rather, we assumed instead a constant initial concentration of glucose and used the model to track its subsequent depletion as it is converted to G6P via phosphorylation; (d) The mechanism of the phosphorylation of glucose was assumed to be the simplified Bi Bi process represented in Figure 4b; (e) The binding of one substrate does not affect the affinity of the binding sites for other substrates, except that the binding of P i at the N-domain reduces the affinity of the C-binding site for G6P; (f) The model allows only one molecule of G6P to bind to an enzyme molecule at a time.

Model Notation
We introduce the following model notation. We write: E : a hexokinase I molecule; G : a glucose molecule; T : an ATP molecule; G6 : a G6P molecule; P : a P i molecule; D : an ADP molecule.
We also added subscripts and superscripts to E, where a subscript denotes a molecule binding to the C-domain of the enzyme and a superscript denotes a molecule binding to its N-domain. Hence, for example, we have:  Figure 7. Diagram of chemical reactions forming all complexes in the mixture. E, 0, 1, 2, and 3 represent the free enzyme, glucose, ATP, G6P, and P i molecules, respectively. The red lines refer to reversible reactions involving glucose binding; the blue: ATP binding; the green: G6P binding; the brown: P i binding. Letters E with a superscript(s) and/or subscript(s) denote complexes of the enzyme, for instance E 01 3 is a complex of the enzyme with one glucose and one ATP molecule bound at the N-binding site, and one P i molecule bound at the C-binding site. We introduce the following notation for the model rate constants. k 0 : the catalytic constant (turnover rate) for hexokinase I; k G , k T , k G6 , k P : the forward rate constants for the binding of glucose, ATP, G6P, and P i , respectively, to their N-binding sites; the reverse rate constants for the dissociation of glucose, ATP, G6P, and P i , respectively, from their N-binding sites; k G , k T , k G6 , k P : the forward rate constants for the binding of glucose, ATP, G6P, and P i , respectively, to their C-binding sites; the reverse rate constants for the dissociation of glucose, ATP, G6P, and P i , respectively, from their C-binding sites; k P G6 , k P −G6 : the forward and reverse rate constants for the binding/unbinding of G6P to/from its C-binding site, when the enzyme has a P i molecule bound to its N-binding site.

The Chemical Reactions
The system has numerous chemical reactions because of the large number of possible bound states for the enzyme, and a small selection of these is given by: There are 147 reactions in all, and these are listed in detail in the Supplementary Material. In Figure 7, we schematically represent all of the chemical reactions producing an enzyme complex. Figure 6 shows the reactions that lead to the production of the product G6P.

Construction of the Governing Ordinary Differential Equations
The complete set of governing equations for the model can be found in the Supplementary Material, but it is not necessary to discuss all of these equations here. However, we do briefly discuss two of them to illustrate how the governing equations are constructed. The equations were developed based on the chemical reactions referred to in the previous subsection and the law of mass action. We begin by considering the equation for the product G6P, which is given by: where: a These terms account for the increase in the concentration of G6P due to the creation of product by enzymatic reactions involving the complexes G E T , G G E T , G E i T , and G G E i T , where i = T, G6, P; see Figure 6; b The increase in the concentration of G6P due to G6P unbinding from the enzyme The increase in the concentration of G6P due to G6P unbinding from the enzyme The increase in the concentration of G6P due to G6P unbinding from the complexes E P G6 , G E P G6 , G E P G6 , and G G E P G6 ; see Figure 7; e The reduction in concentration of G6P due to G6P binding with the species E, G E, G E, Figure 7. Next, consider the equation for the enzyme concentration, given by: where: α This term accounts for the increase in the concentration of the enzyme due to the recovery of the enzyme after the catalytic step has completed; β This gives the increase in the concentration of the free enzyme due to the unbinding from the complexes G E, G E, E i , and E i , where i = T, G6, P; γ This gives the reduction in the concentration of the enzyme due to enzyme binding.

Initial Conditions
The equations described in the previous subsection were solved subject to the initial conditions: where E 0 , G 0 , ATP 0 , and P i0 give the initial constant concentrations of the enzyme, glucose, ATP, and P i , respectively. The initial concentrations for all of the enzyme complexes were taken to be zero.

Computational Methods
In this section, we describe the computational tools used to analyse the model equations. The software developed for this paper was coded using the Python programming language [47].

Numerical Method for Solving the Ordinary Differential Equations
The system of differential equations was numerically integrated using the odeint solver in the module integrate of the SciPy library. SciPy [48] is an open-source Python library that contains numerical routines for applications in science and engineering. The odeint solver [49] uses the LSODA program [50] from the FORTRAN library odepack, and it is capable of solving both stiff and nonstiff systems. Table 1 shows some of the model parameter values, together with their literature sources. We note that the dissociation constant for P i at its C-binding site was taken to be ten-times larger than its value at the N-site [51]. This implies that the higher affinity binding site for P i is in the N-domain, with much weaker binding to the C-site. We recall that the enzyme is active when P i is bound to its N-site (Figure 3e), but inactive when bound to the C-site (Figure 3f). Hence, for low concentrations of P i (few millimolar), the higher affinity N-site dominates and the inhibition of G6P is antagonised. However, for higher P i concentrations, P i binding to the C-site is significant and enzyme activity is inhibited. This behaviour matches experimental findings [4,11,41]. C 220 [51] The Lambda (Λ) and Omega (Ω) methods for approximating kinetic rate constants were discussed in the paper [52]. Rate constants for the current model were estimated using these methods with Λ = 100, Ω = 1.0 and the data displayed in Table 1; see Table 2. The Michaelis-Menten constants for the C-domain binding sites of glucose and ATP are known [37,39]. The corresponding values for the N-domain are unknown and, so, in the absence of other information, were taken here to be the same as their C-domain values.

Global Sensitivity Analysis
A Sobol global sensitivity analysis was implemented to evaluate the importance of the various parameters appearing in the model [56,57]. The model considered in this paper has the structure y = f(q, t), where the inputs are the model parameters q and time t and the output y is a vector that gives the model predictions for the concentrations of the various species at time t. A Sobol sensitivity analysis enabled us to quantify how variations in the model parameters q = (q i ) affect the model output y. This was achieved via the calculation of sensitivity indices.
We now briefly describe the sensitivity indices. The discussion given here largely mirrors that given in Mai et al. [58] and is reproduced here for the convenience of the reader. For the parameter q i , the associated first-order sensitivity index S i is given by: where D i is the variation of the model output with respect to changes in the parameter q i and D is the variation in the model output with respect to changes in all of the model parameters q. For brevity, we do not explicitly display the formulae for D i and D here; the details can be found in [57]. These first-order indices represent the effect of an individual parameter q i on the output without interactions with the other parameters. For the pair of parameters q i and q j (i = j), the associated second-order sensitivity index S ij is given by: where D ij is the variation of the model output with respect to changes in the parameters q i and q j . This index measures the effect of the interaction between the parameters q i and q j on the model output. These ideas generalise in an obvious way for a set of parameters q i , q j , ...., q k , where we define the sensitivity index: and where D ij...k is the variation of the model output with respect to the parameters q i , q j , ...., q k . The index S ij...k measures the effect of the interaction between the parameters q i , q j , ...., q k on the model output.
The total sensitivity index, S tot i , is the sum of all of the indices involving the parameter q i , without repetition. It gives a measure for the total effect of the parameter q i . Rather than display a rather opaque general formula for S tot i [57], it is more instructive here to illustrate the idea with particular examples. If there are three parameters in total q 1 , q 2 , q 3 , then: For four parameters q 1 , q 2 , q 3 , q 4 , we have: S tot 2 = S 2 + S 12 + S 23 + S 24 + S 123 + S 234 + S 124 + S 1234 .
Notice here that we have included S 12 , but not S 21 = S 12 , and so on. If it is found that the values of S tot i and S i are close, then the higher-order indices are small, and this implies that interactions between the parameter q i and the other model parameters do not significantly affect the model output.
The sensitivity analysis in the current study was implemented computationally using the Python package SALib [59,60].

Numerical Results
The Methods Section introduced the mathematical model. It also described the computational methods used to integrate the model equations and included some discussion of the numerical method, the choice of parameter values, and the initial conditions. In the current section, we describe some of the numerical results obtained.
The principal purpose of the numerical solutions displayed here was to gain insight into the cellular phosphorylation of glucose by hexokinase I. To focus attention on the phosphorylation process itself, we made no attempt to model the evolution of intracellular glucose levels. Instead, we simply assumed a constant initial concentration of glucose and then tracked its subsequent conversion via phosphorylation to G6P. For similar reasons, we also made no attempt to model the cellular behaviour of G6P subsequent to its production. Figure 8 shows the time evolution of the concentration of G6P for a range of different initial concentrations of P i and G6P. In the numerical results displayed here, the initial concentration of glucose was taken to be 2.5 mM [54], and the total concentration of the enzyme was taken to be 6.65 × 10 −2 mM [53]. We begin by noting some broad features of the behaviour exhibited in Figure 8. All of the curves are increasing functions of time, as would be expected since G6P levels increase as the available glucose is phosphorylated. The levelling off of the curves corresponds to the exhaustion of the available glucose substrate. It is also noteworthy that the time scale over which the phosphorylation process was completed was of the order of ten seconds, a prediction that is consistent with the literature values [61,62]. In Figure 8, Subplots (a), (b), (c), and (d) correspond to differing initial concentrations of G6P, with (a) having the lowest initial concentration and (d) the highest. It is clear that the time for the phosphorylation process to be completed increased with increasing initial concentration of G6P. This again was as expected since G6P is the species that is responsible for both the allosteric and competitive product inhibition of the enzyme, and so increasing its concentration should slow the phosphorylation process.
The dependence of the phosphorylation behaviour on the initial concentration of inorganic phosphate is more subtle and interesting. Focus, for example, on the subplot Figure 8d, and begin by considering the curve corresponding to [P i ] = 0. This is the curve corresponding to zero initial P i concentration, and it gives a convenient reference. We note that for relatively low P i concentrations (1 mM, 2 mM), the phosphorylation process was faster than the phosphate-free case. However, for the higher concentrations (P i ≥ 10 mM), we note that the phosphorylation rate was slower relative to the phosphate free case. This is in line with experimental findings [4,11,41], which showed that for low concentrations of phosphate (few millimolar), enzyme inhibition by G6P is antagonised and that for higher phosphate concentrations, enzyme activity is inhibited. In the context of the current modelling, this behaviour is explained by recalling that for low concentrations of P i , the higher affinity N-binding site for P i dominates and inhibition by G6P is antagonised. However, for higher P i concentrations, P i binding to the lower affinity C-site is significant and enzyme activity is inhibited. This phenomenon is clearly exhibited in Figure 9, where we show plots of the phosphorylation rate as a function of inorganic phosphate P i concentration and for four different initial concentrations of G6P.  Tables 2 and 3. We note an initial increase in the phosphorylation rate in all cases, followed by a subsequent decrease in the rate. This is discussed further in the main text.

Results of the Global Sensitivity Analysis
The Sobol global sensitivity analysis employed in the current study was described in the Methods Section, and it was implemented using the SALib package [59]. A Sobol analysis enabled us to quantify how variations in the model parameters: affect the model output. The parameters here are the model rate constants and are described in Table 2. The model is of the form y = f(p, t), where the inputs are the parameters p and the time t, and the output y gives the predictions for the concentrations of the model species at time t. In the current study, we confined our attention to the output for the G6P concentration since G6P is the product here. Default settings for the SALib package were used, with one exception: no secondorder indices were calculated [59]. We calculated first-order and total sensitivity indices for the times t = 2, 4, 6, 8, 10 s. The output was then G6P(t = 2, 4, 6, 8, 10 s), and the purpose of the analysis was to evaluate the sensitivity of this output to variations in the parameters using the sensitivity indices. The G6P values were calculated by numerically integrating the governing ordinary differential equations, as previously described. The initial concentration for G6P was taken to be 2.0 mM when numerically integrating the differential equations. Two choices for the initial concentration of phosphate were made, 2.0 mM (low P i concentration) and 10.0 mM (high P i concentration). The initial concentrations for the enzyme, glucose, and ATP are given in Table 3. The remaining species had zero initial concentration. Figure 10 displays first-order sensitivity indices and total sensitivity indices for the model parameters. Figure 10a,b shows values for a low initial phosphate concentration (2.0 mM), while Figure 10c,d gives values for a high initial phosphate concentration (10.0 mM). For convenience, we split the model parameters into three groups: Group I: The sensitivity indices for all of the Group I parameters are small. This means that the rate of G6P production is relatively insensitive to modest variations in the assumed values of these parameters. The parameters in Group I describe, among other things, the rate of binding and unbinding of glucose to both the Nand the C-domains of the enzyme and the rate of binding/unbinding of P i to the N-domain of the enzyme.
We now turn our attention to the parameters in Group III. The first-order sensitivity indices for all of these parameters are relatively large, implying that the G6P production rate is relatively sensitive to variations in the assumed values of these parameters. The Group III parameters determine the turnover rate of the enzyme, the rate of binding/unbinding of ATP to the Cdomain of the enzyme, and the rate of binding/unbinding of G6P to the C-domain of the enzyme when a P i molecule is bound at its N-site. It is noteworthy that the indices S 9 and S −9 are quite sensitive to the phosphate concentration, being significantly larger for lower initial phosphate concentration. Hence, for low phosphate concentrations, the binding of P i to the N-binding site is one of the key regulators of enzyme activity.
The parameters for Group II (k −P , k P ) are also seen to be sensitive to the phosphate concentration, being small for the low phosphate cases and significant for higher phosphate cases. These parameters determine the rate of binding/unbinding of P i to its C-binding site.
It is clear from Figure 10 that the first-order sensitivity indices are close to the total sensitivity indices. Hence, interactions between any of the parameters k i and the other model parameters do not significantly affect the model output.

Further Numerical Results
We now display some further numerical solutions inspired by the results of the sensitivity analysis just presented. In these calculations, the default values used for the k parameters are given in Table 2, and the initial concentrations for the enzyme, glucose and ATP are given in Table 3. The initial concentrations of G6P and P i were taken to be 3 mM and 6 mM, respectively, corresponding to a high phosphate concentration case.
We illustrate how the new numerical results shown in Figure 11 were generated by considering a particular example. Consider the curves displayed in Figure 11a. The solid red curve was generated using the default values for the k parameters. The dashed blue curve was generated using the default values, except that 1.3k 0 was used rather than k 0 , and the dash-dotted green curve was generated using 0.7k 0 rather than k 0 . Hence, the three curves shown in Figure 11a help evaluate the sensitivity of the model output to variations in the parameter k 0 . The remaining subplots in Figure 11 were generated by repeating this process for the parameters k T , k −T , k P −G6 , k G , and k G6 . The results shown in Figure 11 are consistent with the predictions of the sensitivity analysis since the model output is seen to be quite sensitive to the parameters with relatively large sensitivity indices (k 0 , k T , k −T , k P −G6 ), but insensitive to the parameters with small indices (k G , k G6 ).

Model Reduction
Motivated by the results of the sensitivity analysis, we now consider a simplified model (SM). This model was obtained from the full model by setting: In this SM, glucose molecules do not bind to their N-domain site and ATP molecules do not bind to their N-domain site. This reduction removes 46 chemical reactions from the model system and reduces the number of governing equations by 24. We previously saw that the sensitivity indices for all of these parameters are small, and so we anticipated that the output for these models should typically closely match that for the full model. The initial conditions used to generate the numerical solutions for the SM were the same as those used for the full model, with the exception of the initial conditions for G6P and P i , which are specified on the numerical figures. Figure 12   It is seen in Figure 12 that the results of the full model and the SM are close in all cases, suggesting that the full model may reasonably be simplified by dropping the mechanisms of glucose and ATP binding to the N-domains.

Conclusions
In this paper, we developed a comprehensive mathematical model describing the phosphorylation of glucose by the enzyme hexokinase I. Glucose phosphorylation is the first step of the glycolysis pathway, and so, it is carefully regulated by cells. The regulation of hexokinase I is quite complex and includes three inhibitory mechanisms: a competitive product inhibitory mechanism, an allosteric product inhibitory mechanism, and a competitive inhibitory mechanism. We used the mathematical model to help unpack the regulatory behaviour of hexokinase I. We make the following remarks: • Model novelty: The biological evidence underpinning the modelling developed here has been established elsewhere. Nevertheless, the current study represents the first time that these ideas have been synthesised into a detailed mathematical model. The model incorporates numerous bound states for the enzyme, together with their associated activation status; the potential practical significance of this from the point of view of experimental work is explained in the next point. It is also noteworthy that the current modelling does not make any equilibrium assumptions; • The mathematical model and experimental data: The model presented here may be of particular value when combined with an appropriate in vitro experimental approach.
In particular, the mathematical model may aid in identifying the correct biological mechanism from amongst a set of competing biological mechanisms. For example, two competing models for the function of hexokinase I are the Fromm model and the Wilson model, as described in Chapter III of [63]. Since the model developed here is sufficiently refined, both suggestions may be simulated by choosing appropriate parameter values in the model. Comparing both simulated curves with appropriate experimental data may assist in selecting the correct mechanism; see Figure 13; • Biological insight: The model was numerically integrated using the SciPy Python library, and the solutions obtained were found to be consistent with the known behaviour of hexokinase I. The simulations also provided biological insights into subtle aspects of the enzyme behaviour, such as the dependence of the phosphorylation rate on the concentration of inorganic phosphate or the concentration of the enzyme product G6P. Further insight was obtained by carrying out a global sensitivity analysis, which revealed the key mechanisms of hexokinase I regulation. The results of this analysis showed that the rate of phosphorylation is sensitive to the following factors: the turnover rate of the enzyme; the rate of binding/unbinding of ATP to/from the C-domain of the enzyme; the rate of binding/unbinding of G6P to/from the Cdomain of the enzyme with a P i molecule bound at the N-domain for low phosphate concentration; and the rate of binding/unbinding of phosphate to/from the C-domain of the enzyme for high phosphate concentration; • Simplified model: One reduced model was developed based on the results of the sensitivity analysis. This simpler model produced results that closely matched the results of the full model; • Software. The software developed in this study to numerically integrate the governing equations and to implement the sensitivity analysis has been made available in the Supplementary Material and online (https://github.com/vinh-mai/Hexokinase-2019 accessed on 1 August 2021). Although the model developed in the current study is comprehensive and detailed, there is scope for improvement and future work. For example, one avenue for future research is to incorporate the full detail of the Bi Bi mechanism in the modelling. Furthermore, glucose-6-phosphate binding to its N-binding site not only allosterically inhibits the enzyme, but also stimulates enzyme release from mitochondria [64,65], and the incorporation of this effect provides another interesting avenue for future study. Finally, the possible inhibition of hexokinase I by ADP [51] was not explored in the current study.