Mathematical Modeling of Glutathione Status in Type 2 Diabetics with Vitamin B12 Deficiency

Deficiencies in vitamin B12 and glutathione (GSH) are associated with a number of diseases including type 2 diabetes mellitus. We tested newly diagnosed Indian diabetic patients for correlation between their vitamin B12 and GSH, and found it to be weak. Here we seek to examine the theoretical dependence of GSH on vitamin B12 with a mathematical model of 1-carbon metabolism due to Reed and co-workers. We study the methionine cycle of the Reed-Nijhout model by developing a simple “stylized model” that captures its essential topology and whose kinetics are analytically tractable. The analysis shows—somewhat counter-intuitively—that the flux responsible for the homeostasis of homocysteine is, in fact, peripheral to the methionine cycle. Elevation of homocysteine arises from reduced activity of methionine synthase, a vitamin B12-dependent enzyme, however, this does not increase GSH biosynthesis. The model suggests that the lack of vitamin B12–GSH correlation is explained by suppression of activity in the trans-sulfuration pathway that limits the synthesis of cysteine and GSH from homocysteine. We hypothesize this “cysteine-block” is an essential consequence of vitamin B12 deficiency. It can be clinically relevant to appreciate that these secondary effects of vitamin B12 deficiency could be central to its pathophysiology.


The Reed-Nijhout model of 1-carbon metabolism
We use the computational model of 1-carbon metabolism and glutathione synthesis due to Reed et al. (Reed et al., 2008). The version of the Reed-Nijhout model we used was encoded by Lukas Endler (available for download on the European Bioinformatics Institute's models database (Endler, 2008)).
The Reed-Nijhout-Endler model contains a feeding rhythm: parameters named breakfast, lunch, dinner and fasting were used to describe the relative levels of amino acids in the blood by using a time-step function named aa input. We removed the daily feeding rhythm when using this model for our computations.
The equations and parameters of the Reed-Nijhout model ( Fig. 1) we work with in this paper are described in the Supplementary Material of (Reed et al., 2008).

Construction of the reduced methionine cycle model
While we would ideally like to study the methionine cycle by itself, it is not immediately clear how to decouple the methionine cycle from the full network. Hence we show that a reduced methionine cycle model can be constructed by carefully excising the methionine cycle sub-network from the comprehensive Reed-Nijhout model. Fig. 2a identifies the methionine cycle sub-network -the metabolites Hcy, Met, SAM and SAH -that we are interested in isolating. The corresponding equations from the Reed-Nijhout model are reproduced here for completeness: where J M S stands for the reaction velocity corresponding to the enzyme MS. The various reaction velocities are defined as: Parameter values (Reed et al., 2008) are listed in Table 1.
We begin by noting that the "rest" of the network couples to the methionine cycle through the following variables: blood methionine, betaine (BET), glycine (Gly), serine (Ser), glutathione disulfide (GSSG) and 5-methyltetrahydrofolate (c5mf). That is, the kinetics of the metabolites in the methionine cycle is linked to the full network through these six quantities.
However, when we examined these "external variables" in the Reed-Nijhout model we noted that BET, Gly, Ser and GSSG changed very little with V M S max , Fig. 7a. This suggests that a good approximation is to treat BET, Gly, Ser and GSSG as roughly constant, i.e. as parameters in the reduced model. We build the reduced model by setting these variables to the constants shown in Table 2. On the other hand, c5mf varies substantially with V M S max (Fig. 7b). This is because MS, apart from playing a role in the remethylation of Hcy to Met, also acts as an enzyme in the conversion of 5-methyltetrahydrofolate (c5mf) to THF. Fig. 7b shows that the steady state value of c5mf falls with increased activity of MS. We found a good phenomenological fit of the variation of c5mf with V M S max to be: as can be seen in Fig. 7b. The Fig. 8 above shows that the reduced model is a reasonably good approximation of Reed-Nijhout model in modeling the function of the methionine cycle. Fig. 9 gives a numerical perspective of the information shown in Fig. 8. The maximum deviation from the Reed-Nijhout model is only by about 14%, further strengthening the validity of our reduced model. These approximations allow us to decouple the methionine cycle from the full network and construct a stand-alone reduced methionine cycle model. We verified that the reduced model displays steady-state dynamics that closely approximate the full model dynamics, Fig. 8.

Construction of the stylized methionine cycle motif
To analyze the reduced methionine cycle further, we employed the following strategy: We built a simplified methionine cycle motif (Fig. 2b) -an even simpler, stylized model that reflects the topology of the reduced model.
The complex equations of enzyme-mediated fluxes of the reduced model (Fig. 2a) were replaced with simple mass action kinetics in the stylized model (Fig. 2b). For example, the constant flux of methionine from blood into the cytosol is shown as a flux with rate k 1 in the stylized model. The methionine flux from the cytosol into the blood, regulated by the rate constant k outmet met , is assumed to be regulated by the rate constant k −1 in the stylized model. Fluxes due to enzymes MAT1 and MAT3, that regulate the conversion of Met into SAM, are clubbed together and simplified to be shown as a single flux regulated by the rate constant k 3 in the stylized version. Fluxes due to enzymes GNMT and DNMT, that regulate the conversion of SAM into SAH, are also clubbed together and simplified to be shown as a single flux regulated by the rate constant k 4 in the stylized version. SAHH, which regulates the interconversion between Hcy and SAH in the reduced model, is replaced by the rate constants: k 5 to regulate the conversion from SAH to Hcy, and k −5 to regulate the backward conversion from Hcy to SAH in the stylized model. The rate constant k 2 in the stylized model replaces the rate constant regulating the conversion of Hcy into the cystathionine residues. Fluxes due to enzymes MS and BHMT, that regulate the remethylation of Hcy into Met, are clubbed together and simplified to be shown as a single flux regulated by the rate constant k 0 in the stylized model.

A weak methionine efflux explains the homocysteine homeostasis
We used the reduced methionine cycle model and the full Reed-Nijhout model to test the hypothesis that the sensitivity of Hcy to changes in V M S max is dependent on the strength of the methionine efflux via k outmet met . In Fig. 10 we show evidence that the models that support the hypothesis and thus, provide an explanation to the homeostatic behavior of Hcy.
We varied the value of k outmet met in the reduced model over three orders of magnitude and observed the changes in Hcy with respect to V M S max (Fig. 10a). At lower values of k outmet met , Hcy is almost unaffected by changes in V M S max whereas for higher values of k outmet met , Hcy varies significantly with V M S max . In vitamin B 12 deficient scenario (V M S max < 500 µM/Hr), a larger value of k outmet met can be seen to lead to a higher steady state of Hcy resulting in hyperhomocysteinemialike conditions.