A Mathematical Model of the Folate Cycle

A mathematical model is developed for the folate cycle based on standard biochemical kinetics. We use the model to provide new insights into several different mechanisms of folate homeostasis. The model reproduces the known pool sizes of folate substrates and the fluxes through each of the loops of the folate cycle and has the qualitative behavior observed in a variety of experimental studies. Vitamin B12 deficiency, modeled as a reduction in the Vmax of the methionine synthase reaction, results in a secondary folate deficiency via the accumulation of folate as 5-methyltetrahydrofolate (the “methyl trap”). One form of homeostasis is revealed by the fact that a 100-fold up-regulation of thymidylate synthase and dihydrofolate reductase (known to occur at the G1/S transition) dramatically increases pyrimidine production without affecting the other reactions of the folate cycle. The model also predicts that an almost total inhibition of dihydrofolate reductase is required to significantly inhibit the thymidylate synthase reaction, consistent with experimental and clinical studies on the effects of methotrexate. Sensitivity to variation in enzymatic parameters tends to be local in the cycle and inversely proportional to the number of reactions that interconvert two folate substrates. Another form of homeostasis is a consequence of the nonenzymatic binding of folate substrates to folate enzymes. Without folate binding, the velocities of the reactions decrease approximately linearly as total folate is decreased. In the presence of folate binding and allosteric inhibition, the velocities show a remarkable constancy as total folate is decreased.

The folate cycle plays a central role in cell metabolism. Among its important functions are the delivery of one-carbon units to the methionine cycle, for use in methylation reactions, and the synthesis of pyrimidines and purines. Dietary folate deficiency or genetic polymorphisms in folate-metabolizing enzymes are associated with megaloblastic anemia, developmental abnormalities including neural tube defects and Down's syndrome, and various types of cancer, especially those of the gastrointestinal tract and leukemias (1)(2)(3)(4)(5)(6)(7)(8). Elevated homocysteine concentrations, a biomarker of a low folate status, have been implicated in cardiovascular diseases and Alzheimer's disease (6,9,10). Further-more, several enzymes in the cycle are the targets of cancer chemotherapeutic agents (11,12). Because of its importance in human health, folate metabolism has long been the focus of clinical, nutritional, and biochemical investigations. Biochemical studies have provided extensive and detailed information about each of the enzymes and metabolites of the folate cycle. As a consequence, the components and the reaction diagram of the folate cycle are well understood (13).
The structure of the folate cycle is relatively complex (Fig. 1) and consists of several interacting loops. Most of the reactions depend in a nonlinear way on the concentrations of their substrates. Therefore, the behavior of the full cycle in response to dietary and genetic variation cannot be easily deduced from the reaction diagram itself. The global behavior of the cycle can, however, be investigated through a mathematical model of the reactions of and interactions among the components of the cycle.
Our first purpose in this paper is to develop a basic mathematical model of the cytosolic folate cycle based on known biochemistry. The model uses standard reaction kinetics and does not, at this stage, incorporate details such as polyglutamation, compartmentalization, tissue specificity, and long range inhibitions. Despite these simplifications, this basic model reproduces and explains many experimental findings on the global behavior of the system, including system homeostasis in the presence of chemotherapeutic agents and up-regulation of enzymes during the cell cycle. Our second purpose is to investigate the role of folate-binding proteins in the control of folate homeostasis. We show that allosteric inhibition by folate substrates results in a remarkable robustness of reaction velocities as total folate is reduced.
We view this basic model as a platform on which we can build tissue-specific models incorporating additional details such as those mentioned above. Such more detailed models can then be used to investigate specific biological and biochemical questions, such as the roles of polyglutamation, the control of purine and pyrimidine synthesis, and the effects of chemotherapeutic agents (11). In addition, the basic model, when combined with our previous model of the methionine cycle (14), can be used to address questions involving the interactions between these two cycles.
In a series of papers, Jackson (15)(16)(17)(18)(19)(20) showed that mathematical models can be used to develop a better understanding of folate metabolism. His studies inspired others to develop mathematical models to understand the effects of anti-folate drugs (21)(22)(23)(24). These studies are now more than 15 years old and, naturally, were unable to utilize the wealth of recent biochemical and genetic findings about folate metabolism. The present model builds on these previous efforts but incorporates recent biochemical findings to build a more comprehensive model of the folate cycle. In addition, our purposes are more general in that we are not primarily concerned with the effects of antifolate drugs but instead wish to understand the emergent properties of this nonlinear complex interconnected system and the responses of the folate cycle and methionine cycle to dietary and genetic variation. Fig. 1 illustrates the components of the folate cycle treated in this paper. The structure of the cycle and the relevant reactions are based on the reaction scheme of Wagner (13) and represent the cytosolic reactions in mammalian liver.

Structure of the Model-
Most of the reactions in the folate cycle are bimolecular. In the present study, we focus on variation of the folate compounds and therefore treat other substrates as constants. We assume that each of the reactions is a random bimolecular reaction. For the unidirectional reactions, the velocity can be written as follows (25), where F indicates the folate substrate, and S indicates the nonfolate substrate. We take ␥ ϭ 1, which means that each substrate does not affect the binding of the other. Thus, the velocity of each unidirectional reaction can be written in the following simple form.
In the case of reversible reactions we take the analogous form, where S f , F f and S r , F r refer to the substrates of the forward and reverse reactions, respectively. The forward direction is indicated for each reversible reaction in Table II. There is a nonenzymatic interconversion between THF 1 and 5,10-CH 2 -THF that we assume to be pseudo-firstorder mass action in each direction (15,23). We shall denote the rate of each of the reactions in Fig. 1 by the symbol V with a subscript indicating the enzyme that catalyzes the reaction; thus, V MS denotes the rate of the reaction catalyzed by methionine synthase. The rate for the nonenzymatic interconversion between THF and 5,10-CH 2 -THF is denoted by V NE . Each of these V values varies in time, of course, since the current rate depends on current substrate concentrations. Thus, both the basic folate model and the model with enzyme inhibition have the following form.
The basic folate model and the enzyme-binding folate model differ in the form of the V values and the choice of some rate constants.
Parameter Values for the Basic Model- Table II shows that for most of the reactions in Fig. 1 there are a wide range of K m and especially V max values reported in the literature. The great variation in parameter values is probably due to the fact that different authors use different species, different tissues, and different cell lines (cancerous versus noncancerous). Moreover, the V max of some enzymes varies with the cell cycle. Hence, the literature does not give strong guidance on parameter choices for a general model of folate metabolism. Our initial choices of parameter values for the basic model under normal cytosolic physiological conditions were based on the following considerations. First, in our earlier model of the methionine cycle (14), we found the velocity of the reaction catalyzed by methionine synthase, V MS , to be ϳ85 M/h. This implies that at steady-state, V MTHFR is also ϳ85 M/h (see Fig. 1). Second, the ratios of the rates of three reactions that use 5,10-CH 2 -THF as a substrate, V MTHFR , V MTD , and V TS , were found by Green et al. (31) to be ϳ35:1100:1. If V MTHFR is ϳ85 M/h, then V MTD and V TS must be ϳ2670 and 2.42 M/h. At steady state, V DHFR ϭ V TS and V MTCH ϭ V MTD . The sum of the rates of the four reactions that interconvert 10f-THF and THF must also equal V MTD . Finally, The intracellular concentrations of the various folate co-enzymes are well established, and we are assuming constant values for the nonfolate substrate concentrations (Table I) Table  II). Therefore, we provisionally chose reasonable values for these constants within the published ranges.
The basic model does not include intracellular catabolism of folates. In the model, folate enters and leaves the cell as 5mTHF as indicated by the rates F in and F out in the differential equation for 5mTHF. Approximately 0.8% of the total body folate pool is lost by excretion per day (32)(33)(34). We assume that the various body pools of folate are in equilibrium and that therefore the intracellular folate pool is lost at the same rate. Therefore, since our steady-state folate concentration is 20.0 M, we take F in ϭ (20)(0.008)/24 ϭ 0.0067 M/h. If we assume that output from the cytosol is a first-order rate process with rate constant ␣, then at steady-state we must have F out ϭ 0.0067 ϭ ␣[5mTHF]. Since the normal steady-state concentration of 5mTHF is 5.16 M, this determines ␣ ϭ 0.0013/h.

The Basic Folate Model
We show that the basic folate model has the qualitative behavior seen in a variety of experimental studies.
Steady-state Behavior-At steady state, with a total folate pool of 20 M, the concentrations of the folate metabolites and the fluxes are summarized in Table III Tables I and II).
In all cases, including those not illustrated in Table III, the flux control coefficients were much smaller than 1. In general, they were Ͻ0.4 and were only larger than 0.5 for the fluxes of the varied enzyme and the steady-state concentration of its substrate. Perhaps the most interesting overall pattern emerging from the total data set is that the sensitivity to variation in enzymatic parameters is inversely proportional to the number of reactions that interconvert two folate substrates. The detailed analysis of the sensitivities of the folate cycle will be the subject of a separate publication.
The Methyl Trap Hypothesis-It is well known that vitamin B 12 deficiency results in a secondary folate deficiency. This observation is explained by the "methyl trap" hypothesis, which proposes that B 12 deficiency reduces the activity of MS and this leads to the accumulation of 5mTHF at the expense of other folate forms (50 -52). This reduces the rates of other reactions of folate metabolism. In addition, the monoglutama-ted form of 5mTHF leaks out of cells, thus reducing the intracellular folate pool. In our simulations, vitamin B 12 deficiency is modeled by reducing V max MS . The results of a progressive reduction in V max MS are illustrated in Fig. 2, which shows that there is a progressive accumulation of folate as 5mTHF. This is accompanied by a reduction is other folate forms and a concomitant drop in the reaction velocities throughout the folate cycle.
The Dihydrofolate Loop-In the basic model, the rate V TS is very low. Indeed, this reaction is known to be a rate-limiting step for DNA synthesis and varies greatly with the stage of the cell cycle (53-55, 57, 58). The transcription of thymidylate synthase and dihydrofolate reductase genes are co-regulated by the E2F transcription factor (59 -63) and increase greatly at the G 1 /S transition (56). It has been reported that this increase in transcription coincides with a 100-fold increase in the activity of thymidylate synthase (56). To determine the effect of this change on the rest of the folate cycle, we increased the V max values of TS and DHFR 100-fold. The results can be seen in the second column of Table III (TS,DHFR1). Although V TS and V DHFR increased almost 100-fold, the folate pools and the other fluxes changed by less than 7%. Thus, the up-regulation of TS  and DHFR greatly increases dTMP synthesis but has little effect on the rest of the folate cycle.
Suppression of dihydrofolate reductase by methotrexate and of thymidylate synthase by 5-fluorouracil are used in cancer chemotherapy to inhibit pyrimidine synthesis and therefore cell division. These drugs act by binding tightly to these enzymes, so we model their effects by reducing the corresponding V max values. Reducing V max TS to 10% of its normal value reduced V TS , V DHFR , and [DHF] proportionally but had little or no effect on the fluxes and pools elsewhere in the folate cycle. By contrast, reducing V max DHFR has very little effect on V TS , V DHFR , and [DHF] until it falls below 10% of its normal value (Fig. 3). This insensitivity of V TS to a substantial reduction of DHFR activity is in accord with observations in the literature (20) that very large doses of methotrexate are required to lower V TS although methotrexate significantly inhibits DHFR. The model can be used to explain this behavior. Under "normal" conditions, [5,10-CH2-THF] is low and [DHF] is very low (see Table III). To reduce V TS , [5,10-CH 2 -THF] must be substantially reduced, and since it is small to begin with, this only happens when much of the total folate accumulates as DHF. Since [DHF] is normally exceedingly low, this requires almost complete inhibition of DHFR.
Folate Half-life-Although the folate is removed from the system by a first-order process, we cannot calculate the half-life of total folate from the rate constant, ␣, itself. This is because not all folate is in the form 5mTHF, and thus the half-life will depend on all of the reactions of the system, which govern how the various folates are redistributed after some of it is removed.
To calculate the half-life using the model, we set F in ϭ 0 and the total initial cytosolic folate pool to 20 M. The model predicts that 89 days are required to reduce the total initial cytosolic folate pool by half. This corresponds well with total body half-life of ϳ100 days (32)(33)(34)64).
Dependence on Total Folate-The basic model is sensitive to variations in the total folate pool. The steady-state concentrations and most velocities decrease approximately linearly with decrease in the size of the folate pool (Fig. 4). The reason for the approximate linearity of the velocities is that in most cases (except V PGT and V AICART ) the steady state concentration of the folate metabolite is considerably smaller than the associated K m values (Table II).
Normal daily fluctuations of folate intake do not affect cellular folate levels very much, since, on the average, only 0.8% of total folate is consumed and excreted each day (see above). Long term folate deficiency that reduces total folate by one-half reduces the velocities of most of the reactions by approximately one-half. Thus, the basic folate model predicts that folate function should proportionally decrease with total folate, and there-fore one should expect to see functional consequences such as hyperhomocysteinemia and megaloblastic anemia within months. However, long term folate deprivation studies in humans showed a rapid decline in serum folate followed by "weeks or months when the serum folate concentrations is low but there is no other evidence of folate deficiency" (65). Thus, the basic folate model does not account for the observed homeostasis of folate function under long term folate deprivation and gradually declining total folate pools. We shall see under "Folate-Enzyme Binding and Homeostasis" that the inclusion of

Folate-Enzyme Binding and Homeostasis
It has been known for some time that in mammalian liver folates are tightly bound to a number of specific folate-binding proteins (66 -70). Interestingly, these folate-binding proteins have turned out to be the enzymes involved in the folate cycle (26,42,(71)(72)(73)(74). The total concentration of folate binding sites on these proteins exceeds the total concentration of the folate pools, and they bind folates with dissociation constants in the 100 nM range. This binding not only reduces pools of free folates but also inhibits the activities of the enzymes. The enzyme inhibitions by folate identified to date are shown in Table IV. These observations raise the question of how the velocities of the folate reactions are maintained at low free folate levels and how the inhibition of folate enzymes affects the overall functioning of the folate cycle.
Many of the inhibitions shown in Table IV are substrate inhibitions. Below we first illustrate the effect on reaction velocity when an enzyme is inhibited by the noncatalytic binding of its substrate (see "Enzyme Inhibition"). We shall see that substrate inhibition has an important homeostatic effect. Under "Folate Homeostasis," we carry over this analysis to the whole folate cycle. We shall see that the homeostatic effect of folate binding does not depend on the identity of the inhibiting folate and that the homeostatic effect becomes complex when not all enzymes are inhibited.
Enzyme Inhibition-It is well known that enzymes can be inhibited by noncatalytic binding of their substrates (86 -88). The two most straightforward mechanisms are noncompetitive and uncompetitive inhibition in which the substrate binding alters the apparent V max or the K m , respectively, of the enzyme (25,88). We shall assume that the inhibition is noncompetitive so we obtain the following expression for the reaction velocity at steady state, declines, more enzyme is released from inhibition. Folate Homeostasis-In order to see the effect of the above mechanism on the overall operation of the folate cycle, we consider two cases. In the first case, we assume that all reactions in the folate cycle are inhibited noncompetitively by THF. In the second case, we assume that THF inhibits SHMT, MTD, MTCH, FTS, and FTD and that DHF inhibits MTHFR.
In the first case, we multiply the equation for the velocity of each reaction in the basic folate model by the following factor.
Note that this situation is somewhat different from classical substrate inhibition, because THF is acting as a noncompetitive inhibitor also in reactions in which it is not a substrate. We choose A so that the factor equals 1 when [THF] has its steadystate value at a total folate pool of 20 M. Thus, when total folate equals 20 M, all of the steady-state velocities and metabolite concentrations are the same as in the basic folate model (Table III). Note, however, that in this case the "concentrations" of the metabolites include both bound and free forms. The above factor in effect reduces the velocity of each reaction by dividing the V max by 1 ϩ [THF]/K i and scales up the V max to compensate so the folate cycle runs "normally" when total folate (bound plus free) is 20 M. The question is what happens to the velocities when total folate is reduced? We found that in the presence of THF inhibition, the overall fluxes through the folate cycle were remarkably robust to variation in the total folate pool. In Fig. 5 we show the steady-state velocities of some of the principal reactions as functions of the total folate pool for the case K i ϭ 1. As can be seen, the folate pool can be reduced to less than a quarter of its normal value of 20 M without significantly affecting the velocities. The reactions not illustrated showed similar stable plateaus between 20 and 5 M. The stabilizing effect of folate binding was stronger as K i was lowered (stronger inhibition). This is the full folate cycle analogue of the mechanism discussed above for a single enzyme. As total folate decreases, free folate declines, but enzyme inhibition is relieved so as to maintain a nearly constant overall reaction velocity. This suggests that folate binding does not just act as a reservoir for folate but is also a dynamic mechanism for maintaining folate homeostasis.
The results above do not depend on our choice of THF as the universal inhibitor. Suppose, for example, we choose DHF as the universal inhibitor with its K i chosen so that the following is true,  which ensures that the relative strengths of the inhibitions are the same despite the significant difference in the concentrations of these two metabolites. Then the graphs of the velocities as total folate declines are the same as in Fig. 5 (simulation not shown). The reason is that as total folate declines, the individual folate concentrations decline proportionally. This homeostatic mechanism does not require a universal folate inhibitor but works equally well if different enzymes are inhibited by different folates. If the ratios of the folate concentrations to their K i values are as in Equation 12, the graphs of the reaction velocities as a function of the total folate pool are identical to the graphs in Fig. 5 (simulations not shown).
We also investigated the natural question of what the results would be if some but not all of the folate enzymes were inhibited by folates. Thus, in case 2, we included only the inhibition of SHMT, MTD, MTCH, FTS, and FTD by THF and the inhibition of MTHFR by DHF. We chose K i THF ϭ 1 M as above and K i DHF so that Equation 12 holds. The consequences of this pattern of inhibition are illustrated in Fig. 6. As can be seen, each velocity responded differently as total folate was lowered from 20 to 1 M.
V PGT (and also V AICART , V TS , and V DHFR ; not shown) scales down approximately linearly with total folate. This is not surprising, since these reactions were not subject to inhibition. Note, however, that V MS shows a strong homeostatic effect, although it is not inhibited. This is because V MS must equal V MTHFR at steady-state, and MTHFR is inhibited by DHF so that V MTHFR shows a strong homeostatic response.
The behavior of the inner cycle is more complicated. The strongest homeostatic effect is exhibited by SHMT, FTS, and FTD; in fact, the velocities of FTS and FTD go up as total folate is lowered from 20 to 10 M. The shapes of the two velocity curves are different (although the inhibition is the same) because the K m values are different. The homeostatic effect is relatively small for MTD and MTCH although they are just as strongly inhibited. The reasons are as follows. V PGT and V AICART scale down linearly with total folate, and this causes a relative shift of folate from the THF to the 10f-THF pool. This causes V NE , which is linear in [THF], to decrease faster than linearly as total folate is reduced (not shown). Since V SHMT changes little with decreasing folate, the total flux V SHMT ϩ V NE declines as the total folate pool is reduced. At steady state, the fluxes V MTD and V MTCH must be close to V SHMT ϩ V NE , because V TS and V MTHFR are relatively small. Thus, V MTD and V MTCH must also decline as total folate decreases.
In case one, where all of the reactions were equally inhibited, the velocities responded in similar ways to variation in the folate pool (Fig. 5). By contrast, in case two, where only certain reactions were inhibited, a much more complex but understandable pattern of homeostasis resulted. This complexity is not surprising in such a highly interconnected and nonlinear system. DISCUSSION We have constructed a basic model for the folate cycle that incorporates most of the known reactions. This model is able to reproduce many of the known properties of folate metabolism such as the methyl trap, the relative insensitivity to methotrexate inhibition of DHFR, and the long half-life of folate. In the second part of the paper, we extended the basic model to include folate binding and sequestration and examined the consequences for folate homeostasis.
The ultimate goal of our work is not just to understand the biochemistry of folate metabolism but to understand how this complex cycle performs its various biological functions in the presence of genetic and environmental variation. Folate metabolism is not a single cycle but is, in fact, a system of at least three interlocked cycles (Fig. 1). For example, in the central cycle in Fig. 1, SHMT picks up single carbon units from serine that can be used for purine synthesis by AICART, for pyrimidine synthesis in the TS-DHFR cycle, and for the remethylation of homocysteine in the MTHFR-MS cycle. An important biological question is how the necessary changes required by cell function in each of these cycles affect the behavior of the other cycles. In the one case investigated in this paper, we found that dramatic variation of the activity of TS and DHFR with the cell cycle does not affect the rest of folate metabolism (Table III). Indeed, we found that in most cases, enzyme parameter variation only affected local velocities and pool sizes.
We chose parameter values for the model to be within the (sometimes very broad) published ranges and to satisfy the known values of folate pool sizes and fluxes through various portions of the cycle. These constraints do not uniquely determine all of the parameter values, and in some cases we simply chose reasonable values within the published ranges. Our sensitivity analysis, illustrated by three examples under "Dependence on Parameters," shows that the overall qualitative behavior of the cycle is not very sensitive to the exact choices made. It may be that it is the structure of the cycle itself that causes this relative insensitivity to parameter variation and that the broad range of published values for the parameters is not only due to different choices of experimental tissues but also reflects the fact that there is little natural selection for parameter constancy within a robust and stable system.
The basic model does not include certain aspects of folate metabolism. For instance, for simplicity we have not included the reactions catalyzed by the bifunctional enzyme formiminotransferase-cyclodeaminase that converts THF to 5-formiminotetrahydrofolic acid and then to 5,10-CHϭTHF, creating a fourth cycle inside the central cycle (13,89). The overall effect of adding these reactions would be to increase the flux from 5,10-CHϭTHF to 10f-THF to THF and to decrease somewhat the fluxes emanating from 5,10-CH 2 -THF. As a second example, we have not considered the compartmentalization of folate metabolism between mitochondria and cytoplasm. Mitochondria contain a substantial portion of the cellular folate and some of the enzymes of the central cycle, namely SHMT, MTD, MTCH, FTD, and FTS (13,26). It is generally thought that mitochondria import serine from and export formate to the cytoplasm (18,38,90). Our basic model should be thought of as summing the folate cycles of cytoplasm and mitochondria. Separating the cytoplasmic cycle from the mitochondial cycle would probably lower substantially the flux through the central cycle in the cytoplasm without affecting the outer two cycles. This important aspect of folate metabolism as well as the roles of formiminotransferase-cyclodeaminase, polyglutamation, long range inhibitions, and tissue specificity will be investigated in future studies.
Although folate-binding proteins have long been known, the functional significance of folate sequestration by these proteins is still unclear. Most authors assume that folate binding primarily serves as a storage mechanism for folate substrates. Our analysis shows that folate binding may have a substantial homeostatic effect, in that it allows many reaction velocities in the cycle to be maintained as total folate decreases. The kinetics of many of the inhibitions have yet to be established, and all of the possible inhibitions have not yet been identified, so the entire scope of this homeostatic mechanism cannot yet be determined. We considered two cases. In the first case, we assumed that a single folate noncompetitively inhibits each reaction of the folate cycle, and we observed a homeostatic effect on all reactions (Fig. 5). In the second case, we assumed that only a subset of the enzymes in the folate cycle were inhibited, and we saw substantial homeostatic effects, focused largely, but not exclusively, on the inhibited reactions (Fig. 6). In other simulation experiments, we found that the homeostatic effect did not depend on the identity of the folate species involved in the inhibition. As more quantitative information about folate binding becomes available, it can be incorporated into the model, and further hypotheses can be tested.