# The methionine metabolism in human liver cell: steady-state sensitivity analysis by biochemical systems theory.

Introduction

A true understanding of metabolic functions and design is crucially dependent on mathematical and computational methods for characterizing the system in concern. The concept of an effective mathematical tool called sensitivity analysis is widely applied in variety of fields. From the late 1960's, the need for description of organizationally complex systems led Michael Savageau and co-workers to built a powerful framework known as Biochemical Systems Theory (BST), for systems-analysis of biochemical processes ([1] to [6]) that offers an efficient mathematical procedure for the sensitivity analysis of given system.

In understanding the regulatory mechanism of the methionine-homocysteine metabolism an investigation of the full cycle through a mathematical model could be a useful approach. A simple mathematical model consisting of two differential equations for only two metabolites, AdoMet and AdoHcy was initially proposed by Martinov et al. [8] which found insufficient to explain the role of AdoMet concentration in regulating the remethylation reaction [12]. Reed et al. [12] later added several pathways in the model of Martinov et al. and extended the model to describe the concentration changes of all the metabolites in the cycle (Methionine, AdoMet, AdoHcy and Hcy).The mathematical model of Reed et al. [12] can be used to develop better understanding of the methionine cycle through an investigation of sensitivity that efficiently elucidate the characteristics of the metabolic reaction system, since, the behavior of the cycle is very complicated due to the fact that some enzymes are activated and inhibited by several intermediates of the cycle.

In the present work, therefore, we carry out sensitivity analysis of the methionine cycle by the mathematical model of Reed et al. [12] and exploit the concealed property of this system in the framework of BST. For the steady state analysis, we determine logarithmic gains and parameter sensitivities of the methionine cycle from the S-system type differential equations ([1] to [4]; [18]) derived by power-law approximations of the rate laws in the mathematical model of Reed et al. [12] and characterize the system at a nominal steady state.

Theoretical

Mathematical model of methionine metabolism

The metabolic pathway map of the methionine cycle is shown in Fig.1a. To study the salient features of its regulation properties, Reed et al. [12] formulated four differential equations (Eqs.(1)-(4)):

d[Met]/dt = [[upsilon].sub.MS] + [[upsilon].sub.BHMT] + [[upsilon].sub.Metin]-[[upsilon].sub.MATI]- [[upsilon].sub.MATIII] (1)

d[AdoMet]/dt = [[upsilon].sub.MATI] + [[upsilon].sub.MATIII]-[[upsilon].sub.METH]-[[upsilon].sub.GNMT] (2)

d[AdoHcy]/dt = [[upsilon].sub.METH] + [[upsilon].sub.GNMT]-[[upsilon].sub.AH] (3)

d[Hcy]/dt = [[upsilon].sub.AH]-[[upsilon].sub.CBS]-[[upsilon].sub.MS]-[[upsilon].sub.BHMT] (4)

All the rate equations on the right hand side are shown in Table 1 and the values of all other parameters and constants included in the rate equations are as given in Table 2.

[TABLE 1 OMITTED]

BST analysis

S-System representation

BST gives a computational method to systematically and efficiently perform the sensitivity analysis in the steady state and transient state by use of differential equations including nonlinear flux expressions in metabolic reaction systems. The S-system representation (Eq.(5)), in the formalism of BST, aggregates the processes affecting a dependent metabolite, Xi, for its synthesis and degradation separately and a product of power-law functionis used to describe the rate ([V.sub.+i] for synthesis and [V.sub.-i] for degradation) of each aggregate function. The dependent concentration variables ([X.sub.i]) are represented consecutively from [X.sub.1], [X.sub.2], to [X.sub.n] and the independent concentration variables are represented from [X.sub.n+1] to [X.sub.n+m].

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (5)

The parameters in the system are then identified readily, which are two rate constants, [[alpha].sub.i] for synthesis and [[beta].sub.i] for degradation, and two kinetic orders, [g.sub.ij](synthesis) and [h.sub.ij](degradation) associated with each aggregate rate laws (for details see [1] to [4]).

Under a steady state condition, the time derivative of each dependent variable in Eq. (5) becomes equal to zero and then the matrix notation of the resulting algebric equations gives

[y].sub.d]= -[[A].sup.-1.sub.d][[A].sub.i][y].sub.i] + [[A].sup.-1.sub.d] b]=-[M][[A].sub.i][y].sub.i] + [M]b] (6)

where [y].sub.d] represents the vector that contains logarithms of dependent variables ([y.sub.i] = log [X.sub.i]), [A] represents the matrix that contains kinetic orders of dependent and independent variables of n x (n + m) dimension, and b] represents the vector of dimension n where [b.sub.i] = log ([[beta].sub.i]/[[alpha].sub.i]). The steady state fluxes through any dependent variable (or pool), [X.sub.i], are obtained using the known steady state concentrations and the appropriate aggregate rate law [14]. The matrix expressions for the fluxes of S-system representation are thus given by

(log[V.sub.+])]=(log[alpha])] + [G]y] (7)

and

(log[V.sub._])]=(log[beta])] + [H]y] (8)

where [V.sub.+] and [V.sub._] indicate the vectors of the net flux into and from the [X.sub.i] pool, respectively, and [G] and [H] represent the matrices of the kinetic orders, [g.sub.ij] and [h.sub.i,j], respectively. The solutions of equations (6) to (8) give the steady state values of the dependent variables of a system represented in the S-system form.

The change of a system property in response to the relative change in the independent variable is called logarithmic gains. The logarithmic gains simply represent the percentage change in the logarithmic value of the steady state of a dependent variable and net flux or local flux as a result of an infinitesimal change in the logarithmic value of an independent variable [15], which are computed as follows:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (9)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

Similarly, the changes of a system property in responses to the relative changes in the rate constant and kinetic order are called rate-constant sensitivities and kinetic-order sensitivities, respectively.

S([X.sub.i], [[alpha].sub.j]) = [partial derivative] ln [X.sub.i]/[partial derivative] ln [[alpha].sub.j] = [partial derivative] [X.sub.i]/[partial derivative] [[alpha].sub.j] [[alpha].sub.j] [X.sub.i] = [-M.sub.ij] (i, j = 1, ..., n) (11)

S([X.sub.i], [[beta].sub.j]) = [partial derivative] ln [X.sub.i]/[partial derivative] ln [[beta].sub.j] = [partial derivative] [X.sub.i]/[partial derivative] [[beta].sub.j] [[beta].sub.j] [X.sub.i] = [-M.sub.ij] (i, j = 1, ..., n) (12)

Similarly, the kinetic-order sensitivities, defined as the ratio of percentage change in a dependent variable or dependent flux to an infinitesimal percentage change in a kinetic order parameter, are written as:

S([X.sub.i], [g.sub.kp]) = [partial derivative] ln [X.sub.i]/[partial derivative] ln [g.sub.kp] = ([partial derivative[X.sub.i]/[partial derivative][g.sub.kp]])([g.sub.kp]/[X.sub.i]) = [y.sub.i]S([y.sub.i], [g.sub.kp]) (13)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

S-system type equations for methionine cycle

For characterization of the methionine cycle in the framework of BST, it is necessary to assign BST variables to the dependent and independent variables. The metabolite concentrations, [Met], [AdoMet], [AdoHcy], and [Hcy], are the dependent variables and expressed by [X.sub.1], [X.sub.2], [X.sub.3], and [X.sub.4], respectively. The independent variables include the methionine input of the cycle-Metin and the maximum reaction rates, owing to their affects on the pathway and remaining unchanged by the systems dynamics. In the remethylation reaction, [Hcy] acquires a methyl group from the concentration of 5-methyl tetrahydrofolate (5mTHF) and converted into methionine (Fig.1) by methinine synthase. The concentration of 5mTHF thus enters the pathway and remains unchanged by the systems dynamics. Therefore, we include the 5mTHF as an independent variable. All the independent variables are expressed by [X.sub.5] to [X.sub.16] (see Fig.1b and Table 2). According to the assignment of variables on the map in fig.1b, transformation of Eqs.(1)-(4), in the Michaelis-Menten form, using the steady state values leads to the S-system type differential equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (18)

It should be noted that Eqs.(15)-(18) include 4 dependent variables (metabolite concentrations), 12 independent variables (maximum reaction rates), 8 rate constants for both incoming and outgoing fluxes and 19 kinetic orders for incoming and 20 kinetic orders for outgoing fluxes.

[FIGURE 1 OMITTED]

Results

In Eqs.(1)-(4) with [[upsilon].sub.Metin]=200, the derivatives of the metabolite concentrations were set equal to zero and the algebraic equations thus obtained were solved by the Newton-Raphson method to get the steady-state metabolite concentrations. As seen in Table 3, our calculated values are almost identical to those of Reed et al. [12], indicating that our calculation condition is almost the same as that of Reed et al. [12]. Moreover, this steady state is included in the normal physiological range in the experimental values of Martinov et al [8]. Thus, the steady state values were used to transform Eqs.(1)-(4) to the following S-system type differential equations:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (21)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (22)

Although all the calculations in the present study were performed with double precision of accuracy, the rate constants and kinetic orders in Eqs.(19)-(22) are given with six significant digits, which will be sufficient to roughly reproduce our calculated result. Naturally, the steady state values given by setting equal to the derivatives of [X.sub.1]-[X.sub.4] in Eqs.(19)-(22) is equal to those in Table 3.

We investigated the influences of 12 independent variables (11 enzymatic activities and [5mTHF]) on 4 dependent variables (metabolite concentrations) at a nominal steady state in Eqs.(19)-(22). The result is shown in Fig.2, where each bar of the three-dimensional plot represents the absolute value of the logarithmic gain, [absolute value of L([X.sub.i], [X.sub.j]]. The fact that L([X.sub.1], [X.sub.5]) = 0.013 indicates that the Met concentration ([X.sub.1]) hardly responds to the independent variable, [X.sub.5], i.e., the flux of the extra cellular Met into the cycle, [[upsilon].sub.Metin]. This result is identical to the calculated result of Reed et al. [12]. The [X.sub.1] pool has the maximum logarithmic gain to an infinitesimal perturbation of [X.sub.9] (= [V.sup.METH.sub.max]), i.e., L([X.sub.1],[X.sub.9]) = 0.69. Likewise, the concentrations of AdoMet, AdoHcy, and Hcy ([X.sub.2], [X.sub.3], [X.sub.4], respectively) markedly respond to changes in [X.sub.9]. These results imply that all the metabolite concentrations in the methionine cycle are influenced by the independent variables that directly determine the flux [[upsilon].sub.meth]. The largest absolute value of the logarithmic gain is L([X.sub.4],[X.sub.1]5) = -0.94 for the response of Hcy ([X.sub.4]) to an infinitesimal change in the rate constant [X.sub.1]5 ([[psi].sub.1]) for the reaction ([[upsilon].sub.CBS]) that is catalized by cystathionine [beta]-synthase. This value approximately implies that an 1% increase in [[psi].sub.1] decreases the steady state value of Hcy ([X.sub.4]) by 0.94%. Thus, it is clear that the [X.sub.4] pool is the most sensitive part in the methionine cycle. From Fig. 2, however, it is noticeable that the steady-state logarithmic gains of the methionine cycle model remain below unity. Hence, the system may not response drastically to the change in independent variables.

[FIGURE 2 OMITTED]

Robustness of the model

We calculated the rate-constant and kinetic-order sensitivities to investigate the robustness of the methionine cycle model. Figure 3.a) shows the three dimensional plot of the absolute values of the rate-constant sensitivities which indicates that the model is considerably robust in response to any amplification in the rate constant parameters, as no rate constant sensitivity exceeds a value of 2.5. The prominent rate constant sensitivities are however occur for the metabolites Hcy and AdoHcy in response to the rate-constant parameter, [[alpha].sub.3] and [[alpha].sub.4].

Fig.3.b) shows the three-dimensional plot of the absolute values of the kineticorder sensitivity values. The infinitesimal change in 19 [g.sub.ij] and 20 [h.sub.ij] give 156 kinetic-order sensitivities for the metabolite concentrations. It is obvious that Met ([X.sub.1]), which is the first metabolite, is not so sensitive to many kinetic orders, because the absolute values of the sensitivities are less than 0.2 other than S([X.sub.1], [g.sub.2,1]) = -3.84, S([X.sub.1], [g.sub.2,7]) = -4.57, and S([X.sub.1], [g.sub.2,9]) =5.97. On the other hand, AdoMet ([X.sub.2]), which is the second metabolite, is remarkably influenced by many kinetic orders; Of 39, 6 absolute values are in the range of 4-7, 13 are in the range of 1-4, and the others are less than 1. The largest sensitivity is S([X.sub.2], [h.sub.2,9]) = -6.72. On the other hand, AdoHcy ([X.sub.3]) and Hcy ([X.sub.4]) respond more sensitively to the kinetic orders. These include sensitivity values of 8-14; S([X.sub.3], [g.sub.3,9]) = 12.1, S([X.sub.3], [h.sub.3,3]) = -11.8, S([X.sub.3], [g.sub.4,3]) = 9.34, S([X.sub.4], [h.sub.2,9]) = -9.20, S([X.sub.4], [g.sub.3,9]) = 13.9, S([X.sub.4], [h.sub.3,3]) = -13.6, S([X.sub.4], [g.sub.4,3]) = 13.6, and S([X.sub.4], [g.sub.4,11]) = -8.47. The perturbations of [g.sub.39] and [g.sub.43] significantly increase [X.sub.4], whereas the perturbation of [h.sub.33] remarkably decreases [X.sub.4]. It is seen that 34% of the absolute values of the sensitivities are kept less than 0.1, implying that the present system remains entirely insensitive in response to perturbations of one third of the kinetic orders. On the other hand, the present system is strongly influenced by 39% of the kinetic orders. Moreover, the influences tend to appear in the rear metabolites in the cycle. This is due to the fact that, the incoming and outgoing fluxes with respect to the rear metabolites are regulated by themselves and the forward metabolites. Thus, from the standpoint of the kinetic-order sensitivities, the present model is rather sensitive and not sufficiently robust.

The results of Figs.3 suggest that the Hcy ([X.sub.4]) and AdoHcy ([X.sub.3]) concentrations are the vulnerable positions of the system that will response to any minor environmental perturbation, such as pathological consequences.

[FIGURE 3 OMITTED]

Conclusion

In the present work, we applied the computational algorithm in the formalism of BST to analyze a mathematical model developed for the methionine metabolism in human liver cell. Our method efficiently transformed the Michaelis-Menten rate laws into aggregates of power law functions or S-system equations that provided means to steady-state solution of the system. The overall steady-state sensitivity results suggest that the model is self-consistent and the steady state is locally stable. Consequently the model is mathematically guaranteed to tolerate a perturbation only if the perturbation is infinitesimally small. The analysis based on the logarithmic gains and parameter sensitivities indicate that the methionine cycle model is relatively robust.

References

[1] Savageau, M.A., 1969a, "Biochemical systems analysis I. Some mathematical Properties of the rate law for the component enzymatic reactions," J. theor. Biol., 25, pp.365-369.

[2] Savageau, M.A., 1969b, "Biochemical systems analysis II. The steady state solutions for an n-pool system using a power-law approximation," J. theor. Biol., 25, pp.370-379.

[3] Savageau, M.A., 1970, "Biochemical systems analysis III. Dynamic solutions using a power law approximation," J. theor. Biol., 26, pp.215-226.

[4] Savageau, M.A., 1971, "Concepts relating the behavior of biochemical systems to their underlying molecular properties," Arch. Biochem. Biophys., 145, pp.612-621.

[5] Savageau, M.A., 1972, "The behavior of intact biochemical control systems," Current Topics Cell. Reg., 6, pp.63-130.

[7] Finkelstein, J.D., 1990, "Methionine metabolism in mammals," J. Nutr.Biochem, 1, pp.228-237.

[8] Martinov, M.V., Vitvitsky, V.M., Mosharov, E.V., Banerjee R and Ataullakhanov F.I, 2000, "A substrate switch: A new mode of regulation in the methionine metabolic pathway," J. Theor. Biol., 204, pp.521-532.

[9] Cabrero, C., Puerta, J. and Alemany, S., 1987, "Purifications and comparisons of two forms of S-adenosyl-L-methionine synthase from rat liver," Eur. J. Biochem., 170, pp. 299- 304.

[10] Cabrero, C. and Alemani, S., 1988, "Conversion of rat liver S-adenosyl-L methionine synthase from high-Mr form to low-Mr form by LiBr," Biochim. Biophys. Acta, 952, pp.277-281.

[11] Ogawa, H. and Fujioka, M., 1982, "Purification and properties of glysine N-methyltransferase from rat liver," J. Biol. Chem., 257, pp. 3447-3452.

[12] Reed, M. C., Nijhout, H.F., Sparks, R., Ulrich, C.M., 2004, "A mathematical model of methionine cycle," J.Theor. Biol., 226, pp. 33-43.

[13] Finkelstein, J.D. and Martin, J.J., 1984, "Inactivation of betaine-homocysteine methyltransferase by adenosylmethionine and adenosylhomocysteine," Biochem. Biophys. Res. Commun, 118, pp.14-19.

[14] Shiraishi, F. and Savageau, M.A., 1992a, "The tricarboxylic acid cycle in Dictyostelium discoideum: II. Evaluation of model consistency and robustness," J. Biol. Chem. 267, pp. 22919-22925.

[15] Shiraishi, F. and Savageau, M.A., 1992b, "The tricarboxylic acid cycle in Dictyostelium discoideum: III. Analysis of steady-state and dynamic behavior," J. Biol. Chem, 267, pp.22926-22933.

(1) Jarin Akhter and (2) Fumihide Shiraishi

(1) Institute of Forestry and Environmental Sciences, University of Chittagong, Chittagong 4331, Bangladesh E-mail: jarinakhter@yahoo.com

(2) Dept. of Bioscience and Biotechnology, Graduate School of Bioresource and Bioenvironmental Sciences, Kyushu University, 6-10-1, Higashi-ku, Hakozaki, Fukuoka 812-8581, Japan
```Table 2: Metabolite concentrations, enzymatic activities and the
constants employed in the model. The S-system symbols [X.sub.1] to
[X.sub.4] represent the dependent variables and [X.sub.5] to
[X.sub.14] represent the independent variables in the methionine
cycle model.

Variable                    Symbol       Unit                  Value

[Met]                       [X.sub.1]    [mu]M              53.5
[Hcy]                       [X.sub.4]    [mu]M              0.88
[[upsilon].sub.Metin]       [X.sub.5]    [mu]M [h.sup.-1]   200.0
[V.sup.MATI.sub.max]        [X.sub.6]    [mu]M [h.sup.-1]   561
[K.sup.MATI.sub.i]                       [mu]M              41
[K.sup.MATI.sub.m]                       [mu]M              50
[V.sup.MATIII.sub.max]      [X.sub.7]    [mu]M [h.sup.-1]   22870
[K.sup.MATIII.sub.m1]                    [mu]M              19100
[K.sup.MATIII.sub.m2]                    [mu]M              21.1
[V.sup.GNMT.sub.max]        [X.sub.8]    [mu]M [h.sup.-1]   10600
[K.sup.GNMT.sub.max]                     [mu]M              4500
[K.sup.GNMT.sub.m]                       [mu]M              20
[V.sup.METH.sub.max]        [X.sub.9]    [mu]M [h.sup.-1]   4521
[K.sup.METH.sub.m1]                      [mu]M              17
[K.sup.METH.sub.m2]/[A] *                [mu]M              10
[[omega].sub.1]             [X.sub.10]   [h.sup.-1]         100
[[omega].sub.2]             [X.sub.11]    --                10
[V.sup.BHMT.sub.max]        [X.sub.12]   [mu]M [h.sup.-1]   2500
[V.sup.BHMT.sub.m]                       [mu]M              12
[V.sup.MS.sub.max]          [X.sub.13]   [mu]M [h.sup.-1]   500
[K.sup.MS.sub.m,5Hcy]                    [mu]M              0.1
[K.sup.MS.sub.m,5mTHF]                   [mu]M              25

[K.sup.MS.sub.d]                                               1
[5mTHF]                     [X.sub.14]   [mu]M              5.2
[[phi].sub.1]               [X.sub.15]   [mu][M.sup.-1]     1.7
[h.sup.-1]
[[phi].sub.2]               [X.sub.16]   [h.sup.-l]            30

* The symbol [A] represents the concentration of the substrate for
methylation.

concentrations of the dependent variables of the model compared with
the calculated steady state values of Reed et al., 2004 and
experimental values given in Martinov et al., 2000.

Variable    Calculated   Calculated *1   Observed *2       flux

[X.sub.1]    53.48393        53.5           45-80      [V.sub.+1]=
[X.sub.2]    137.9079        137.6         20-100      [V.sub.+2]=
[X.sub.3]    12.84791        13.2           1-30       [V.sub.+3]=
[X.sub.4]    0.883841        0.88           0.1-1      [V.sub.+4]=
400.9

*1) [12]
*2) [8]
```