Correcting for the effects of natural abundance in stable isotope resolved metabolomics experiments involving ultra-high resolution mass spectrometry

Background Stable isotope tracing with ultra-high resolution Fourier transform-ion cyclotron resonance-mass spectrometry (FT-ICR-MS) can provide simultaneous determination of hundreds to thousands of metabolite isotopologue species without the need for chromatographic separation. Therefore, this experimental metabolomics methodology may allow the tracing of metabolic pathways starting from stable-isotope-enriched precursors, which can improve our mechanistic understanding of cellular metabolism. However, contributions to the observed intensities arising from the stable isotope's natural abundance must be subtracted (deisotoped) from the raw isotopologue peaks before interpretation. Previously posed deisotoping problems are sidestepped due to the isotopic resolution and identification of individual isotopologue peaks. This peak resolution and identification come from the very high mass resolution and accuracy of FT-ICR-MS and present an analytically solvable deisotoping problem, even in the context of stable-isotope enrichment. Results We present both a computationally feasible analytical solution and an algorithm to this newly posed deisotoping problem, which both work with any amount of 13C or 15N stable-isotope enrichment. We demonstrate this algorithm and correct for the effects of 13C natural abundance on a set of raw isotopologue intensities for a specific phosphatidylcholine lipid metabolite derived from a 13C-tracing experiment. Conclusions Correction for the effects of 13C natural abundance on a set of raw isotopologue intensities is computationally feasible when the raw isotopologues are isotopically resolved and identified. Such correction makes qualitative interpretation of stable isotope tracing easier and is required before attempting a more rigorous quantitative interpretation of the isotopologue data. The presented implementation is very robust with increasing metabolite size. Error analysis of the algorithm will be straightforward due to low relative error from the implementation itself. Furthermore, the algorithm may serve as an independent quality control measure for a set of observed isotopologue intensities.


Background
Application of mass spectrometry to stable isotope tracing experiments for the elucidation of glucose dates back to at least the early 1980's [1,2]. The general scheme for these experiments is to supply a labeled precursor such as uniformly-labeled 13 C glucose ([U-13 C]glucose) to a bacterial culture, tissue culture, or a whole multicellular organism and then extract a set of cellular or excreted metabolites for analysis [3,4]. For identified metabolites, specific patterns of isotopologues are usually observed, which are then interpreted within the context of known cellular metabolic pathways [3][4][5]. Recently, we applied this technique to elucidate specific aspects of lipid metabolism [6].
The ultra-high resolution capability of Fourier transform-ion cyclotron resonance-mass spectrometry (FT-ICR-MS) makes it possibility to identify simultaneously hundreds, if not thousands, of metabolites from crude cell extracts without the need for chromatographic separation [6]. The better than 1 ppm mass accuracy of state-of-the-art FT-ICR-MS is often high enough to provide mass-to-charge ratios (m/z) down to the 3 rd and 4 th decimal place for metabolites less than a few thousand Daltons. This is accurate enough to distinguish relativistic mass differences between expected isotopes of CHONPS elements and unambiguously determine the isotope-specific molecular formula of an individual peak. Furthermore, the FT-ICR-MS's high mass resolution allows for the direct detection or deconvolution of individual isotopologues or mass-equivalent sets of isotopomers for a given metabolite.
Isotopologue identification and quantification of thousands of metabolites in these metabolomic experiments can provide a wealth of data for modeling the flux through metabolic networks. But before isotopologue intensity data can be properly interpreted, the contributions from isotopic natural abundance must be factored out (deisotoped). This is a computationally expensive and analytically intractable problem for data from lower mass resolution spectrometers where individual isotopically-resolved isotopologues cannot be distinguished [7]. In these instances, numerical methods have been employed to approximate and subtract the contributions from isotopic natural abundance [4,[7][8][9]. Some of these calculations are aimed at a different deisotoping problem, namely identifying the related isotopologues and calculating the monoisotopic mass from its isotopic mass distribution [10,11]. Fortuitously, with the isotoperesolved isotopologue peaks from FT-ICR-MS histograms, we can pose a similar but distinct problem that allows for the derivation of a computationally tractable analytical solution. In addition, isotopologues derived from the same molecule (or very similar set of molecules) neatly handle peak intensity referencing issues by providing a natural internal reference.

Derivation of the analytical solution
. , represents the relative distribution of carbon isotopologues from natural abundance only, as a sum of multinomial coefficients multiplied by the intensity of I M+0 , the theoretically untainted 12 C monoisotopic peak. The terms being summed are similar in form to those presented in Snider, 2007. I M+i;NA is the expected intensity of the i th isotopologue peak representing i additional nucleons. NAx C is the fractional natural abundance of the X C isotope. C Max is the number of carbons in the molecule. The multinomial coefficients, derived from the multinomial theorem with 3 variables represent the number of possible isotopomers of identical mass for a molecule with C Max carbons given 3 isotopes of carbon: 12 C, 13 C, and 14 C.
Isotopologue peaks containing 14 C are typically not observed, since the isotope is very rare. Moreover, due to the very high mass resolution in FT-ICR-MS histograms, isotopologue peaks representing molecules comprised exclusively of the major isotope of CHONPS elements (expected elements for biological systems) along with 13 C, are completely resolved/deconvoluted and identified. Thus, we can ignore the contributions from 14 C and from minor isotopes of all other elements excluding carbon. This simplifies the calculation to a single term with a binomial coefficient (binomial term) shown in Equation 2, where NA13 C ≈ 0.01109. The binomial coefficient represents the number of possible isotopomers of identical mass for a molecule with C Max carbons given only 2 isotopes of carbon: 12 C and 13 C.
At natural abundance, each peak, I M+i;NA , is directly related to the theoretically untainted 12 C monoisotopic peak, I M+0 , that has a fractional intensity of 1 when dividing by the sum of isotopologue intensities. However, once 13 C is incorporated into the molecule from a labeling source, the calculation of the contributions from natural abundance becomes more complex [8,9]. The effects of 13 C natural abundance now depend on the amount of 13 C label already present. With each 12 C/ 13 C isotopologue resolved in the FT-ICR-MS histogram, we can use a series of binomial terms to accurately describe and correct for 13 C natural abundance. Equation 3 shows the basic form of these binomial terms as B C (n, k) where k represents the total number of 13 C carbons present, n represents the number of 13 C carbons due to incorporation from a labeling source, and k-n is the number of 13 C carbons due to natural abundance. The binomial coefficient in Equation 3 enumerates the number of ways that k-n 13 C carbons can be incorporated into the molecule when n carbons are already labeled with 13 C. Equation 4 shows the first series needed in the correction, B C sum(n) which represents the fraction of I M+i intensity that is converted to other isotopologues due to the effects of natural abundance.
shows the full correction as the original isotopologue intensity minus natural abundance contributions based on lower mass untainted isotopologue intensities. Division by the fractional intensity, 1 -B C sum(i), compensates for natural abundance effects that lower the intensity of the given isotopologue. As illustrated in Table 1, Equation 5 must be applied in a sequential fashion starting with i = 0, since the results of each step are needed in subsequent steps. In other words, the natural abundance corrected intensities of isotopologues with lower 13 C incorporation from labeling are needed to calculate the natural abundance correction of isotopologues with higher 13 C incorporation from labeling.
Since 15 N incorporation can be distinguished from 13 C incorporation due to the very high mass resolution in FT-ICR-MS histograms, it takes only a trivial conversion of Equations 3, 4, and 5 to handle labeling in 14 N/ 15 N isotopologues. We simply replace N Max for C Max and NA15 N for NA13 C . However, handling all of the mixed 14 N/ 15 N/ 12 C/ 13 C isotopologues that arise from simultaneous 13 C and 15 N labeling requires a series of two binomial terms multiplied together as shown in Equations 6 and 7. Given the peaks are isotopically resolved, there are CMax * NMax separate observable isotopologues, whose intensities are represented by I M+i, j;NA . The multiplied binomial terms, B C (x, i) * B N (y, j), describe the combined effects from both carbon and nitrogen natural abundance.
A version of each equation in larger fonts is available in Additional file 1.

Implementation of the algorithm
We implemented Equations 3, 4, and 5 as an iterative algorithm in the Perl programming language [Additional file 2]. Iteration allows the algorithm to partially compensate for missing (zero intensity) isotopologues. The algorithm ( Figure 1) starts with C Max and the observed 12 C/ 13 C isotopologue intensities contaminated by contributions from 13 C natural abundance. Based on C Max , the algorithm precalculates the binomial coefficients needed in later steps using Equations 3 and 4. During each iteration, the algorithm performs three steps. In step 1, the algorithm calculates the set of uncontaminated 12 C/ 13 C isotopologue intensities using Equation 5 and the observed intensities supplemented with calculated contaminated intensities for missing isotopologues. From Equation 5, it is apparent that this must be done in ascending mass order starting with I M+0 . Sometimes, small negative uncontaminated intensities arise from errors in the observed intensities. These negative intensities are flattened to zero, since they have no basis in reality. Next, the algorithm renormalizes the uncontaminated intensities based on the sum of observed intensities. This is required since missing isotopologues were supplemented with calculated values and negative intensities are flattened to zero. In step 2, the algorithm calculates the set of contaminated intensities based on the uncontaminated set by solving for I M+i;NA in Equation 5. In step 3, the algorithm calculates the absolute difference between observed and calculated contaminated intensities. If this difference decreases, the algorithm performs another iteration until no more improvement is seen. Finally, the algorithm prints the results and ends.

Testing the implementation
We created several sets of simulated isotopologues (test sets) with varying levels of 13 C-labeling and added the expected contributions (contamination) from 13 C natural abundance by solving for I M+i;NA in Equation 5. We then tested the implementation with these test sets. Figure 2 shows the results for three of these test sets of a hypothetical metabolite with 20 carbon atoms. The 13 C Table 1 Sequential correction of 13 C natural abundance effects in a four-carbon example natural abundance contaminated intensities are in red and corrected intensities in green. The red bars in Figure 2A represents the expected observed isotopologue intensities when no 13 C-labeling is present. This naturally collapses into a single green 12 C monoisotopic peak with correction. Figure 2B shows the contaminated and corrected isotopologue intensities when equal amounts of 13 C-labeling for 8, 10, and 12 carbons are present. There is a tapering phenomenon observed in the contaminated intensities due to the fact that the number of carbons affecting the intensities decreases with increasing amounts of 13 C-labeling. Equation 3 captures this phenomenon within its binomial coefficient where it is further demonstrated in Figure 2C with natural abundance having no effect on a metabolite with 100% 13 C-labeling.
The implementation is also quite efficient even in an interpreted programming language like Perl. 10,000 repetitions of this algorithm for all 3 simulated test sets took only 17 seconds on a single core of an Intel T7200 Core 2 Duo mobile processor with 2GB of RAM and running release 5.3 of the RedHat Enterprise Linux Figure 2 Simulated 12 C/ 13 C isotopologues of 20 carbon metabolite. In all 3 charts, the red bars represent the relative isotopologue intensities with contributions (contamination) from 13 C natural abundance. The green bars represent the corrected isotopologue intensities with this contamination removed. Chart A shows the expected isotopologues of a metabolite with no additional 13 C-labeling present. Chart B shows the expected isotopologues of a metabolite with equal amounts of 13 C-labeling for 8, 10, and 12 carbons. Chart C shows the expected isotopologues of a metabolite with 100% 13 C-labeling. Figure 1 Flowchart of 13 C natural abundance correction algorithm. Starting with isotopologue intensities and C Max given as input, the algorithm precalculates needed binomial coefficients. Next, the algorithm calculates the corrected intensities and uses them to calculate the natural abundance contaminated intensities. Then the algorithm compares the observed and calculated contaminated intensities and only continues for another iteration if an improvement is made.
operating system. The implementation is also very accurate. Given the perfect data in these three simulated test sets, the largest error was 4.12 × 10 -16 seen in the I M+1 corrected intensity for the test set representing no 13 Clabeling ( Figure 2A). Furthermore, the implementation appears quite robust since the relative error actually decreases as the number of carbons (C Max ) increases. At a C Max = 100, the relative error is 6.77 × 10 -17 . This implementation does have some numerical limitations, for example, the C Max must be less than 270 carbons due to all numerical quantities being represented as double precision (64 bit) floating point numbers. However, this limitation is easily overcome by using higher precision floating point numbers.
Application to phosphatidylcholine 34:1 observed isotopologue intensities Figure 3 shows the two sets of 12 C/ 13 C isotopologue intensities for phosphatidylcholine 34:1 (34 carbons in 2 fatty acid chains with only 1 double bond), with 13 C natural abundance contaminated intensities in red and corrected intensities in green. The algorithm converged within 8 iterations to produce the corrected intensity results. In comparing the contaminated and corrected intensities, the most significant changes are seen in isotopologues 0-4 and 16-20. The drastic drop in I M+1 , I M+2 , and I M+4 isotopologues make the incorporation of 13 C-labeled glycerol much clearer. Also, the drop in I M+16 , I M+18 , and I M+20 isotopologues supports the expected incorporation of 13 C-labeled acetyl groups in the fatty acid chain biosynthesis.

Discussion and Conclusions
Overall, correcting for the effects of natural abundance makes interpretation of isotopologue intensities from stable isotope tracing experiments easier within the context of cellular metabolism. Such a correction is required before using more quantitative methods of interpretation. Since the relative error is virtually zero with perfectly simulated data and the algorithm is very robust with increasing C Max , the accuracy of this correction is really only limited by the error in the isotopologue intensities themselves. Thus, the propagation of data error through this algorithm should be straightforward to analyze and quantify. Moreover, from Equation 5 it is evident that effects from natural abundance significantly link together groups of observed isotopologue intensities. This difference between calculated and observed intensities should be highly sensitive to the error in a set of isotopologue intensities. Therefore, this difference should be usable as an independent check on the quality of the observed set of isotopologue intensities. Such a quality control check would be especially useful when it is not possible or practical to repeat experiments or to determine whether additional experiments are necessary.

Cell Culture and FT-ICR-MS
We separated glycerophospholipids from crude cell extracts derived from MCF7-LCC2 cells in tissue culture after 24 hours of labeling with uniformly labeled 13 Cglucose. We analyzed the sample on a hybrid linear ion trap 7T FT-ICR mass spectrometer (Finnigan LTQ FT, Thermo Electron, Bremen, Germany) equipped with a TriVersa NanoMate ion source (Advion BioSciences, Ithaca, NY) as described elsewhere [6].