Folate system correlations in DNA microarray data

Mathematical model of Morrison and Allegra

The folate cycle model of Morrison and Allegra [1] has the following mathematical form:

These equations restate the system configuration information of Figure 1A, i.e. they state that the rate at which a metabolite concentration increases equals the sum of the synthesis reaction fluxes (arrows into a node) minus the sum of the degradation reaction fluxes (arrows leaving a node). The r _i in these equations are:

where DHF^T-DHF (total DHF minus free DHF) is the concentration of DHF bound to DHFR, x_i is the concentration of i-glutamated MTX, and all folates are assumed to be penta-glutamated. Not shown are 10 additional differential equations for up to penta-glutamation of MTX either free or bound to DHFR. These 10 equations are irrelevant when MTX = 0 as in the microarray data analyses below; they were used here only to validate the current implementation of the model against its previously published responses, see Figure 2[1].

Model limitations

Since the model has only one compartment, the cytosol, it cannot handle changes in the mitochondrial enzymes MTHFD2 and SHMT2, nor can it handle changes in the extra-cellular folate hydrolase gene FOLH1 (the gene that codes for prostate-specific membrane antigen, PSMA), so these genes were ignored in the analyses. Further, the folate genes GGH (polyglutamate hydrolase), FPGS (polyglutamate synthase), and RFC (reduced folate carrier), were not considered in the microarray analyses because these reactions are included in the model only for MTX and not folate, and because MTX = 0 for the microarray data.

Model modifications

Flux boundary conditions for dUMP and GAR synthesis in the original model [1] were replaced by downstream concentration boundary conditions set to their initial values. This was done because steady state flux differences across patients would otherwise be nullified; e.g. if the flux into GAR were fixed, the steady state flux through GART would also be fixed, and artificially then, there would be no variation in predictions of this flux across patients.

Model-data linkage

For steady state flux predictions of leukemia patient diagnostic samples, MAS5 microarray measurements of Ross et al [6] and Yeoh et al [7] were normalized by dividing by the mean of the leukemia subtype medians. Step functions from 1 (for t < 0) to the resulting metrics (for t ≥ 0) were then used as modulators of the baseline folate model V_max values, i.e. microarray data normalization values were equated to the steady state of the model. Individualized patient steady states were then computed as simulation endpoints 40 hours after the V_max perturbations; all time courses were inspected visually to assure settled steady states at 40 hours. For radiation response data [8], initial measurements were equated to the steady state of the model, so the data was normalized by its values at t = 0. Linear interpolations of the normalized data were then used as time-varying V_max modulators. For both the steady state leukemia data and the time course radiation response data, proportionality between mRNA levels and protein levels was assumed. This assumption was motivated by simplicity and a lack of better alternatives. As proteomic-transcriptomic combined dynamic response data (e.g. [9]) accrues, it will likely be replaced by a set of gene specific lead-lag filter [10] assumptions. In the meantime, it can be viewed and used as a first order approximation, or as a temporary crutch.

Computational details

The computational environment R [11] was used with R packages of Bioconductor [12] to implement this study. Specifically, the package Biobase was used to manipulate microarray data as expression set (class eset) objects and SBMLR [13, 14] was used to simulate a systems biology markup language (SBML) [15–17] representation of the folate metabolism model [1]. R scripts reproducing figures 2, 3, 4, 5, 6, 7, 8 are included with SBMLR as an illustrative example of package use. For convenience, array data used in this study have been repackaged as eset objects in R data packages available from the author's website [14]. Throughout, genes with multiple probe sets were represented by the set with the highest average value.

Folate system correlations across childhood leukemias

Childhood acute lymphoblastic leukemia (ALL) microarray data of Ross et al [6] and Yeoh et al [7] is shown in Figures 3 and 4. Several points can be made regarding this steady state diagnostic bone marrow data. Firstly, since TYMS and DHFR (similarly MTHFD1, GART and ATIC) operate in series, it makes sense that the system would attempt to match these throughput capabilities as closely as possible to avoid the costs of maintaining unneeded excess "equipment." Thus, positive correlations within the DNPS and DNTS branches are expected. Secondly, growing cells require commensurate increases in both DNTSand DNPS, so positive correlations between these cycles are also expected. Finally, DNPS genes are higher in T cell leukemic cells than in B-cell leukemic cells, consistent with measured DNPS fluxes being three fold higher in T cells than in B cells [18].

To assess the added value of the folate model, since MTHFD1 and TYMS are the gatekeepers of DNPS and DNTS, respectively, correlation plots for these genes were compared to corresponding flux predictions in Figure 5. The plots show that, for the most part, model predicted fluxes are more correlated than measured MTHFD1 vs. TYMS mRNA. To estimate the amount of correlation attributable to steady state flux constraints alone, 1000 uncorrelated normally distributed (μ = 1; σ = 0.30) random numbers were applied to the model as V_max modulators. The amount of correlation (r = 0.18, Figure 6) is more than that induced by the model (beyond gatekeeper correlations) for BCR-ABL and T cell leukemias (Figure 5), but less than induced for TEL-AML1 leukemias. Thus, for TEL-AML leukemias additional information must have been contributed by the other folate genes inputted into the model, suggesting more folate system coordination in this more curable leukemia subtype.

Folate system analysis of radiation time course data

Folate system correlations in radiation response time course data [8] were also investigated. The data (Fig. 7) shows that TS, DHFR, GARFT and MTHFD each have a dose-dependent transient increase after irradiation, consistent with radiation induced DNA damage causing a transient rise and fall in P53 activity [19] with subsequent induction of ribonucleotide reductase subunit P53R2 [20] causing a transient increase in de novo deoxynucleotide synthesis for DNA repair; a 17-fold increase in R2 protein 24 h after irradiation [21] is radioprotective [22], further supporting this conjecture. The data also shows a steady decline in many of the gene expression time courses, possibly due to radiation induced cell cycle arrest. At 10 gray, inspection of TS, DHFR, GARFT and MTHFD further suggests that RNA repair driven DNPS (peak at 2 hours) precedes DNA repair driven DNTS (peak at 6 hours).

In Figure 7, gene expression time courses are more likely to be signals if they differ between 3 and 10 gray only in terms of minor time shifts and dose ordered amplitude changes. Based on this, MTHFR and SHMT1 were dismissed as noise. The remaining six gene expression time courses were normalized by their values at t = 0 and applied to the folate model as linearly interpolated time-varying modulators of corresponding V _max parameters. The resulting model-based predicted time courses are shown in Figure 8. These plots affirm the "dose-dependent spike resting on a decreasing backdrop" response of DNPS and DNTS that was qualitatively inferred from Figure 7. In contrast to the data itself, however, at 10 gray, RNA repair driven DNPS (peak at 12 hours) and DNA repair driven DNTS (peak at 10 hours) are reverse ordered and delayed. Since NTP production tends to be in high gear full time, compared to dNTP production, it likely has a greater ability to respond rapidly to oxidative damage. Further, a dNTP synthesis flux peak at 6 hours compared to 10 hours is more consistent with a P53 spike at 5 hours [19]. Thus, until measurements of DNPS and DNTS suggest otherwise, the data-based inference that RNA repair precedes DNA repair is tentatively more credible than its model-based counterpart.