Figure 1 shows the biochemical pathways in the hepatic cellular model used in this paper. Rectangular boxes represent the substrates that can vary in the model, and the ellipses contain the acronyms of the enzymes that catalyze particular reactions. There is one differential equation for each substrate that says that the rate of change of the concentration of the substrate is the sum of the reaction velocities (μM/hr) that produce it minus the sum of the reaction velocities that use it. Full names for all the enzymes and substrates are given in Additional File 1. Non-boxed substrates are taken to be constant or indicate products of reactions. This model extends the model for one carbon metabolism in [20] by adding the transsulfuration pathway, the synthesis of glutathione, and its transport into the blood. Here we discuss the main ideas involved in modeling the transsulfuration pathway, referring the reader to [20] (and its online supplementary material) for a discussion of the other parts of the model. A complete description of the full mathematical model and the values of all parameters are given in the online Additional File 1.

One of the most interesting features of one carbon metabolism is that reaction velocities are often affected not just by substrate and product concentrations and enzyme activities but also by the concentrations of substrates in distant parts of the reaction network that act as allosteric activators or inhibitors of the enzyme. As a result, the formulas for reaction velocities are often complicated functions of many variables. This is well illustrated by the velocity equation for the CBS reaction, the first step in the transsulfuration pathway:

$\begin{array}{l}{V}_{CBS}=\left(\frac{{V}_{max}[Hcy][Ser]}{({K}_{m}^{hcy}+[Hcy])({K}_{m}^{ser}+[Ser])}\right)\hfill \\ \cdot \left(\frac{(1.2)C{([SAM]+[SAH])}^{2}}{{30}^{2}+{([SAM]+[SAH])}^{2}}\right)\left(\frac{{k}_{a}+[{H}_{2}{O}_{2}]}{{k}_{a}+[{H}_{2}{O}_{2}]{}_{norm}}\right)\hfill \end{array}$

The first factor is simply Michaelis-Menten kinetics for the CBS reaction that uses homocysteine (Hcy) and serine (Ser) as substrates. We take the Michaelis constants from the literature. The second term is the activation of CBS by SAM and SAH that was discovered by Finkelstein and Martin [22]. The form of the activation was derived by nonlinear regression on the data in [23, 24]. The constant C is chosen so that the second term equals one in the normal steady state (see below). The last term is the activation of CBS by oxidative stress [25, 26], which is represented by the concentration of H_{2}O_{2.}

The differential equation for the concentration of cytosolic cysteine (Cys) is straightforward:

$\begin{array}{l}\frac{d}{dt}[Cys]=-{V}_{GCS}([Cys],[Glu],[GSH],[{H}_{2}{O}_{2}])\hfill \\ +{V}_{CGTL}([Cysta])+{V}_{cysin}([bCys])-\left(\frac{(0.35){[Cys]}^{2}}{200}\right)\hfill \end{array}$

The first term on the right is the production of cysteine from cystathionine (Cysta) by CGTL and the second term is the import of cysteine into the cell, which depends on the concentration of cysteine in the blood ([bCys]). The third term is the loss of cysteine in the reaction catalyzed by GCS that makes glutamyl-cysteine and the fourth term represents the loss of cysteine to other pathways (for example to sulfate and taurine). Cysteine also used for protein synthesis and is produced by protein catabolism; in the model we assume that these two rates balance. The form for the fourth term was chosen because the data in [27] indicate that at normal cysteine concentrations (approximately 200 μM) most of the flux away from cysteine is toward GSH and only a moderate amount towards other pathways. However, as [Cys] rises an increasing fraction is sent towards the synthesis of taurine. Formulas for V_{CGTL} and V_{bCYSc} (the transport of cysteine into the cell from the blood) appear in Additional File 1.

The first step in the synthesis of GSH is the formation of γ-glutamyl-cysteine (GluCys) from the constituent amino acids glutamate and cysteine by the enzyme γ-glutamyl-cysteine synthetase (GCS, also called glutamate cysteine ligase (GCL)). The reaction is reversible and GSH is a competitive inhibitor of GCS against glutamate [28–30].

The third factor in the following formula is the activation of GCS by H

_{2}O

_{2} [

25,

26].

$\begin{array}{l}{V}_{GCS}=\left(\frac{{k}_{a}+[{H}_{2}{O}_{2}]}{{k}_{a}+{[{H}_{2}{O}_{2}]}_{ss}}\right)\cdot \hfill \\ \left(\frac{{V}_{max}([Glu][Cys]-[GluCys]/{K}_{e})}{{K}_{m}^{cys}{K}_{m}^{glu}+{K}_{m}^{cys}[Glu]+{K}_{m}^{glu}[Gys](1+{\scriptscriptstyle \frac{[GSH]}{{K}_{i}}}+{\scriptscriptstyle \frac{[Glu]}{{K}_{m}^{glu}}})+{\scriptscriptstyle \frac{[GluCys]}{{K}_{p}}}+{\scriptscriptstyle \frac{[GSH]}{{K}_{i}}})}\right).\hfill \end{array}$

The second step in the synthesis of GSH is the addition of the glycine residue to GluCys by the enzyme glutathione synthase (GS). We follow [

30,

31] and use a reversible bi-reactant Michaelis-Menten mechanism.

${V}_{GS}=\left(\frac{{V}_{max}([GluCys][Gly]-[GSH]/{K}_{e})}{{K}_{m}^{glucys}{K}_{m}^{gly}+{K}_{m}^{gly}[GluCys]+{K}_{m}^{glucys}[Gly](1+{\scriptscriptstyle \frac{[GluCys]}{{K}_{m}^{glucys}}})+{\scriptscriptstyle \frac{[GSH]}{{K}_{p}}})}\right).$

The differential equation for cytosolic GSH is:

$\begin{array}{c}\frac{d}{dt}[GSH]={V}_{GS}([Gly],[GluCys]-{V}_{gsh\_out}([GSH])\\ -2{V}_{GPX}([GSH],[{H}_{2}{O}_{2}])\\ +2{V}_{GR}([GSSG],[NADPH])-(0.002)[GSH].\end{array}$

The first term is the synthesis of GSH from glycine and glutamyl-cysteine. The second term is the transport of GSH out of the liver cell into the blood. V_{gsh_out} is actually the sum of two terms, one for the high affinity transporter and one for the low affinity transporter [32]. We ignore the canalicular transport into the bile because it is a relatively small percentage of total export [33]. The third term is the rate of production of oxidized GSSG from GSH via the enzyme GPX and the fourth term is the conversion of GSSG back to GSH via the enzyme GR. We ignore other reactive oxygen species besides H_{2}O_{2}, or, put a different way, H_{2}O_{2} represents them all. The number two occurs in both these terms because two molecules of GSH combine to make one GSSG. In the fifth term we are assuming that 0.2% of the GSH is removed each hour in detoxification reactions that form conjugates [12].

The kinetics of sinusoidal efflux of GSH has been well studied in the perfused rat liver. The major part of the flux is carried by the low affinity transporter, which has sigmoidal kinetics with a V_{max} in the range 900–1400 μM/hr, a K_{m} of approximately 3000 μM, and a Hill coefficient of approximately 3 [34–36]. In our model we use a V_{max} of 1100 μM/hr, a K_{m} of 3000 μM, and a Hill coefficient of 3. We use standard Michaelis-Menten kinetics for the high affinity GSH transporter and for the two GSSG transporters.

We track five variables in the blood, [GSH], [GSSG], [Gly], [Cys] and [Glut]. Glycine, glutamate, and cysteine enter blood from intestinal absorption at rates that we vary in various experiments with the model; the normal rates are 630 μM/hr, 273 μM/hr and 70 μM/hr, respectively. We assume for convenience that the volume of the blood is the same as the volume of the liver. Glycine, cysteine, and glutamate leave the blood by transport into liver cells (depending on their concentrations) and they are also formed in the blood by the breakdown of GSH and GSSG into their component amino acids [37, 38]. We also assume that normally 10% of the cysteine, glycine, and glutamate, in the blood is taken up per hour by other cells and that an additional 25% of cysteine is converted to cystine. Under normal conditions a large percentage of blood GSH and GSSG is broken down into the component amino acids and a small amount is taken up by other cells or otherwise leaves the system. As above, full details and formulas appear in Addition File 1.

For each *in silico* computation, the values of various constants (like H_{2}O_{2}) are given, as are the methionine and serine levels in the blood, and the rates of input of cysteine glutamate, and glycine into the blood. These are the "inputs" to the model. The differential equations are then solved to determine the steady-state values of the concentrations of all the variables and the steady state rates of all the reactions. Of course, if the inputs are different the steady state will be different. We experiment with the model by changing the inputs or changing parameters (for example, a parameter that gives the strength of a particular allosteric interaction) and determine what the effect is. By removing interactions we can take the model apart piece by piece so that we can understand how and why glutathione metabolism works the way it does. We also allow the inputs to vary as functions of time (for example the amino acid input will vary because of meals) and compute the time course of each concentration and reaction rate. This allows us to investigate the homeostatic mechanisms that protect the system against fluctuations in the inputs.

A number of substrate concentrations are fixed in the model and in all the simulations reported below. These include: cytosolic GAR (10), NADPH (50), betaine (50), formaldehyde (500), dUMP (20), and total cellular folate (20). All concentrations are in μM.

### Limitations of the model

This model was designed to allow us to study various regulatory mechanisms in the transsulfuration pathway and the effects of oxidative stress, particularly as applied to Down syndrome and autism. No mathematical model can track all of the variables that might affect a complex biochemical system such as glutathione metabolism. This is also true, of course, in biological experimentation. This model is no exception. We ignore canalicular excretion of GSH. We use K_{m} values in the ranges determined experimentally but there is much less information on V_{max} values. Often we choose V_{max} values so that the steady state concentrations of substrates and products lie within the normal published ranges. Cellular amino acid concentrations are increased by feeding and protein degradation and decreased by protein synthesis, growth and use in one-carbon metabolism. In this model we assume that protein synthesis and degradation are in balance and that no amino acids are used for growth. The consequences of this assumption are outlined in the discussion.

One-carbon metabolism and the transsulfuration pathway contain many allosteric interactions by which substrates in one part of the pathway affect the activity of distant enzymes. We use experimentally determined forms for these allosteric interactions but sometimes the details of the kinetics are not known, forcing us to make reasonable educated guesses. Similarly, many effects of oxidative stress on the enzymes of one carbon metabolism and the transsulfuration pathways are known but detailed kinetics are not available.

In this paper we are mainly interested in intracellular liver metabolism, so we take a somewhat simple view of the fates glutathione and its metabolites in the blood. Future work will include a more detailed model of the blood compartment and inter-organ regulation of glutathione and its component amino acids. Thus, we do not expect that our model will make perfect quantitative predictions. Rather, we want to use it to investigate the qualitative features of glutathione metabolism in the normal state and in various disease states.